2

I used rasterio to scale a raster I have to a new resolution. It was at (5.0, 5.0) resolution, and I wanted to make it (15.0, 15.0).

However, the resulting resolution is not square anymore.

Here's my code (it's pretty standard, pretty much just from the rasterio documentation):

with rio.open(f_path) as ds:


        orig_res = ds.res[0] # 5.0
        des_res = 15.0


        # Calculate scale factor for this specific input layer
        upscale_factor = orig_res / des_res
        print("upscale factor: ", upscale_factor)

        # Scale while reading in array
        scaled_arr = ds.read(
            out_shape=(
                ds.count,
                int(ds.height * upscale_factor),
                int(ds.width * upscale_factor)                
            ),
            resampling=Resampling.bilinear
        )

        # Also have to update the transform to be scaled
        scaled_transform = ds.transform * ds.transform.scale(
            (ds.height / scaled_arr.shape[1]),
            (ds.width / scaled_arr.shape[2])
        )

        scaled_profile = ds.profile
        scaled_profile.update(transform=scaled_transform, height=scaled_arr.shape[1], width=scaled_arr.shape[2])



        with rio.open(f_out_path, 'w', **scaled_profile) as scaled_ds:
            scaled_ds.write(scaled_arr)

            print(f"Orig arr shape: ({ds.count}, {ds.shape[0]}, {ds.shape[1]})")

            print(f"Scaled ds shape: ({scaled_ds.count}, {scaled_ds.shape[0]}, {scaled_ds.shape[1]})")


            print("Orig ds res: ", ds.res)
            print("Scaled ds res: ", scaled_ds.res)

Here is the output printed:

upscale factor:  0.3333333333333333
Orig arr shape: (1, 9240, 2671)
Scaled ds shape: (1, 3080, 890)
Orig ds res:  (5.0, 5.0)
Scaled ds res:  (15.0, 15.00561797752809)

So is it bad that my new raster does not have square pixels anymore? (15.0, 15.00561797752809)

Will this affect anything with how I use the raster, or will it be okay? I also have other rasters that I'm up/down sampling to all be (15.0, 15.0), and they are also not exactly 15.0, so will it cause issues if my rasters have very slightly different resolutions?

Or is there a way to remedy this?

2
  • Try going the other way - multiply the resolution by a scale factor (i.e. 3 in your case) and divide the dimensions - see gist.github.com/lpinner/13244b5c589cda4fbdfa89b30a44005b
    – user2856
    Mar 13, 2020 at 4:20
  • @user2856 that worked perfectly! If you wanted to post that as answer with a short reason for why that way is better, that would be great!
    – rasen58
    Mar 13, 2020 at 18:23

1 Answer 1

1

I'm not sure why you're getting different cell sizes, but your issue is likely related to floating point error. Your upscale factor is 1/3 which can't be represented exactly in either binary (0.010101...) or decimal (0.333...).

I suggest going the other way - multiply the resolution by a scale factor (i.e. 3 in your case) so the cell sizes are multiples of the original and divide the dimensions.

with rio.open(f_path) as ds:


        orig_res = ds.res[0] # 5.0
        des_res = 15.0


        # Calculate scale factor for this specific input layer
        upscale_factor = des_res / orig_res  # 3
        print("upscale factor: ", upscale_factor)

        # Scale while reading in array
        scaled_arr = ds.read(
            out_shape=(
                ds.count,
                int(ds.height / upscale_factor),
                int(ds.width / upscale_factor)                
            ),
            resampling=Resampling.bilinear
        )

        # Also have to update the transform to be scaled
        scaled_transform = ds.transform * ds.transform.scale(
            (ds.height / scaled_arr.shape[1]),
            (ds.width / scaled_arr.shape[2])
        )


        scaled_profile = ds.profile
        scaled_profile.update(transform=scaled_transform, height=scaled_arr.shape[1], width=scaled_arr.shape[2])



        with rio.open(out_path, 'w', **scaled_profile) as scaled_ds:
            scaled_ds.write(scaled_arr)

            print(f"Orig arr shape: ({ds.count}, {ds.shape[0]}, {ds.shape[1]})")

            print(f"Scaled ds shape: ({scaled_ds.count}, {scaled_ds.shape[0]}, {scaled_ds.shape[1]})")


            print("Orig ds res: ", ds.res)
            print("Scaled ds res: ", scaled_ds.res)

And here is another example from a GitHub Gist (Apache 2.0 license)

from contextlib import contextmanager

import rasterio
from rasterio import Affine, MemoryFile
from rasterio.enums import Resampling


@contextmanager
def resample_raster(raster, out_path=None, scale=2):
    """ Resample a raster
        multiply the pixel size by the scale factor
        divide the dimensions by the scale factor
        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 500m and dimensions of (512, 512)
        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 125m and dimensions of (2048, 2048)
        returns a DatasetReader instance from either a filesystem raster or MemoryFile (if out_path is None)
    """
    t = raster.transform

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

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

    data = raster.read(
            out_shape=(raster.count, height, width),
            resampling=Resampling.bilinear,
        )

    if out_path is None:
        with write_mem_raster(data, **profile) as dataset:
            del data
            yield dataset

    else:
        with write_raster(out_path, data, **profile) as dataset:
            del data
            yield dataset


@contextmanager
def write_mem_raster(data, **profile):
    with MemoryFile() as memfile:
        with memfile.open(**profile) as dataset:  # Open as DatasetWriter
            dataset.write(data)

        with memfile.open() as dataset:  # Reopen as DatasetReader
            yield dataset  # Note yield not return


@contextmanager
def write_raster(path, data, **profile):

    with rasterio.open(path, 'w', **profile) as dataset:  # Open as DatasetWriter
        dataset.write(data)

    with rasterio.open(path) as dataset:  # Reopen as DatasetReader
        yield dataset

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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