5

I have a GeoTIFF which is at 30m resolution. I need to downsample it to 250m resolution (using some aggregating function such as mean, mode, sum). My stack is python (rasterio and GDAL bindings)

I can find a way to do this directly using numpy writing some nasty double loops and iterating over the matrices directly but looking for a cleaner way.

Note that the main issue is that 250 is not a multiple of 30 and there will be some overflow. For instance if the aggregating function is the sum the algorithm would be to use:

  1. 8 X 8 block which is 240m X 240m in base resolution
  2. Chose a partial pixel strip at 10m X 10m to complete the sum to 250m
  3. The value for 2. would be area weighted and added to the sum from 1.

Please let me know if I was unclear. Again I can find a brute force way to do it using matrix manipulations but was wondering if there is an elegant way. I am open to using R also and read about the projectRaster and aggregate methods but they seem to work on integer cut offs also.

In general it may be really cumbersome to do this for an arbitrary resolution say 31.567m or so.

|improve this question
  • Perhaps gdalwarp with our favorite resampling method gdal.org/gdalwarp.html? If the partial pixels make problem oversample to 10 m pixel size first. – user30184 Nov 18 '17 at 20:50
  • there is no way to calculate the sum? – Digd May 13 at 14:38
  • Unfortunately, No – nickves May 13 at 14:47
  • This does not provide an answer to the question. To critique or request clarification from an author, leave a comment below their post. - From Review – Vince May 13 at 16:47
  • If you have a new question, please ask it by clicking the Ask Question button. Include a link to this question if it helps provide context. - From Review – whyzar May 13 at 17:25
4

I'd use gdalwarp to resample (downsample in this case) data:

gdalwarp -tr 250 250 -r average input_r30.tif output_r250.tif

Instead, using the GDAL Python bindings:

from osgeo import gdal

inDs = gdal.Open('input_r30.tif')
outDs = gdal.Warp('output_r250.tif', inDs,
                  format = 'GTiff',
                  xRes = 250, yRes = 250,
                  resampleAlg = gdal.GRA_Average)

inDs = None
outDs = None

In both examples, I've used the average resampling method, but you can choose it among the supported ones:

GRA_NearestNeighbour = 0
GRA_Bilinear = 1
GRA_Cubic = 2
GRA_CubicSpline = 3
GRA_Lanczos = 4
GRA_Average = 5
GRA_Mode = 6
GRA_Gauss = 7

Because 250 is not a multiple of 30, the output will have a little bigger spatial extent than the input.

UPDATE: Thanks to @MikeT, who has opened a ticket in GDAL bug tracker, "the extra GRA constants will be available for convenience through the SWIG Python bindings in the next GDAL 2.2.3 release, soon to be available." These are the relatively new constants:

GRA_Max = 8
GRA_Min = 9
GRA_Med = 10
GRA_Q1 = 11
GRA_Q3 = 12

UPDATE: GDAL 2.2.3 has been released today (2017/11/23):
https://lists.osgeo.org/pipermail/gdal-dev/2017-November/047782.html

|improve this answer
  • 1
    GDAL 2.x brought in a few more -r resampling method options, such as min, max, med, q1, q3; see details here. – Mike T Nov 18 '17 at 21:14
  • 1
    Thanks a lot Antonio. Few more questions 1) How does it treat partial pixels, in other words if the input resolution is say 35.576m, and the output is 100m, does it combine 3 pixels into 1 or fractional pixels? 2) Also why would be extent be larger? I also notice that it tends to give different pixel size if I don't constrain it using the -tr parameter, why is that, In other words if I dont specify constraints on the output how does it chose the pixel size while warping (say from aea proj to utm )? – ganesh reddy Nov 18 '17 at 21:50
  • 1
    @AntonioFalciano the extra resampling methods are available for Python too, just a bit hidden. – Mike T Nov 19 '17 at 20:33
  • 1
    @dr_times GRIORA should mean GDAL Raster IO Resample Algorithm, while GRA GDAL Resample Algorithm. The correspondent constants (e.g. GRIORA.Average and GRA.Average) should have the same value (e.g. 5) and so are equivalent operationally. – Antonio Falciano Nov 20 '17 at 15:25
  • 1
    @AntonioFalciano the extra GRA constants will be available for convenience through the SWIG Python bindings in the next GDAL 2.2.3 release, soon to be available. – Mike T Nov 20 '17 at 21:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.