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))

Answer 1

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)

Answer 2

I found the answer by user2856 most helpful, but struggled with rounding errors if target resolution was not a multiple of source resolution, so rasterio.warp.reproject ended up being the better solution for me. Not needed for the problem of the OP, but posting here b/c it may be useful to some.

The rasterio.warp.reproject solution also offers to produce "target aligned pixels" (-tap option in the gdal commands)!

import rasterio
from rasterio.enums import Resampling
target_res = (2048, 2048)
with rasterio.open("myfile.tif") as src:
     data, transform = rasterio.warp.reproject(source=src.read(),
                                src_transform=src.transform,
                                src_crs=src.crs,
                                dst_crs=src.crs,
                                dst_nodata=src.nodata,
                                dst_resolution=target_res,
                                resampling=Resampling.bilinear)
     profile = src.profile
     profile.update(transform=transform, driver='GTiff',
                    height=data.shape[1], width=data.shape[2])
 with rasterio.open("resampled.tif", 'w', **profile) as dst:
                dst.write(data)

or, to get "target aligned pixels":

with rasterio.open("myfile.tif") as src:
 transform, width, height = rasterio.warp.aligned_target(
                               src.transform,
                               src.width,
                               src.height,
                               target_res)

 data, transform = rasterio.warp.reproject(source=src.read(),
                               destination=np.zeros((src.count, height, width)),
                               src_transform=src.transform,
                               dst_transform=transform,
                               src_crs=src.crs,
                               dst_crs=src.crs,
                               dst_nodata=src.nodata,
                               resampling=Resampling.bilinear)
 profile = src.profile
 profile.update(transform=transform, driver='GTiff',
                height=data.shape[1], width=data.shape[2])

 with rasterio.open("resampled_tap.tif", 'w', **profile) as dst:
                dst.write(data)