On Friday 17 July 2015 11:29:11 Mike Toews wrote: > On 17 July 2015 at 11:13, Even Rouault <[email protected]> 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 _______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
