Using the GridR’s Grid Resampling Chain - Basic

This guide demonstrates how to effectively use GridR’s basic_grid_resampling_chain that wraps the core array_grid_resampling method. It covers key aspects of its operation, helping you understand how to manage your data and resources efficiently.

Here’s what we’ll explore:

  • Managing I/O Datasets: How to properly handle input and output datasets when working with the chaining function.

  • Memory Resource Management: Techniques for controlling memory usage during computations.

  • Extended Computational Features: A comparison of the chain’s capabilities versus the core array_grid_resampling method.

First, let’s address “why basic”? This chain is prefixed “basic” because it provides users with direct control over several memory usage-related parameters, allowing for fine-tuned resource management. Currently, this memory management is not automatic, requiring users to adapt these parameters for different use cases; automatic management is identified as a future enhancement for improved usability.

Setting things up

First proceed with some import : standard, community tiers and grid packages.

We also import the well known mandrill raster image.

from gridr.misc.mandrill import mandrill

In order to work with the basic_grid_resampling_chain we will need at least :

  • a raster image to read

  • a resampling grid to read

First we are going to generate the raster image from our 3 channels mandrill image.

RASTER_IN = './grid_resampling_chain_001_raster_in.tif'
GRID_IN_F64 = './grid_resampling_chain_001_grid_in_f64.tif'

Before we dive into the main topic of this guide, we’ll also need to define some utility functions. These will help us write our data and create a regular grid.

import numpy as np
import rasterio
import warnings

warnings.filterwarnings("ignore", category=rasterio.errors.NotGeoreferencedWarning)

def write_array(array, dtype, fileout):
    if array.ndim == 3:
        with rasterio.open(fileout, "w", driver="GTiff", dtype=dtype,
                height=array.shape[1], width=array.shape[2], count=array.shape[0],
                ) as array_in_ds:
            for band in range(array.shape[0]):
                array_in_ds.write(array[band].astype(dtype), band+1)
            array_in_ds = None
    elif array.ndim == 2:
        with rasterio.open(fileout, "w", driver="GTiff", dtype=dtype,
                height=array.shape[0], width=array.shape[1], count=1,
                ) as array_in_ds:
            array_in_ds.write(array.astype(dtype), 1)
            array_in_ds = None

def shape2(array):
    if array.ndim == 3:
        return (array.shape[1], array.shape[2])
    else:
        return array.shape

# write mandrill as tif
write_array(mandrill, dtype=mandrill.dtype, fileout=RASTER_IN)
def create_grid(nrow, ncol, origin_pos, origin_node, v_row_y, v_row_x, v_col_y, v_col_x, dtype):
    """
    """
    x = np.arange(0, ncol, dtype=dtype)
    y = np.arange(0, nrow, dtype=dtype)
    xx, yy = np.meshgrid(x, y)
    xx -= origin_pos[0]
    yy -= origin_pos[1]
    yyy = origin_node[0] + yy * v_row_y + xx * v_col_y
    xxx = origin_node[1] + yy * v_row_x + xx * v_col_x
    return yyy, xxx

Let’s create the grid we’ll be working with throughout this guide.

nrow = 50
ncol = 40
origin_pos = np.array((0.3,0.2))
origin_node = np.array((0., 0.))
v_row_y = 5.2
v_row_x = 1.2
v_col_y = -2.7
v_col_x = 7.1
grid_row_f64, grid_col_f64 = create_grid(nrow, ncol, origin_pos, origin_node, v_row_y, v_row_x, v_col_y, v_col_x, dtype=np.float64)

write_array(np.array([grid_row_f64, grid_col_f64]), dtype=np.float64, fileout=GRID_IN_F64)

grid_on_image.png

I/O Datasets

The GridR’s “chain”-called methods work on rasterio DatasetReader and DatasetWriter objets. We need to create a context with all the required datasets.

As far as the inputs are concerned, its quite easy : we just have to call open on read mode.

Concerning the output rasters, its quite more tedious : we have to define opening arguments such as the height, width, count, driver, dtype, … Let’s start simple :

  • We are going to apply the grid transform with a (1, 1) resolution with no other arguments (no window, no mask related options, …).

  • We will use the Geotiff format with a float64 data type.

  • We will work on the first band only

Given a (1, 1) resolution, the output raster’s shape corresponds to the grid’s shape.

The output raster opening arguments can be defined as following :

output_shape = grid_row_f64.shape

raster_out_open_args = {
    'driver': "GTiff",
    'dtype': np.float64,
    'height': output_shape[0],
    'width': output_shape[1],
    'count': 1
}

Now, let’s call the basic_grid_resampling_chain function within its context.

from gridr.chain.grid_resampling_chain import basic_grid_resampling_chain

output_raster_path = "./grid_resampling_chain_001_output_raster.tif"

with rasterio.open(GRID_IN_F64, 'r') as grid_in_ds, \
        rasterio.open(RASTER_IN, 'r') as array_src_ds, \
        rasterio.open(output_raster_path, "w", **raster_out_open_args) as array_out_ds:

    basic_grid_resampling_chain(
            grid_ds = grid_in_ds,
            grid_row_coords_band = 1,
            grid_col_coords_band = 2,
            grid_resolution = (1, 1),
            array_src_ds = array_src_ds,
            array_src_bands = 1,
            array_out_ds = array_out_ds,
            interp = "cubic",
            nodata_out = 0,
        )

Let’s open the output and show the image.

with rasterio.open(output_raster_path, "r") as ds:
    #print(ds.profile) 
    plot_im(ds.read(1), prefix='resampling_1_1')

resampling_1_1.png

That’s a tiny image result due to the subsampling the grid is performing !

Nevertheless, we can notice that the nodata_out value has been set where input image was not available.

Let’s apply a zoom by increasing the resolution to be (10, 10).

Now it’s interesting : we have to give to rasterio the output shape for that resolution. It is not quite complicated to compute from scratch but we will use here the grid_full_resolution_shape utility method.

from gridr.core.grid.grid_commons import grid_full_resolution_shape

grid_resolution = (8, 8)
output_shape = grid_full_resolution_shape(shape=grid_row_f64.shape, resolution=grid_resolution)

