[gdal-dev] Sum pixel values in selected area of an image

[Lukasz Tracewski]

Unless I did not correctly understand the advice or the 'gdalwarp' manual page, the first point ("make sure both raster have same dimensions") is deviating from what I want to achieve here. Now I see that more detailed information might be needed - my apologies for the confusion.
Raster with forest are huge images, with 30 meters spatial resolution, and size of 36000 x 36000 pixels (lets call them "forest rasters"). The other raster describes area over which I would like to calculate the forest area. Those are bilevel images, with "1" meaning that the pixel should be included in calculations and "0" otherwise. Physically they represent habitats of certain animal species and spawn from a couple of hundreds km2 to hundreds of thousands km2. We will call them "animal rasters".
Here is an example how such habitat for a certain species looks like in ASCII format:
ncols 6nrows 5xllcorner 916000yllcorner 496000cellsize 1000.000000NODATA_value -99990 0 0 0 0 0 0 1 1 1 1 0 1 1 1 1 1 0 0 0 1 1 0 0 0 0 0 0 0 0 
The scale is 1000 meters. The forest raster is 30 meters.
If I now follow the proposed workflow, I have to rescale 36,000 x 36,000 pixels forest raster to 6x5 pixels, loosing therefore all resolution. That is of course assuming I understood correctly what Even proposed.
What I would like to get is something as follows:1. With gdalinfo read pixel size of "animal raster", in this case 0.008950249595932.2. Rescale the "forest raster" to match the scale: gdalwarp -tr 0.008950249595932 0.008950249595932 -r average -tap forest.tif forest_scaled.tif . Now they are at the same scale, but obviously the dimensions are different. Purpose of this part is to make the calculations faster.3. Multiply the two rasters, taking into account dimensionality. Since both images are goereferenced, I would expect that there is something in GDAL that takes that into account. If they are at the same scale, one could superimpose them and take nearest neighbour for multiplication.
I do not mind at all going to C++ or Python to achieve this. 

Thanks for help!Lucas

> From: even.rouault at spatialys.com
> To: gdal-dev at lists.osgeo.org
> Subject: Re: [gdal-dev] Sum pixel values in selected area of an image
> Date: Thu, 2 Oct 2014 17:00:46 +0200
> CC: lukasz.tracewski at outlook.com
> 
> Le jeudi 02 octobre 2014 09:02:00, Lukasz Tracewski a écrit :
> > Hi,
> > I only recently started my adventure with GDAL and GIS in general, so
> > please accept my apologies for perhaps maybe the most accurate formulation
> > of the problem. I also did my best to find the answer already on this
> > mailing list and outside, but to no avail - again possibly due to lack of
> > experience. Last year Hansen et al. prepared a detailed map, 30 meters
> > resolution, of forest cover for the whole planet. The data is publicly
> > accessible, both for online viewing and download. Both can be found here:
> > http://earthenginepartners.appspot.com/science-2013-global-forest The
> > forest cover images are GeoTIFF and have pixel values in range [0, 100]
> > that describes percentage of forest cover. I have some georeferenced
> > images that essentially are composed of ones and zeros, e.g.:0 1 1 0 00 1
> > 1 1 00 1 0 0 0... Those images are in 1000 meters scale and can spawn over
> > whole continents.  My aim is to calculate forest cover for them: wherever
> > the value is "1", this pixel should be added to the forest cover. Say the
> > Hansen image looks like this:10 20 30 40 5010 80 80 0 00 0 0 0 0... Then
> > rows and column of this images should be multiplied by respective rows and
> > columns of my own image, producing: 0 * 10 + 20 * 1 + 30 * 1 + 40 * 0 + 50
> > * 0 + 10 * 0 + 80 * 1 + ... = 210 Mind that scale of both images is
> > different, so the example above is actually not accurate. The real
> > calculations should average pixel values, then sum them and then convert
> > to area. Can anyone point me where to start? Maybe you know of any
> > examples that could give me a hint?
> 
> Lukasz,
> 
> There are likely many ways of doing that. Here's a potential workflow
> 
> 1) With gdalwarp -r averge -ts you could make sure both raster have same 
> dimensions
> 2) Use gdal_calc.py to do the multiplication pixel by pixel.
> 3) Use gdalinfo -stats to compute statistics on the raster generated by 2). 
> Multiply the MEAN value by the raster width and raster height, and that will 
> give you the sum.
> 
> Even
> 
> > 
> > Thanks,Lucas
> 
> -- 
> Spatialys - Geospatial professional services
> http://www.spatialys.com
 		 	   		  
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20141005/948f5dd0/attachment.html>