Hi everyone!

I tried to convert gdal_translate command to rasterio and faced some problem I described here. The original command is:

gdal_translate -oo OVERVIEW_LEVEL=1 -projwin 10816494.68930465 3653956.75283452 10819454.39626612 3650628.7941839 /vsicurl/https://raw.githubusercontent.com/satellogic/telluric/master/tests/data/google.xml gdalized.tif

Eventually I came up with the following solution that produces raster aligned with one produced by gdal_translate:

import rasterio
from rasterio.windows import Window
from math import floor
with rasterio.open(
    "/vsicurl/https://raw.githubusercontent.com/satellogic/telluric/master/tests/data/google.xml",
    OVERVIEW_LEVEL=1,
) as src:
    bounds = (10816494.68930465, 3650628.7941839, 10819454.39626612, 3653956.75283452)
    xoff = floor((bounds[0] - src.transform.c) / src.transform.a) - 1
    yoff = floor((bounds[3] - src.transform.f) / src.transform.e) - 1
    xsize = int((bounds[2] - bounds[0]) / src.transform.a + 0.5)
    ysize = int((bounds[1] - bounds[3]) / src.transform.e + 0.5)
    out_window = Window(xoff, yoff, xsize, ysize)
    height = int(out_window.height)
    width = int(out_window.width)
    out_kwargs = src.profile
    out_kwargs.update(
        {
            "driver": "GTiff",
            "height": height,
            "width": width,
            "transform": src.window_transform(out_window),
        }
    )
    with rasterio.open("aligned.tif", "w", **out_kwargs) as out:
        out.write(src.read(window=out_window, out_shape=(src.count, height, width),))

It works fine but I still have couple of questions:

• Why I have to subtract 1 to get right values of xoff and yoff?
• How can I adjust spatial extent of raster to match extent of raster produced by gdal_translate? Now there is 1 px difference between them (but it is definitely not shift because all overlapped pixels have the same value).

Any advice would be appreciated.

if window is not None and window.is_polygon:
                    # Make sure we operate with the same coordinates
                    window = window.transform(profile['crs'].to_epsg())

                    # Get the coordinates from raster and window bounds
                    s_xmin, s_ymin, s_xmax, s_ymax = raster_bounds.bounds
                    w_xmin, w_ymin, w_xmax, w_ymax = window.bounds

                    # Get the source pixel width and height
                    x_pixel = profile['transform'].a
                    y_pixel = profile['transform'].e

                    # Calculate the shift vector from source bounds to target bounds
                    shift_xmin = w_xmin - s_xmin
                    shift_ymin = w_ymin - s_ymin
                    shift_xmax = w_xmax - s_xmax
                    shift_ymax = w_ymax - s_ymax

                    # Normalize the shift vector by the source pixel width and height
                    shift_xmin_norm = ceil(shift_xmin / x_pixel) * x_pixel
                    shift_ymin_norm = floor(shift_ymin / y_pixel) * y_pixel
                    shift_xmax_norm = floor(shift_xmax / x_pixel) * x_pixel
                    shift_ymax_norm = ceil(shift_ymax / y_pixel) * y_pixel

                    # Create the new boundaries of the target raster
                    new_xmin = s_xmin + shift_xmin_norm
                    new_ymin = s_ymin + shift_ymin_norm
                    new_xmax = s_xmax + shift_xmax_norm
                    new_ymax = s_ymax + shift_ymax_norm

                    # Calculate width & height (Pixel) for the target raster
                    new_width = abs(int((new_xmax - new_xmin) / x_pixel))
                    new_height = abs(int((new_ymax - new_ymin) / y_pixel))
                    if new_width == 0 or new_height == 0:
                        raise Exception(f'New raster has invalid measure (width or height is 0 pixels)')

                    # Update profile with windows dimensions
                    profile.update(width=new_width)
                    profile.update(height=new_height)

                    # Create new affine and update profile with it
                    new_affine = rasterio.transform.from_origin(
                        new_xmin,
                        new_ymax,
                        new_width,
                        new_height
                    )
                    profile.update(transform=new_affine)

But i have replaced all this now by:

# Make sure we operate with the same coordinates
    # Make sure we operate with the same coordinates
    window = window.transform(profile['crs'].to_epsg())

    # Create a new window and affine from the window extent
    raster_window = source.window(*window.bounds)
    raster_window = raster_window.round_offsets(op='ceil').round_lengths(op='floor')

    # Update target raster with the new affine
    new_affine = source.window_transform(raster_window)
    profile.update(transform=new_affine)

    # Update target raster with new windows shape (Make sure raster stays within extent)
    c_xmin, c_ymin, c_xmax, c_ymax = source.window_bounds(raster_window)
    w_xmin, w_ymin, w_xmax, w_ymax = window.bounds
    x_size, y_size = source.res
    corrected_width = floor(abs(c_xmin - w_xmax) / x_size)
    corrected_height = floor(abs(w_ymin - c_ymax) / y_size)
    profile.update(width=corrected_width)
    profile.update(height=corrected_height)

Both methods actually created and affine with the same values, that is why I replaced my own code with the rasterio functions.

Thanks Thomas. I've tried your approach but it still gives me a shifted output compared to gdal_translate.

I also still fight with the shifted by one pixel problem for 1 extent out of 5 extents. 4 return perfectly aligned clipped subsets of the source raster, 1 does not.

The issue has been fixed via https://github.com/mapbox/rasterio/commit/15d8dd759256f12d412433f5b9824c8d16b33f8d.