Create raster aligned with one created by gdal_translate

Denis Rykov
 

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.

Join main@rasterio.groups.io to automatically receive all group messages.