Skip to content

rasterio.merge incompatible with WarpedVRT in 1.4.0 #3196

@adamjstewart

Description

@adamjstewart

Background

TorchGeo is a PyTorch domain library for geospatial data, and wraps around rasterio for file I/O. A common use case is to take a collection of arbitrary raster files and use windowed reading to load small image patches from larger raster scenes to then train ML models.

  • In the case of neighboring files, we automatically merge these scenes (rasterio.merge.merge)
  • In the case of mixed-CRS, we automatically reproject them to a common CRS (rasterio.vrt.WarpedVRT)
  • In the case where we iterate over these scenes in a grid-like fashion and the scene size is not a multiple of the patch size, we sample outside the bounds (boundless=True)

Expected behavior and actual behavior.

In prior rasterio releases, rasterio.merge.merge supported rasterio.vrt.WarpedVRT objects. However, rasterio 1.4.0 seems to have (unintentionally?) broken support for this. Instead of reprojecting, merging, and windowed reading, we now see the following error:

Traceback (most recent call last):
  File ".../test.py", line 33, in <module>
    merge([x_vrt, y_vrt])
  File ".../rasterio/merge.py", line 431, in merge
    temp_src = src.read(
               ^^^^^^^^^
  File "rasterio/_warp.pyx", line 1326, in rasterio._warp.WarpedVRTReaderBase.read
ValueError: WarpedVRT does not permit boundless reads

This occurs even when boundless reads are not requested. Since we would also like boundless reads, I'm wondering if this was an intentional change, and what the alternative might be. Should we give up on merging? Or on reprojecting? Or on boundless reads? Or should we manually perform all of these operations instead of leaving it up to rasterio.merge.merge?

Steps to reproduce the problem.

The following code works with rasterio 1.3 but not 1.4:

import numpy as np
import rasterio as rio
from affine import Affine
from rasterio.crs import CRS
from rasterio.merge import merge
from rasterio.vrt import WarpedVRT

size = 2
crs = CRS.from_epsg(3857)

x = np.zeros((size, size))
y = np.zeros((size, size))

x_transform = Affine(1, 0, 0, 0, -1, 2)
y_transform = Affine(1, 0, 2, 0, -1, 2)

profile = {
    'width': size,
    'height': size,
    'count': 1,
    'dtype': np.uint8,
    'crs': CRS.from_epsg(4326),
}

with rio.open('x.tif', 'w', transform=x_transform, **profile) as x_dst:
    x_dst.write(x, 1)

with rio.open('y.tif', 'w', transform=y_transform, **profile) as y_dst:
    y_dst.write(y, 1)

with rio.open('x.tif') as x_src, rio.open('y.tif') as y_src:
    with WarpedVRT(x_src, crs=crs) as x_vrt, WarpedVRT(y_src, crs=crs) as y_vrt:
        merge([x_vrt, y_vrt])

Environment Information

rasterio info:
  rasterio: 1.4.0
      GDAL: 3.9.2
      PROJ: 9.4.1
      GEOS: 3.13.0
 PROJ DATA: /Users/Adam/spack/opt/spack/darwin-sequoia-m2/apple-clang-16.0.0/proj-9.4.1-jz4g7qt3jp6pbuwhkre2mtb4sfdrakyt/share/proj
 GDAL DATA: None

System:
    python: 3.11.9 (main, Sep 23 2024, 12:36:15) [Clang 16.0.0 (clang-1600.0.26.3)]
executable: /Users/Adam/spack/opt/spack/darwin-sequoia-m2/apple-clang-16.0.0/python-venv-1.0-mcw7wh5o6eq6ztq244pxsof6xdozpwk3/bin/python3
   machine: macOS-15.0-arm64-arm-64bit

Python deps:
    affine: 2.1.0
     attrs: 23.1.0
   certifi: 2023.07.22
     click: 8.1.7
     cligj: 0.7.2
    cython: None
     numpy: 2.1.1
click-plugins: None
setuptools: 69.2.0

Installation Method

Installed from source using Spack.

@calebrob6

Metadata

Metadata

Assignees

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions