Aggregating rasters with GDAL

Question

I have a set of files that have been mosaiced into a global VRT in a particular projection and resolution. However, I have a need to report statistics on a (say) one degree grid, where the reported statistics will be e.g. count, mean, mode or some percentile of the input pixels that fall within the grid cell. To make things more fun, there are NODATA instances too.

My solution is to create the mosaiced VRT in WGS84 Lat/Long, enforcing an output pixel size so that 1.0 % output_pixel_size == 0 (in other words, I can get an integer number of pixels into a 1 degree cell). Then, you run gdaladdo on the VRT and extract what you want. This should work for averages, but percentiles and some other things are a bit more complicated.

I could do it in Python, but I was wondering if there's an off-the-shelf solution I could use?

Answer 1

In the end, I have solved the problem in the following fashion:

1. I generate a global VRT in the source projection using gdalbuildvrt
2. I use gdalwarp to "reproject" the first mosaicked VRT into WGS84 Lat/Long, with pixel spacing being the closest to the original projection area at the equator.
3. I have written a script in Python that reads in a chunk of pixels under the grid cell that I want using the GDAL Python bindings.
4. The result of 3 is a numpy array, and you can do whatever you want with it.

The code is quite specific as I'm actually reading in to datasets to do some extra processing, but it's fairly simple code.

Answer 2

I'm not sure if you are using ArcGIS...? But this is relatively easy using the Mosaic Dataset structure - it automatically handles data in different projections and at different resolutions. Especially if this is not just a one-time project -- If you'll have to maintain and update this over time, the Mosaic Dataset makes it very easy to manage and generate a variety of outputs (e.g. resampled to 1 degree grid, or a wide range of functions that can calculate new products on-the-fly)...