3

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:

Raster A and B

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.

rasterio blocks structure

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:

wrong output written "vertically"

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:

raster A full extent

zoomed to see the stripes:

enter image description here

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)
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)
9
  • Thanks for your answer! I get an error when suing blockysize and blockxsize inside the bounds method: KeyError('blockysize'). As this is the key of the source dataset and it's not defined in your code, I thought it would be set by default to 256 but apparently it's not.
    – umbe1987
    Jan 31 '20 at 15:18
  • Solved using rasterio.readthedocs.io/en/latest/topics/profiles.html (kwds = src.profile, kwds['blockxsize'] = 256, kwds["blockysize"]). However, I get another error: TypeError('from_bounds() argument after ** must be a mapping, not tuple')
    – umbe1987
    Jan 31 '20 at 15:28
  • Got it: it has to be *src_bounds not **src_bounds. Trying your code, will let you know asap about the result.
    – umbe1987
    Jan 31 '20 at 15:33
  • Your code seems to work fine except that 1) the resulting raster has strips (as in EDIT part of my answer), and 2) pixels of input and output raster are shifted by exactly one row and column. Still much better than my code, hopefully I will achieve good results adjusting it (or maybe you have some ideas?)
    – umbe1987
    Jan 31 '20 at 15:59
  • 1
    Regarding the shift, I am not sure why this happens, but it might be a floating point issue? Try logging the bounds and dst_window. Your final window might have an floating point width and height which could result the shift? You could try to create a window with rounded values based on the results you get from dst_window
    – Thomas
    Jan 31 '20 at 16:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.