7

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
5
  • Perhaps gdalwarp with our favorite resampling method gdal.org/gdalwarp.html? If the partial pixels make problem oversample to 10 m pixel size first. Commented Nov 18, 2017 at 20:50
  • there is no way to calculate the sum? Commented May 13, 2019 at 14:38
  • Unfortunately, No Commented May 13, 2019 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 Commented May 13, 2019 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 Commented May 13, 2019 at 17:25

1 Answer 1

Reset to default
8

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
13
  • 1
    GDAL 2.x brought in a few more -r resampling method options, such as min, max, med, q1, q3; see details here. Commented Nov 18, 2017 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 )? Commented Nov 18, 2017 at 21:50
  • 1
    @AntonioFalciano the extra resampling methods are available for Python too, just a bit hidden. Commented Nov 19, 2017 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. Commented Nov 20, 2017 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. Commented Nov 20, 2017 at 21:16

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.