Search code examples
pythongdalresamplinggeotiff

How to Resample GeoTIFF on Python?

I am trying to resample a GeoTIFF file from pixels representing 30m² to pixels representing 1000m² on Python.

I have a geotiff that has the data according to those coordinates. The spatial resolution is 30m². I would like to change that to 1000m².

I was originally able to do this on QGIS while projecting (warping) the CRS from W84 to EPSG 7390. However, I would like to do the resampling without changing the CRS. When I tried it on QGIS I got an error ("creating 0x0 dataset).

I tried doing it on Python using GDAL but ran into the same error. Here is the code I used (from this Stack Exchange post:

from osgeo import gdal

infn = '/path/to/source.tif'
outfn = '/path/to/target.tif'

xres=1000
yres=1000
resample_alg = 'near'

ds = gdal.Warp(outfn, infn, xRes=xres, yRes=yres, resampleAlg=resample_alg)
ds = None

This is the GDAL error:

ERROR 1: Attempt to create 0x0 dataset is illegal,sizes must be larger than zero.

Also, just to check that the problem wasn't in opening the files, I also ran:

referenceFile = infn
reference = gdal.Open(referenceFile, 0)  # this opens the file in only reading mode
referenceTrans = reference.GetGeoTransform()

To which the output was:

(-58.168519761062974, 0.0002694945852358564, 0.0, -17.16626609035358, 0.0, -0.0002694945852358564)

I'm not entirely sure what I'm doing wrong here. I am very new to geographical data and would really appreciate some help. How do I

Note: I wasn't sure whether to post this here or on StackExchange, so please let me know if it doesn't belong here.

Solution

  • You need to provide the resolution in the unit of the output projection. So if you don't change it, and the projection is EPSG:4326 (WGS84), the units are degrees. So you are asking GDAL to warp to a resolution of 1000 degrees, which as the error shows results in an invalid output dimension of 0x0 pixels.

    A resolution of 1000 meters is more or less ~0.009 degrees (1km / 111km per degree latitude). There isn't a 1-on-1 mapping between resolution in meters and degrees, because it depends on the location itself. You can see the table in the link below for some examples:
    https://en.wikipedia.org/wiki/Longitude#Length_of_a_degree_of_longitude

    So using 0.009° as the resolution should be closer to what you expect, but the exact resolution you need might depend on the actual application.