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

I am not 100 sure, but I now simply rounded all windows (origin + offsets) and this seems to have fixed it.

Thank you for the info. I saw shorter values tan in QGIS and those came from Python's float limitations. When using the Decimal class I also get the same over-precise values.

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.

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

When opening my raster with QGIS, I can see the Pixel Size of my raster as:

0.0008333333299974727193,-0.000833333330020215703

When I read he same raster with rasterio, I get pixel size values with a lower precision and I guess this causes some problems when calculating new affines for another extent

0.0008333333299974727 0.0008333333300202157

I have tries to pass the decimal_precision=X when opening my raster with rasterio, but this has no effect. I guess this is a general float problem on my system. My suggestion therefore would be to use python Decimals or any other type representing real floats.

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.

Hi! I had faced the similar shifting issue yesterday: https://github.com/mapbox/rasterio/issues/1932

• include the source raster (within)
• lay within the source raster bounds (contained)
• or partially overlap with the source raster bounds

1. I have calculated the exact pixel aligned bounds for the target raster
2. I have the exact pixel aligned bounds from the source raster
3. I create an intersection polygon from both
   intersection = INTERSECT(source.bounds, target.bounds)
4. I use the intersection to create a read-window for the source
   with rasterio.open('source.tif', 'r') as source:
       read_window = source.window(*intersection.bounds))
5. I use the intersection to create a write-window for the target
   with rasterio.open('target.tif', 'w', **profile) as target:
       write_window = target.window(*intersection.bounds))
6. Then I read from the source and write to the target
   target.write(source.read(window=read_window), window=write_window)