Is having rectangular cell sizes in a raster bad (after resampling)?

Question

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?

Answer 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