print(f"output shape : {output_shape}")
output shape : (393, 313)
raster_out_open_args = {
    'driver': "GTiff",
    'dtype': np.float64,
    'height': output_shape[0],
    'width': output_shape[1],
    'count': 1
}

# We overwrite on the previous output raster - no need to keep it
with rasterio.open(GRID_IN_F64, 'r') as grid_in_ds, \
        rasterio.open(RASTER_IN, 'r') as array_src_ds, \
        rasterio.open(output_raster_path, "w", **raster_out_open_args) as array_out_ds:

    basic_grid_resampling_chain(
            grid_ds = grid_in_ds,
            grid_row_coords_band = 1,
            grid_col_coords_band = 2,
            grid_resolution = grid_resolution,
            array_src_ds = array_src_ds,
            array_src_bands = 1,
            array_out_ds = array_out_ds,
            interp = "cubic",
            nodata_out = 0,
        )

resampling_8_8.png

Some data within the output has been set to nodata_out. In addtition, we can generate a binary output validity mask to indicate valid data points. This mask will adhere to GridR’s validity convention (where 1 represents valid data and 0 represents invalid data). To achieve this, we need to pass a dedicated output dataset for the mask.

mask_out_open_args = {
    'driver': "GTiff",
    'dtype': np.uint8, # <= save as unsigned int
    'height': output_shape[0],
    'width': output_shape[1],
    'count': 1, # <= the mask is common for all bands
    'nbits': 1, # <= GDAL option to save as true binary for less disk usage
}

output_mask_path = "./grid_resampling_chain_001_output_mask.tif"

# We overwrite on the previous output raster - no need to keep it
with rasterio.open(GRID_IN_F64, 'r') as grid_in_ds, \
        rasterio.open(RASTER_IN, 'r') as array_src_ds, \
        rasterio.open(output_raster_path, "w", **raster_out_open_args) as array_out_ds, \
        rasterio.open(output_mask_path, "w", **mask_out_open_args) as mask_out_ds:

    basic_grid_resampling_chain(
            grid_ds = grid_in_ds,
            grid_row_coords_band = 1,
            grid_col_coords_band = 2,
            grid_resolution = grid_resolution,
            array_src_ds = array_src_ds,
            array_src_bands = 1,
            array_out_ds = array_out_ds,
            interp = "cubic",
            nodata_out = 0,
            mask_out_ds = mask_out_ds, # <= set the dedicated dataset for mask out
        )

mask_out.png

Using a raster mask for the grid

Here we use an 8 bit raster having the same shape of the grids in order to devalidate points in the grid.

The basic_grid_resampling_chain function offers enhanced flexibility compared to the core array_grid_resampling method:

  • The data type used must be an 8-bit integer, which can be either unsigned or signed.

  • The user can specify the value representing valid points, with the constraint that this value must be positive. Any values different from this specified positive value will be considered invalid and subsequently masked.

grid_mask_in_path = './grid_resampling_chain_001_grid_mask_in_u8_1_0.tif'

grid_mask_in_band = 1
grid_mask_in_valid_value = 1
grid_mask_in_invalid_value = 0
grid_mask_dtype = np.uint8

# Set the invalid points to lie in a window
roi = np.array(((10,40), (5,100)))
grid_mask = np.full(grid_row_f64.shape, grid_mask_in_valid_value, dtype=np.uint8)
grid_mask[np.logical_and(
        np.logical_and(grid_row_f64 >= roi[0][0], grid_row_f64 <= roi[0][1]),
        np.logical_and(grid_col_f64 >= roi[1][0], grid_col_f64 <= roi[1][1]))] = grid_mask_in_invalid_value

write_array(grid_mask, dtype=grid_mask_dtype, fileout=grid_mask_in_path)

grid_mask_input.png

The figure above illustrates the grid mask:

  • Red nodes indicate invalid points as defined by the grid mask.

  • Orange nodes represent points that fall outside the image domain, which will also be considered invalid for resampling.

The grid mask is provided to the method using three arguments:

  • grid_mask_in_ds: The opened dataset corresponding to the mask.

  • grid_mask_in_unmasked_value: The specific value within the mask that should be considered valid (unmasked).

  • grid_mask_in_band: The band channel within the dataset that contains the mask information. This is particularly useful when the grid values and the mask are combined within a single dataset.

# We overwrite on the previous output raster - no need to keep it
with rasterio.open(GRID_IN_F64, 'r') as grid_in_ds, \
        rasterio.open(RASTER_IN, 'r') as array_src_ds, \
        rasterio.open(grid_mask_in_path, "r") as grid_mask_in_ds, \
        rasterio.open(output_raster_path, "w", **raster_out_open_args) as array_out_ds, \
        rasterio.open(output_mask_path, "w", **mask_out_open_args) as mask_out_ds:

    basic_grid_resampling_chain(
            grid_ds = grid_in_ds,
            grid_row_coords_band = 1,
            grid_col_coords_band = 2,
            grid_resolution = grid_resolution,
            array_src_ds = array_src_ds,
            array_src_bands = 1,
            array_out_ds = array_out_ds,
            interp = "cubic",
            nodata_out = 400, # <= just to illustrate another value
            mask_out_ds = mask_out_ds,
        
            grid_mask_in_ds = grid_mask_in_ds,
            grid_mask_in_unmasked_value = grid_mask_in_valid_value, # <= give the validity value
            grid_mask_in_band = 1, # <= give the band index
        )

grid_mask_output.png

Using different mask type and convention

Now, we’ll demonstrate the flexibility regarding mask data types and values.

grid_mask_in_path_i8 = './grid_resampling_chain_001_grid_mask_in_i8_0_m10.tif'

grid_mask_in_band_i8 = 1
grid_mask_in_valid_value_i8 = 0
grid_mask_in_invalid_value_i8 = -10
grid_mask_dtype_i8 = np.int8

roi = np.array(((10,40), (5,100)))
grid_mask_i8 = np.full(grid_row_f64.shape, grid_mask_in_valid_value_i8, dtype=grid_mask_dtype_i8)
grid_mask_i8[np.logical_and(
        np.logical_and(grid_row_f64 >= roi[0][0], grid_row_f64 <= roi[0][1]),
        np.logical_and(grid_col_f64 >= roi[1][0], grid_col_f64 <= roi[1][1]))] = grid_mask_in_invalid_value_i8

