7

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? Downsampled image along with the original image in QGIS

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

2 Answers 2

8

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)
4
  • Oh yes, you are totally right. I was downsampling instead of upsampling so I thought I would reverse the operator to match the documentation example. But I did not change the corresponding transform. I guess I just saw the division symbols and thought they should be symmetric. I have thanked you three times today and this is the fourth time. I am definitely going to see if I can get your example incorporated into the documentation.
    – krishnab
    Jul 25, 2019 at 4:30
  • @krishnab or just pass in a fraction for the scale parameter.
    – user2856
    Jul 25, 2019 at 4:33
  • Yep, I will just keep your code as is. P.S. I am in the process of opening the issue on github for rasterio with the updated documentation idea.
    – krishnab
    Jul 25, 2019 at 4:35
  • I added a new issue to the rasterio github page based on our conversations today. github.com/mapbox/rasterio/issues/1732 . I suggested the gist that you created as a better example. If the maintainers agree, then I can create the pull request for the documentation.
    – krishnab
    Jul 25, 2019 at 5:15
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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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