Should resampling/downsampling a raster (using Rasterio) cause the coordinates to change?

Question

I was able to successfully downsample a raster file, after getting Stack Exchange help. I loaded the new downsampled file into QGIS along with the original file, to check that things went well. The tricky thing was that the new image seems a bit shifted and smaller. So it seems like the coordinates in the image changed along with its size.

I was not sure what to expect after the operation, but I imagined that the two files should have the same coordinates--because the CRS was preserved in the downsampling. I figured that the two images would be the same size, but one would be blurrier than the other because its resolution was downsampled.

Is this image below the way that things should look after downsampling?

Here is the code I used to resize.

from contextlib import contextmanager  
import rasterio
from rasterio import Affine
from rasterio.enums import Resampling

dat = 'original_image.tif'

@contextmanager
def resample_raster(raster, scale=2):
    t = raster.transform

# rescale the metadata
    transform = Affine(t.a / scale, t.b, t.c, t.d, t.e / scale, t.f)
    height = raster.height / scale
    width = raster.width / scale

profile = src.profile
    profile.update(transform=transform, driver='GTiff', height=height, width=width, crs=src.crs)

data = raster.read( # Note changed order of indexes, arrays are band, row, col order not row, col, band
            out_shape=(int(raster.count), int(height), int(width)),
            resampling=Resampling.cubic)

with rasterio.open('image-15cm.tif','w', **profile) as dst:
        dst.write(data)
        yield data

with rasterio.open(dat) as src:
    with resample_raster(src, 3.5) as resampled:
        print('Orig dims: {}, New dims: {}'.format(src.shape, resampled.shape))
        print(repr(resampled))

user2856 · Accepted Answer

You changed the following lines from my original code from:

transform = Affine(t.a / scale, t.b, t.c, t.d, t.e / scale, t.f) # <== division
height = raster.height * scale                                   # <== multiplication
width = raster.width * scale

To:

transform = Affine(t.a / scale, t.b, t.c, t.d, t.e / scale, t.f)   # <== division
height = raster.height / scale                                     # <== division again!!!
width = raster.width / scale

But did not change the transform to match.

My code uses the scale factor to multiply the dimensions of a raster, i.e. given a pixel size of 250m, dimensions of (1024, 1024) and a scale of 2, the resampled raster would have an output pixel size of 125m and dimensions of (2048, 2048).

If you want to use my original code as is to decrease the resolution of a raster - pass in a fraction as the scale, i.e given a pixel size of 250m, dimensions of (1024, 1024) and a scale of 0.5, the resampled raster would have an output pixel size of 500m and dimensions of (512, 512)..  Note you will have to change height and width to ensure integer output.

height = int(raster.height * scale)
width = int(raster.width * scale)

If you want to use the scale to divide the resolution of a raster, you then need to also change the operator you use to modify the transform

i.e.

# rescale the metadata
transform = Affine(t.a * scale, t.b, t.c, t.d, t.e * scale, t.f)  # <== multiplication
height = int(raster.height // scale)                              # <== division
width = int(raster.width // scale)

Should resampling/downsampling a raster (using Rasterio) cause the coordinates to change?

One Answer

Add your own answers!

Ask a Question