Hi all,

I am trying to rewrite some legacy code that was using gdal_translate -projwin to use rasterio. The goal of the code is just to crop an image, given geographical bounds in the CRS of the image.

I struggled to understand what gdal_translate was doing when the projection of the projwin bounds into image coordinates did not align perfectly with the pixels, but I think I have it figured out now:

• the coordinates of the upper left corner of projwin are "floored" to match integer pixel values
• the pixel width and height derived from projwin are "rounded" to match integer values

(the projwin documentation doesn't really explain what happens in that case, I had to find through experiments).

I now have a python function doing nearly exactly the same thing as what gdal_translate -projwin does:

import subprocess

import rasterio


def rio_crop(src_path, dst_path, ulx, uly, lrx, lry):
    with rasterio.open(src_path) as src:
        bounds = (ulx, lry, lrx, uly)

        # Get the pixel coordinates of the bounds in src_path
        window = src.window(*bounds)

        # Do a 'ceil' operation on the offsets
        window = window.round_offsets()

        # Do a 'round' operation on the lengths
        # window.round_lengths() only offers to do a 'ceil' or 'floor', not a 'round'
        width = round(window.width)
        height = round(window.height)
        window = rasterio.windows.Window(window.col_off, window.row_off, width, height)

        profile = src.profile
        profile.update(
            {
                "driver": "GTiff",
                "compress": "deflate",
                "height": height,
                "width": width,
                "transform": src.window_transform(window),
            }
        )

        with rasterio.open(dst_path, "w", **profile) as out:
            out.write(src.read(window=window))


def gdal_crop(src_path, dst_path, ulx, uly, lrx, lry):
    cmd = [
        "gdal_translate",
        src_path,
        dst_path,
        "-of",
        "GTiff",
        "-projwin",
        str(ulx),
        str(uly),
        str(lrx),
        str(lry),
    ]
    subprocess.check_output(cmd)


if __name__ == "__main__":
    src_path = "/vsicurl/https://storage.googleapis.com/gcp-public-data-landsat/LC08/01/176/039/LC08_L1TP_176039_20180411_20180417_01_T1/LC08_L1TP_176039_20180411_20180417_01_T1_B8.TIF"
    ulx = 318753
    uly = 3319235
    lrx = 321313
    lry = 3316675
    for delta in range(0, 30, 6):
        print(delta)
        rio_crop(src_path, f"rio_{delta}.tif", ulx, uly, lrx + delta, lry + delta)
        gdal_crop(src_path, f"gdal_{delta}.tif", ulx, uly, lrx + delta, lry + delta)
        with rasterio.open(f"rio_{delta}.tif") as rio:
            with rasterio.open(f"gdal_{delta}.tif") as gdal:
                assert rio.bounds == gdal.bounds
                print(rio.tags()["AREA_OR_POINT"])
                print(gdal.tags()["AREA_OR_POINT"])

My questions are:

1. Why is there no way to do a 'round' (not a 'ceil' or a 'floor') with the window.round_lengths() function? If it existed, it would make my code simpler.
2. When I compare the images output by rio_crop and by gdal_crop, they differ by their value of AREA_OR_POINT. GDAL preserves the original value (POINT), whereas rasterio changes it to AREA. I am always confused by the meaning of this variable, but I was wondering why rasterio doesn't behave like GDAL in that case?

Thanks!

Guillaume Lostis

• Why is there no way to do a 'round' (not a 'ceil' or a 'floor') with the window.round_lengths() function? If it existed, it would make my code simpler.

I've wondered the same (and think even asked about it within another issue? though maybe I just thought about doing that). I'd be happy to add support for `round` if there's no opposition from Sean or other maintainers. 

Haha, then I didn't my due diligence either. :)