Unexpected result using rasterio windowed writing - Access window out of range in RasterIO()

Question

I am trying to write (copy and paste) a raster into another raster using raserio's blocks and windows.

I have two rasters (A and B):

A has a greater extent than B and it completely contains B. The two rasters share the same cell size (roughly 30m) and CRS (EPSG:4326).

The lower and right boundaries of A and B overlap like below:

I am having hard times understanding the problem of my code. I think the issue may be in the order I am iterating through each block when using block_windows() method, but I might be wrong.

Looking at the picture in the rasterio's documentation and printing the block indexes at each iteration, it seems the blocks should write east to west (let's say "horizontally"), and once the upper part is all written, the process should jump to the next lower row.

However, when I run my code, it seems the writing process happens "vertically" (north to south), and then jump to the next column). For this problem, when I reach the end of the first column, the code tries to write into a block that is outside A extent and I receive the error:

RasterioIOError('Read or write failed. /home/umberto/test/output/mosaic2.tif: Access window out of range in RasterIO().  Requested (42379,56159) of size 256x256 on raster of 85196x56160.')

Here is what I see in QGIS:

The black rectangle represents raster A where I am trying to paste B block by block. The whitish part in the middle of the raster represents the 256X256 blocks in the first column of B written to A.

Unfortunately, the writing process cannot jump to the next column as the the indexes of rows makes the code trying to access a part of the raster which is below the raster extent itself.

How am I doing it

The idea is to write each window of raster B into raster A. However, as the window object column and row offsets are different when considering the two rasters, I am translating for each block the offsets from B to A (window and new_window in the below code).

Here is my code (simplified but still relevant):

with rasterio.open(inputfile, 'r') as src:
    with rasterio.open(
        os.path.join(dest_dir, 'mosaic2.tif'),
        'w',
        driver='GTiff',
        compress="DEFLATE",
        height=56160,
        width=85196,
        count=1,
        dtype=np.float64,
        crs='+proj=latlong',
        transform=out_transform,
    ) as mosaic_raster:
        for ji, window in src.block_windows(1):
            r = src.read(1, window=window)
            left, top = src.xy(
                window.col_off, window.row_off, offset='ul')
            right, bottom = src.xy(
                window.col_off + window.width, window.row_off + window.height, offset='ul')
            new_window = from_bounds(
                left, bottom, right, top, mosaic_raster.transform, window.height, window.width)
            mosaic_raster.write(r, 1, window=Window(int(new_window.col_off), int(
                new_window.row_off), int(new_window.width), int(new_window.height)))

EDIT

I have made some interesting improvements:

I am now able to write the full extent of raster A correctly over raster B, having the same pixel values in the same locations for both rasters.

The only thing left before having a complete solution is that the output raster has some (apparently regularly spaced) stripes, sign that there might be some adjustments to do in the window coordinates translation.

Here is a picture of the updated situation:

zoomed to see the stripes:

Finally, the updated code:

with rasterio.open(inputfile, 'r') as src:
    with rasterio.open(
        os.path.join(dest_dir, 'mosaic2.tif'),
        'w',
        driver='GTiff',
        compress="DEFLATE",
        height=56160,
        width=85196,
        count=1,
        dtype=np.float64,
        crs='+proj=latlong',
        transform=out_transform,
    ) as mosaic_raster:
        for ji, window in src.block_windows(1):
            r = src.read(1, window=window)
            left, top = src.xy(
                row=window.row_off,
                col=window.col_off)
            right, bottom = src.xy(
                row=window.row_off + window.height,
                col=window.col_off + window.width)
            new_window = from_bounds(
                left=left,
                bottom=bottom,
                right=right,
                top=top,
                transform=mosaic_raster.transform,
                height=window.height,
                width=window.width)
            mosaic_raster.write(r, 1, window=new_window)

Answer 1

The window describes the relative location of the block within your file.

If CRS of both files are identical, you can use the bounds of the input window to compute its relative location within output window

import numpy as np
import rasterio
from rasterio import Affine
from rasterio.windows import bounds, from_bounds

# adopt as needed
input_file = "/some/file.tif"
output_file = "/another/file.tif"
out_transform = Affine(1, 0, 0, 0, -1, 0)

profile = {"driver": 'GTiff',
           "compress": "DEFLATE",
           "height": 56160,
           "width": 85196,
           "count": 1,
           "dtype": np.float64,
           "crs": '+proj=latlong',
           "transform": out_transform,
           "tiled": True}

with rasterio.open(input_file, 'r') as src:
    with rasterio.open(output_file, 'w', **profile) as dst:
        for _, src_window in src.block_windows(1):
            data = src.read(1, window=src_window)
            # convert relative input window location to relative output window location
            # using real world coordinates (bounds)
            src_bounds = bounds(src_window, transform=src.profile["transform"])
            dst_window = from_bounds(*src_bounds, transform=dst.profile["transform"])

            dst.write(data, window=dst_window, indexes=1)