write_array(grid_mask_i8, dtype=grid_mask_dtype_i8, fileout=grid_mask_in_path_i8)

with rasterio.open(GRID_IN_F64, 'r') as grid_in_ds, \
        rasterio.open(RASTER_IN, 'r') as array_src_ds, \
        rasterio.open(grid_mask_in_path_i8, "r") as grid_mask_in_ds, \
        rasterio.open(output_raster_path, "w", **raster_out_open_args) as array_out_ds, \
        rasterio.open(output_mask_path, "w", **mask_out_open_args) as mask_out_ds:

    basic_grid_resampling_chain(
            grid_ds = grid_in_ds,
            grid_row_coords_band = 1,
            grid_col_coords_band = 2,
            grid_resolution = grid_resolution,
            array_src_ds = array_src_ds,
            array_src_bands = 1,
            array_out_ds = array_out_ds,
            interp = "cubic",
            nodata_out = 0,
            mask_out_ds = mask_out_ds,
        
            grid_mask_in_ds = grid_mask_in_ds,
            grid_mask_in_unmasked_value = grid_mask_in_valid_value_i8, # <= give the validity value
            grid_mask_in_band = 1, # <= give the band index
        )

grid_mask_out_int8.png

Using a raster mask for the image

Here we use an 8 bit raster having the same shape of the input array in order to devalidate points in the input array.

Basically if input data is devalidate, interpolation that involves those data must be considered as not valid resulting in output valued to nodata_out and the optional array_out_mask set to 0 instead of 1.

Similar to the grid masking, the basic_grid_resampling_chain method provides greater flexibility than the core function:

  • The data type must be an 8-bit integer, which can be either unsigned or signed.

  • Both valid and invalid values can be defined by the user

array_src_mask_validity_valid = 1
array_src_mask_validity_invalid = 0

# Create the mask raster - initiate to be valid
array_in_mask = np.full(mandrill[0].shape, array_src_mask_validity_valid, dtype=np.uint8)

# Define a region to be masked
masked_pos = slice(50,71), slice(150, 171) 
array_in_mask[masked_pos] = array_src_mask_validity_invalid

raster_mask_in_path_u8 = './grid_resampling_chain_001_raster_mask_in_u8_1_0.tif'

write_array(array_in_mask, dtype=np.uint8, fileout=raster_mask_in_path_u8)

array_in_mask_input.png

The input array mask is shown in the figure above as a small red square near the mandrill’s left eye.

Please note that throughout this guide, we will progressively add and demonstrate features to show how they can be used together. However, each feature can also be used independently.

The input array mask is provided to the method using three arguments:

  • array_src_mask_ds: The opened dataset corresponding to the mask.

  • array_src_mask_validity_pair: The tuple containing the values to consider as valid and invalid.

  • array_src_mask_band: The band channel within the dataset that contains the mask information. This is particularly useful when the image and the mask are combined within a single dataset.

with rasterio.open(GRID_IN_F64, 'r') as grid_in_ds, \
        rasterio.open(RASTER_IN, 'r') as array_src_ds, \
        rasterio.open(grid_mask_in_path, "r") as grid_mask_in_ds, \
        rasterio.open(raster_mask_in_path_u8, "r") as raster_mask_in_ds, \
        rasterio.open(output_raster_path, "w", **raster_out_open_args) as array_out_ds, \
        rasterio.open(output_mask_path, "w", **mask_out_open_args) as mask_out_ds:

    basic_grid_resampling_chain(
            grid_ds = grid_in_ds,
            grid_row_coords_band = 1,
            grid_col_coords_band = 2,
            grid_resolution = grid_resolution,
            array_src_ds = array_src_ds,
            array_src_bands = 1,
            array_out_ds = array_out_ds,
            interp = "cubic",
            nodata_out = 0, # <= just to illustrate
            mask_out_ds = mask_out_ds,
            grid_mask_in_ds = grid_mask_in_ds,
            grid_mask_in_unmasked_value = grid_mask_in_valid_value,
            grid_mask_in_band = 1,
        
            array_src_mask_ds = raster_mask_in_ds,
            array_src_mask_band = 1,
            array_src_mask_validity_pair = (array_src_mask_validity_valid,
                                            array_src_mask_validity_invalid),
        )

array_in_mask_output.png

Understanding I/O and Memory Management

In the previous section, we worked with rasterio Datasets. As you may have noticed, we did not perform any explicit read() or write() calls on these Datasets; these operations are handled internally by the basic_grid_resampling_chain method.

To give you more control over these I/O operations and memory usage, the method provides several parameters that can be fine-tuned:

  • io_strip_size

  • io_strip_size_target

  • tile_shape

  • ncpu (Note: Multiprocessing is not yet implemented for this parameter.)

Adjusting these parameters allows you to manage the memory footprint of your operations.

The basic_grid_resampling_chain Method’s Main Loop

The basic_grid_resampling_chain method processes the output raster in sequential line strips. For each strip, it computes all columns (or those defined by an optional window) for a contiguous set of rows. The maximum number of rows in a strip is controlled by the io_strip_size parameter.

