[gdal-dev] Gdal_calc.py with more than 26 input files

[Even Rouault]

On Friday 17 July 2015 11:29:11 Mike Toews wrote:
> On 17 July 2015 at 11:13, Even Rouault <even.rouault at spatialys.com> wrote:
> > But when you have more than 26 input files, it is dubious that you really
> > want to individually identify them with a particular letter/name. As in
> > one of the examples, you likely just want to apply a global operation on
> > them, like sum, mean, .... So they would be better passed as a lis
> > t so as to be able to build a 3 dimension numpy array, so you can do
> > things
> > like:
> > 
> > --input foo.tif bar.tif ... --calc "rasters.mean(axis=0)"
> > 
> > And you could also still reference a particular raster with rasters[i]
> > syntax
> > 
> > (To be clear: the above is not something currently available)
> 
> It's actually not too far off. For instance, stack several single-band
> var*.tif rasters in a multi-band VRT, then calculate the mean along
> the 0th axis:
> 
> gdalbuildvrt -separate var_stacked.vrt var*.tif
> gdal_calc.py -A var_stacked.vrt --outfile=avg_var.tif
> --calc="A.mean(axis=0)"
> 
> However, the result avg_var.tif is not as expected due to the use of
> blocks. Is this ticket worthy?

In your above example, only the first band of var_stacked.vrt will be used. So 
basically you're computing the average of all lines of a block (possibly a 
single line depending on the dataset), and this average line will be written 
in the first line of each block.

Perhaps there should be a check on the shape of the result to verify it 
matches the shape of the input blocks and error out if not, if that's what you 
meant ?
Unless there would be some new option to mean that a multiband raster should 
be loaded as a 3d numpy array. The current behaviour of using a single band is 
a documented one, so we cannot change that without breaking existing valid 
uses.


> -Mike

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com