You can define io_strip_size in two distinct modes:

  • GridRIOMode.INPUT: The target output strip size is calculated by multiplying the io_strip_size` value by the input grid’s row resolution.

  • GridRIOMode.OUTPUT: The io_strip_size parameter directly defines the target output strip size.

The main loop processes these output strips independently and sequentially:

  • Only the required input grid region to compute the current output strip is loaded into memory.

  • During each iteration, two single shared buffers are used. These buffers are allocated before the loop to match the larger shape of their respective strips: one for the input grid data and one for the output raster data.

It’s important to note that setting the io_strip_size parameter to 0 will trigger a unique strip computation, meaning the entire output raster will be processed at once.

Internal Strip Tiling

The tile_shape parameter adds another level of control, allowing a single strip to be processed in smaller tiles. This can significantly optimize memory usage for very large rasters, primarily by enhancing the efficiency of memory handling for the input image and its associated mask.

For each tile, the method calls the array_compute_resampling_grid_geometries function, located in the gridr.core.grid.grid_utils subpackage. This function returns the minimal bounding box required to read data from the input array_src_ds and ensures all necessary pixels are available for the resampling process.

Using tiles offers several key advantages:

  • It minimizes read operations to only the data required to process each tile.

  • If the data for a tile is completely unavailable or fully masked by the grid mask, the tile is skipped and filled with nodata_out.

  • Since read operations are rectangular, using smaller tiles limits the amount of unused data that is read, which is particularly beneficial when the grid is diagonal to the input image.

A key point to remember is that tiles do not share memory for the input data. This means the same input area can be read multiple times. Therefore, avoid using a tile_shape that is too small to prevent redundant reads and an excessive number of iterations.

It is important to note that setting the tile_shape parameter to None will disable tiling, causing the full strip to be processed in a single operation.

The “Basic” in basic_grid_resampling_chain

The “basic” prefix highlights a current limitation: a lack of automation. Setting key parameters can be a tedious process because the true memory cost is determined by the grid values, not by the parameters themselves. While future updates will likely automate this setup, it’s currently essential to set these parameters based on your knowledge of the grid values.

Debugging and Logging

To access verbose debugging information about the stripping and tiling process, you can provide a debug-enabled logger.Logger() object through the optional logger parameter.

import logging
from gridr.io.common import GridRIOMode

log_file = './gridr_resampling_chain_001_gridr_log.log'
fhandler = logging.FileHandler(filename=log_file, mode='w')
formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
fhandler.setFormatter(formatter)
gridr_logger = logging.getLogger('gridr')
gridr_logger.setLevel(logging.DEBUG) 
gridr_logger.addHandler(fhandler)
with rasterio.open(GRID_IN_F64, 'r') as grid_in_ds, \
        rasterio.open(RASTER_IN, 'r') as array_src_ds, \
        rasterio.open(grid_mask_in_path, "r") as grid_mask_in_ds, \
        rasterio.open(raster_mask_in_path_u8, "r") as raster_mask_in_ds, \
        rasterio.open(output_raster_path, "w", **raster_out_open_args) as array_out_ds, \
        rasterio.open(output_mask_path, "w", **mask_out_open_args) as mask_out_ds:

    basic_grid_resampling_chain(
            grid_ds = grid_in_ds,
            grid_row_coords_band = 1,
            grid_col_coords_band = 2,
            grid_resolution = grid_resolution,
            array_src_ds = array_src_ds,
            array_src_bands = 1,
            array_out_ds = array_out_ds,
            interp = "cubic",
            nodata_out = 0, # <= just to illustrate
            mask_out_ds = mask_out_ds,
            grid_mask_in_ds = grid_mask_in_ds,
            grid_mask_in_unmasked_value = grid_mask_in_valid_value,
            grid_mask_in_band = 1,
            array_src_mask_ds = raster_mask_in_ds,
            array_src_mask_band = 1,
            array_src_mask_validity_pair = (array_src_mask_validity_valid,
                                            array_src_mask_validity_invalid),
            io_strip_size = 200, # Here we know the total number or rows is 393 => we will have 2 strips
            io_strip_size_target = GridRIOMode.OUTPUT,
            tile_shape = (200, 160), # A strip will be maxed at 200 LINES x 313 COLUMNS => A (100, 160) tile will create 4 tiles per strip

            logger = gridr_logger,
        )

2026-03-13 13:52:02,013 - gridr - DEBUG - Grid shape : 50 rows x 40 columns
2026-03-13 13:52:02,013 - gridr - DEBUG - Grid mask shape : 50 rows x 40 columns
2026-03-13 13:52:02,013 - gridr - DEBUG - Grid mask unmasked value : 1
2026-03-13 13:52:02,013 - gridr - DEBUG - Grid mask dtype : uint8
2026-03-13 13:52:02,013 - gridr - DEBUG - Computed strip size (number of rows) : 200
2026-03-13 13:52:02,013 - gridr - DEBUG - Computed full output shape : 393 rows x 313 columns
2026-03-13 13:52:02,013 - gridr - DEBUG - No window given : the full grid is considered as window
2026-03-13 13:52:02,013 - gridr - DEBUG - Window : [[  0 392]
 [  0 312]]
2026-03-13 13:52:02,013 - gridr - DEBUG - Shape out : [393 313]
2026-03-13 13:52:02,014 - gridr - DEBUG - Read grid buffer shape : [26 40]
2026-03-13 13:52:02,014 - gridr - DEBUG - Read grid buffer shape 3 : [ 2 26 40]
2026-03-13 13:52:02,014 - gridr - DEBUG - Create float64 computation array buffer with shape [  1 200 313]
2026-03-13 13:52:02,014 - gridr - DEBUG - Create float64 computation array buffer with shape [  1 200 313] DONE
2026-03-13 13:52:02,015 - gridr - DEBUG - Create write mask buffer with shape [200 313]
2026-03-13 13:52:02,015 - gridr - DEBUG - Create write mask buffer with shape [200 313] DONE
2026-03-13 13:52:02,015 - gridr - DEBUG - Initializing interpolator optimized_bicubic
2026-03-13 13:52:02,015 - gridr - DEBUG - Required margins for interpolator optimized_bicubic : [[2 2]
 [2 2]]
2026-03-13 13:52:02,015 - gridr - DEBUG - Chunk 0 - chunk_win: [[  0 199]
 [  0 312]]
2026-03-13 13:52:02,016 - gridr - DEBUG - Chunk 0 - shape : (np.int64(200), np.int64(313))
2026-03-13 13:52:02,016 - gridr - DEBUG - Chunk 0 - grid read shape : (np.int64(26), np.int64(40))
2026-03-13 13:52:02,016 - gridr - DEBUG - Chunk 0 - buffer slices : (slice(0, 200, None), slice(0, 313, None))
2026-03-13 13:52:02,016 - gridr - DEBUG - Chunk 0 - buffer slices as win : [[  0 199]
 [  0 312]]
2026-03-13 13:52:02,016 - gridr - DEBUG - Chunk 0 - target write window : [[  0 199]
 [  0 312]]
2026-03-13 13:52:02,016 - gridr - DEBUG - Chunk 0 - resampling starts...
2026-03-13 13:52:02,016 - gridr - DEBUG - Chunk 0 - Tiled processing with tiles of 200 x 160
2026-03-13 13:52:02,016 - gridr - DEBUG - Chunk 0 - Number of tiles to process : 2
2026-03-13 13:52:02,016 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) - preparing args...
2026-03-13 13:52:02,016 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) - tile's chunk origin : [0 0]
2026-03-13 13:52:02,016 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) - tile window : [[  0 199]
 [  0 159]]
2026-03-13 13:52:02,017 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) - tile full-resolution window within the chunk's grid : [[  0 199]
 [  0 159]]
2026-03-13 13:52:02,017 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - Required margins for interpolator optimized_bicubic : [[2 2]
 [2 2]]
2026-03-13 13:52:02,017 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - Computing grid metrics... 
2026-03-13 13:52:02,017 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - SAFECHECK_SOURCE_BOUNDARIES : Computing source boundaries... 
2026-03-13 13:52:02,017 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - SAFECHECK_SOURCE_BOUNDARIES : GeometryBoundsF64(xmin=-2.37, xmax=169.62999999999997, ymin=-54.230000000000004, ymax=129.77)
2026-03-13 13:52:02,017 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - dst low res bounds : GeometryBoundsUsize(xmin=0, xmax=20, ymin=0, ymax=25) 
2026-03-13 13:52:02,017 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - src bounds : GeometryBoundsF64(xmin=-2.37, xmax=169.62999999999997, ymin=-54.230000000000004, ymax=129.77) 
2026-03-13 13:52:02,017 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - dst win : [[ 0 25]
 [ 0 20]]
2026-03-13 13:52:02,017 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - src read win (preliminary) : [[-55 130]
 [ -3 170]]
2026-03-13 13:52:02,017 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - src read win (preliminary) overflow : [[55  0]
 [ 3  0]]
2026-03-13 13:52:02,018 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - src read win before margin : [[  0 130]
 [  0 170]]
2026-03-13 13:52:02,018 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - src read win after margin : [[ -2 132]
 [ -2 172]]
2026-03-13 13:52:02,018 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - src read win required pad : [[2 0]
 [2 0]]
2026-03-13 13:52:02,018 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - src read win read : [[  0 132]
 [  0 172]] with shape (np.int64(133), np.int64(173))
2026-03-13 13:52:02,018 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - Computing required memory for array source read
2026-03-13 13:52:02,018 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - memory needed for array_src_win + margin band 1: 181632 bytes
2026-03-13 13:52:02,018 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - Computing required memory for array source mask read
2026-03-13 13:52:02,018 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - Memory required for array_src_mask_win + margin : 22704 bytes
2026-03-13 13:52:02,018 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - Memory required for array_src_mask + margin : 22704 bytes
2026-03-13 13:52:02,018 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - Total memory required for array_src_win and mask : 204336 bytes
2026-03-13 13:52:02,018 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - tile read buffer shape : [  1 135 175]
2026-03-13 13:52:02,018 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - Reading tiles...
2026-03-13 13:52:02,018 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - Reading source window for source band 1 - source window : [[  0 132]
 [  0 172]] - target indices : (0, slice(np.int64(2), np.int64(135), None), slice(np.int64(2), np.int64(175), None)) ...
2026-03-13 13:52:02,019 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - Reading source window for source band 1 - source window : [[  0 132]
 [  0 172]] - target indices : (0, slice(np.int64(2), np.int64(135), None), slice(np.int64(2), np.int64(175), None)) [DONE]
2026-03-13 13:52:02,019 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - Reading source window for source mask - source window : [[  0 132]
 [  0 172]] - target indices : (slice(np.int64(2), np.int64(135), None), slice(np.int64(2), np.int64(175), None)) ...
2026-03-13 13:52:02,019 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([  0, 160])) -  - Reading source window for source mask - source window : [[  0 132]
 [  0 172]] - target indices : (slice(np.int64(2), np.int64(135), None), slice(np.int64(2), np.int64(175), None)) [DONE]
2026-03-13 13:52:02,022 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) - preparing args...
2026-03-13 13:52:02,022 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) - tile's chunk origin : [0 0]
2026-03-13 13:52:02,022 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) - tile window : [[  0 199]
 [160 312]]
2026-03-13 13:52:02,022 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) - tile full-resolution window within the chunk's grid : [[  0 199]
 [160 312]]
2026-03-13 13:52:02,022 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - Required margins for interpolator optimized_bicubic : [[2 2]
 [2 2]]
2026-03-13 13:52:02,022 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - Computing grid metrics... 
2026-03-13 13:52:02,022 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - SAFECHECK_SOURCE_BOUNDARIES : Computing source boundaries... 
2026-03-13 13:52:02,022 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - SAFECHECK_SOURCE_BOUNDARIES : GeometryBoundsF64(xmin=139.62999999999997, xmax=304.53, ymin=-105.53000000000002, ymax=75.77000000000001)
2026-03-13 13:52:02,022 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - dst low res bounds : GeometryBoundsUsize(xmin=20, xmax=39, ymin=0, ymax=25) 
2026-03-13 13:52:02,023 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - src bounds : GeometryBoundsF64(xmin=139.62999999999997, xmax=304.53, ymin=-105.53000000000002, ymax=75.77000000000001) 
2026-03-13 13:52:02,023 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - dst win : [[ 0 25]
 [20 39]]
2026-03-13 13:52:02,023 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - src read win (preliminary) : [[-106   76]
 [ 139  305]]
2026-03-13 13:52:02,023 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - src read win (preliminary) overflow : [[106   0]
 [  0   0]]
2026-03-13 13:52:02,023 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - src read win before margin : [[  0  76]
 [139 305]]
2026-03-13 13:52:02,023 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - src read win after margin : [[ -2  78]
 [137 307]]
2026-03-13 13:52:02,023 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - src read win required pad : [[2 0]
 [0 0]]
2026-03-13 13:52:02,023 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - src read win read : [[  0  78]
 [137 307]] with shape (np.int64(79), np.int64(171))
2026-03-13 13:52:02,023 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - Computing required memory for array source read
2026-03-13 13:52:02,023 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - memory needed for array_src_win + margin band 1: 106080 bytes
2026-03-13 13:52:02,023 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - Computing required memory for array source mask read
2026-03-13 13:52:02,023 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - Memory required for array_src_mask_win + margin : 13260 bytes
2026-03-13 13:52:02,023 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - Memory required for array_src_mask + margin : 13260 bytes
2026-03-13 13:52:02,023 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - Total memory required for array_src_win and mask : 119340 bytes
2026-03-13 13:52:02,024 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - tile read buffer shape : [  1  81 171]
2026-03-13 13:52:02,024 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - Reading tiles...
2026-03-13 13:52:02,024 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - Reading source window for source band 1 - source window : [[  0  78]
 [137 307]] - target indices : (0, slice(np.int64(2), np.int64(81), None), slice(np.int64(0), np.int64(171), None)) ...
2026-03-13 13:52:02,024 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - Reading source window for source band 1 - source window : [[  0  78]
 [137 307]] - target indices : (0, slice(np.int64(2), np.int64(81), None), slice(np.int64(0), np.int64(171), None)) [DONE]
2026-03-13 13:52:02,024 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - Reading source window for source mask - source window : [[  0  78]
 [137 307]] - target indices : (slice(np.int64(2), np.int64(81), None), slice(np.int64(0), np.int64(171), None)) ...
2026-03-13 13:52:02,024 - gridr - DEBUG - Chunk 0 - tile ((0, np.int64(200)), array([160, 313])) -  - Reading source window for source mask - source window : [[  0  78]
 [137 307]] - target indices : (slice(np.int64(2), np.int64(81), None), slice(np.int64(0), np.int64(171), None)) [DONE]
2026-03-13 13:52:02,026 - gridr - DEBUG - Chunk 0 - write for full strip...
2026-03-13 13:52:02,026 - gridr - DEBUG - Chunk 0 - write ends.
2026-03-13 13:52:02,026 - gridr - DEBUG - Chunk 0 - write mask for full strip...
2026-03-13 13:52:02,026 - gridr - DEBUG - Chunk 0 - write mask ends.
2026-03-13 13:52:02,027 - gridr - DEBUG - Chunk 1 - chunk_win: [[200 392]
 [  0 312]]
2026-03-13 13:52:02,027 - gridr - DEBUG - Chunk 1 - shape : (np.int64(193), np.int64(313))
2026-03-13 13:52:02,027 - gridr - DEBUG - Chunk 1 - grid read shape : (np.int64(25), np.int64(40))
2026-03-13 13:52:02,027 - gridr - DEBUG - Chunk 1 - buffer slices : (slice(0, 193, None), slice(0, 313, None))
2026-03-13 13:52:02,027 - gridr - DEBUG - Chunk 1 - buffer slices as win : [[  0 192]
 [  0 312]]
2026-03-13 13:52:02,027 - gridr - DEBUG - Chunk 1 - target write window : [[200 392]
 [  0 312]]
2026-03-13 13:52:02,027 - gridr - DEBUG - Chunk 1 - resampling starts...
2026-03-13 13:52:02,027 - gridr - DEBUG - Chunk 1 - Tiled processing with tiles of 200 x 160
2026-03-13 13:52:02,027 - gridr - DEBUG - Chunk 1 - Number of tiles to process : 2
2026-03-13 13:52:02,027 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) - preparing args...
2026-03-13 13:52:02,028 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) - tile's chunk origin : [200   0]
2026-03-13 13:52:02,028 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) - tile window : [[  0 192]
 [  0 159]]
2026-03-13 13:52:02,028 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) - tile full-resolution window within the chunk's grid : [[  0 192]
 [  0 159]]
2026-03-13 13:52:02,028 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - Required margins for interpolator optimized_bicubic : [[2 2]
 [2 2]]
2026-03-13 13:52:02,028 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - Computing grid metrics... 
2026-03-13 13:52:02,028 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - SAFECHECK_SOURCE_BOUNDARIES : Computing source boundaries... 
2026-03-13 13:52:02,028 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - SAFECHECK_SOURCE_BOUNDARIES : GeometryBoundsF64(xmin=27.63, xmax=198.42999999999998, ymin=75.77000000000001, ymax=254.57)
2026-03-13 13:52:02,028 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - dst low res bounds : GeometryBoundsUsize(xmin=0, xmax=20, ymin=0, ymax=24) 
2026-03-13 13:52:02,028 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - src bounds : GeometryBoundsF64(xmin=27.63, xmax=198.42999999999998, ymin=75.77000000000001, ymax=254.57) 
2026-03-13 13:52:02,028 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - dst win : [[ 0 24]
 [ 0 20]]
2026-03-13 13:52:02,028 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - src read win (preliminary) : [[ 75 255]
 [ 27 199]]
2026-03-13 13:52:02,029 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - src read win (preliminary) overflow : [[0 0]
 [0 0]]
2026-03-13 13:52:02,029 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - src read win before margin : [[ 75 255]
 [ 27 199]]
2026-03-13 13:52:02,029 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - src read win after margin : [[ 73 257]
 [ 25 201]]
2026-03-13 13:52:02,029 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - src read win required pad : [[0 0]
 [0 0]]
2026-03-13 13:52:02,029 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - src read win read : [[ 73 257]
 [ 25 201]] with shape (np.int64(185), np.int64(177))
2026-03-13 13:52:02,029 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - Computing required memory for array source read
2026-03-13 13:52:02,029 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - memory needed for array_src_win + margin band 1: 259072 bytes
2026-03-13 13:52:02,029 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - Computing required memory for array source mask read
2026-03-13 13:52:02,029 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - Memory required for array_src_mask_win + margin : 32384 bytes
2026-03-13 13:52:02,029 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - Memory required for array_src_mask + margin : 32384 bytes
2026-03-13 13:52:02,030 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - Total memory required for array_src_win and mask : 291456 bytes
2026-03-13 13:52:02,030 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - tile read buffer shape : [  1 185 177]
2026-03-13 13:52:02,030 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - Reading tiles...
2026-03-13 13:52:02,030 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - Reading source window for source band 1 - source window : [[ 73 257]
 [ 25 201]] - target indices : (0, slice(np.int64(0), np.int64(185), None), slice(np.int64(0), np.int64(177), None)) ...
2026-03-13 13:52:02,030 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - Reading source window for source band 1 - source window : [[ 73 257]
 [ 25 201]] - target indices : (0, slice(np.int64(0), np.int64(185), None), slice(np.int64(0), np.int64(177), None)) [DONE]
2026-03-13 13:52:02,030 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - Reading source window for source mask - source window : [[ 73 257]
 [ 25 201]] - target indices : (slice(np.int64(0), np.int64(185), None), slice(np.int64(0), np.int64(177), None)) ...
2026-03-13 13:52:02,031 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([  0, 160])) -  - Reading source window for source mask - source window : [[ 73 257]
 [ 25 201]] - target indices : (slice(np.int64(0), np.int64(185), None), slice(np.int64(0), np.int64(177), None)) [DONE]
2026-03-13 13:52:02,033 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) - preparing args...
2026-03-13 13:52:02,033 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) - tile's chunk origin : [200   0]
2026-03-13 13:52:02,033 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) - tile window : [[  0 192]
 [160 312]]
2026-03-13 13:52:02,033 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) - tile full-resolution window within the chunk's grid : [[  0 192]
 [160 312]]
2026-03-13 13:52:02,034 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - Required margins for interpolator optimized_bicubic : [[2 2]
 [2 2]]
2026-03-13 13:52:02,034 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - Computing grid metrics... 
2026-03-13 13:52:02,034 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - SAFECHECK_SOURCE_BOUNDARIES : Computing source boundaries... 
2026-03-13 13:52:02,034 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - SAFECHECK_SOURCE_BOUNDARIES : GeometryBoundsF64(xmin=169.62999999999997, xmax=333.33, ymin=24.47, ymax=200.57)
2026-03-13 13:52:02,034 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - dst low res bounds : GeometryBoundsUsize(xmin=20, xmax=39, ymin=0, ymax=24) 
2026-03-13 13:52:02,034 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - src bounds : GeometryBoundsF64(xmin=169.62999999999997, xmax=333.33, ymin=24.47, ymax=200.57) 
2026-03-13 13:52:02,034 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - dst win : [[ 0 24]
 [20 39]]
2026-03-13 13:52:02,034 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - src read win (preliminary) : [[ 24 201]
 [169 334]]
2026-03-13 13:52:02,034 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - src read win (preliminary) overflow : [[0 0]
 [0 0]]
2026-03-13 13:52:02,034 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - src read win before margin : [[ 24 201]
 [169 334]]
2026-03-13 13:52:02,034 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - src read win after margin : [[ 22 203]
 [167 336]]
2026-03-13 13:52:02,034 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - src read win required pad : [[0 0]
 [0 0]]
2026-03-13 13:52:02,034 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - src read win read : [[ 22 203]
 [167 336]] with shape (np.int64(182), np.int64(170))
2026-03-13 13:52:02,035 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - Computing required memory for array source read
2026-03-13 13:52:02,035 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - memory needed for array_src_win + margin band 1: 244712 bytes
2026-03-13 13:52:02,035 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - Computing required memory for array source mask read
2026-03-13 13:52:02,035 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - Memory required for array_src_mask_win + margin : 30589 bytes
2026-03-13 13:52:02,035 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - Memory required for array_src_mask + margin : 30589 bytes
2026-03-13 13:52:02,035 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - Total memory required for array_src_win and mask : 275301 bytes
2026-03-13 13:52:02,035 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - tile read buffer shape : [  1 182 170]
2026-03-13 13:52:02,035 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - Reading tiles...
2026-03-13 13:52:02,035 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - Reading source window for source band 1 - source window : [[ 22 203]
 [167 336]] - target indices : (0, slice(np.int64(0), np.int64(182), None), slice(np.int64(0), np.int64(170), None)) ...
2026-03-13 13:52:02,035 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - Reading source window for source band 1 - source window : [[ 22 203]
 [167 336]] - target indices : (0, slice(np.int64(0), np.int64(182), None), slice(np.int64(0), np.int64(170), None)) [DONE]
2026-03-13 13:52:02,035 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - Reading source window for source mask - source window : [[ 22 203]
 [167 336]] - target indices : (slice(np.int64(0), np.int64(182), None), slice(np.int64(0), np.int64(170), None)) ...
2026-03-13 13:52:02,036 - gridr - DEBUG - Chunk 1 - tile ((0, np.int64(193)), array([160, 313])) -  - Reading source window for source mask - source window : [[ 22 203]
 [167 336]] - target indices : (slice(np.int64(0), np.int64(182), None), slice(np.int64(0), np.int64(170), None)) [DONE]
2026-03-13 13:52:02,038 - gridr - DEBUG - Chunk 1 - write for full strip...
2026-03-13 13:52:02,038 - gridr - DEBUG - Chunk 1 - write ends.
2026-03-13 13:52:02,038 - gridr - DEBUG - Chunk 1 - write mask for full strip...
2026-03-13 13:52:02,038 - gridr - DEBUG - Chunk 1 - write mask ends.

Extended features

Separate datasets for grid row and columns

As you may have notices, the basic_grid_resampling_chain provides two distinct arguments to set the grids dataset, namely grid_ds and grid_col_ds.

If grid_col_ds is set to None or to the same value as grid_ds, the same dataset will be used for both the row and column grid values. But you can adress two distincts files for the row and column values.

Using geometry masks

The basic_grid_resampling_chain introduces a significant new feature compared to the core array_grid_resampling method: the ability to define ‘geometry’ masks.

Essentially, these geometry masks are rasterized at a tile level before the array_grid_resampling method is called. All masks related to the source image array are then merged into a single, unified mask, which is subsequently passed to the core resampling function.

Let’s illustrate how to implement this.

from shapely.geometry import Polygon

geometry_valid = Polygon([(40, 10), (300, 30), (270, 200), (40, 200)])
# call an internal function to create a star polygon
invalid_polygon = create_star_polygon(100, 100, 50)

geometry_mask_input.png

The figure above displays the new valid geometry (the green shape) and invalid geometry (the red star).

with rasterio.open(GRID_IN_F64, 'r') as grid_in_ds, \
        rasterio.open(RASTER_IN, 'r') as array_src_ds, \
        rasterio.open(grid_mask_in_path, "r") as grid_mask_in_ds, \
        rasterio.open(raster_mask_in_path_u8, "r") as raster_mask_in_ds, \
        rasterio.open(output_raster_path, "w", **raster_out_open_args) as array_out_ds, \
        rasterio.open(output_mask_path, "w", **mask_out_open_args) as mask_out_ds:

    basic_grid_resampling_chain(
            grid_ds = grid_in_ds,
            grid_row_coords_band = 1,
            grid_col_coords_band = 2,
            grid_resolution = grid_resolution,
            array_src_ds = array_src_ds,
            array_src_bands = 1,
            array_out_ds = array_out_ds,
            interp = "cubic",
            nodata_out = 0, # <= just to illustrate
            mask_out_ds = mask_out_ds,
            grid_mask_in_ds = grid_mask_in_ds,
            grid_mask_in_unmasked_value = grid_mask_in_valid_value,
            grid_mask_in_band = 1,
            array_src_mask_ds = raster_mask_in_ds,
            array_src_mask_band = 1,
            array_src_mask_validity_pair = (array_src_mask_validity_valid,
                                            array_src_mask_validity_invalid),
            array_src_geometry_origin=(0., 0.),
            array_src_geometry_pair = [geometry_valid, invalid_polygon],
        )

geometry_mask_output.png

Shifting the grid coordinates

The array_grid_resampling function includes the array_in_origin parameter, which enables shifting the coordinate system’s origin. From an implementation perspective, this shift is directly applied to the full resolution calculated target grid coordinates.

The array_in_origin parameter is primarily used in the basic_grid_resampling_chain to establish the origin point for the current tile subregion. However, it’s important to note that this parameter doesn’t influence the grid geometry computations. Consequently, it is not used to apply a global bias to the current grid strip within the basic_grid_resampling_chain method.

For global coordinate adjustments, the basic_grid_resampling_chain offers the grid_shift parameter. This parameter accepts an optional (shift_row, shift_col) integers or floating points tuple, allowing specific shifts to be applied to row and column coordinates before performing any grid calculations.

with rasterio.open(GRID_IN_F64, 'r') as grid_in_ds, \
        rasterio.open(RASTER_IN, 'r') as array_src_ds, \
        rasterio.open(grid_mask_in_path, "r") as grid_mask_in_ds, \
        rasterio.open(raster_mask_in_path_u8, "r") as raster_mask_in_ds, \
        rasterio.open(output_raster_path, "w", **raster_out_open_args) as array_out_ds, \
        rasterio.open(output_mask_path, "w", **mask_out_open_args) as mask_out_ds:

    basic_grid_resampling_chain(
            grid_ds = grid_in_ds,
            grid_row_coords_band = 1,
            grid_col_coords_band = 2,
            grid_resolution = grid_resolution,
            array_src_ds = array_src_ds,
            array_src_bands = 1,
            array_out_ds = array_out_ds,
            interp = "cubic",
            nodata_out = 0,
            win = None,
            grid_shift = (50.5, 30.),  # <= set the shift through the `grid_shift` parameter
            mask_out_ds = mask_out_ds,
            grid_mask_in_ds = grid_mask_in_ds,
            grid_mask_in_unmasked_value = grid_mask_in_valid_value,
            grid_mask_in_band = 1,
            array_src_mask_ds = raster_mask_in_ds,
            array_src_mask_band = 1,
            array_src_mask_validity_pair = (array_src_mask_validity_valid,
                                            array_src_mask_validity_invalid),
            array_src_geometry_origin=(0., 0.),
            array_src_geometry_pair = [geometry_valid, invalid_polygon],
        )

grid_shift_output.png

As observed, the bias adjustment is applied exclusively to the grid coordinates, leaving the grid’s mask unchanged.

Limit computation to a window

Like the array_grid_resampling method , the basic_grid_resampling_chain provides the win parameter in order to set a window to limit the region to compute.

The indices used in the win definition refers to indices considering the full resolution output.

win = np.asarray([[120, 280], [70, 240]])

apply_window_input.png

The target window is the region defined by the orange rectangle in the figure above.

When we set the window, we need to adjust the output shape of both the image and its corresponding mask to reflect this new, smaller area.

raster_out_open_args = {
    'driver': "GTiff",
    'dtype': np.float64,
    'height': win[0][1] - win[0][0] + 1,
    'width': win[1][1] - win[1][0] + 1,
    'count': 1
}

mask_out_open_args = {
    'driver': "GTiff",
    'dtype': np.uint8,
    'height': win[0][1] - win[0][0] + 1,
    'width': win[1][1] - win[1][0] + 1,
    'count': 1,
    'nbits': 1,
}
with rasterio.open(GRID_IN_F64, 'r') as grid_in_ds, \
        rasterio.open(RASTER_IN, 'r') as array_src_ds, \
        rasterio.open(grid_mask_in_path, "r") as grid_mask_in_ds, \
        rasterio.open(raster_mask_in_path_u8, "r") as raster_mask_in_ds, \
        rasterio.open(output_raster_path, "w", **raster_out_open_args) as array_out_ds, \
        rasterio.open(output_mask_path, "w", **mask_out_open_args) as mask_out_ds:

    basic_grid_resampling_chain(
            grid_ds = grid_in_ds,
            grid_row_coords_band = 1,
            grid_col_coords_band = 2,
            grid_resolution = grid_resolution,
            array_src_ds = array_src_ds,
            array_src_bands = 1,
            array_out_ds = array_out_ds,
            interp = "cubic",
            nodata_out = 0,
            win = win, # <= set the window through the `win` parameter
            mask_out_ds = mask_out_ds,
            grid_mask_in_ds = grid_mask_in_ds,
            grid_mask_in_unmasked_value = grid_mask_in_valid_value,
            grid_mask_in_band = 1,
            array_src_mask_ds = raster_mask_in_ds,
            array_src_mask_band = 1,
            array_src_mask_validity_pair = (array_src_mask_validity_valid,
                                            array_src_mask_validity_invalid),
            array_src_geometry_origin=(0., 0.),
            array_src_geometry_pair = [geometry_valid, invalid_polygon],
        )

apply_window_output.png

The output you see now shows only the area within the target window, with all masks correctly applied to that specific region.