diff --git a/rasterio/merge.py b/rasterio/merge.py index 7b84797e..0f870d46 100644 --- a/rasterio/merge.py +++ b/rasterio/merge.py @@ -11,7 +11,7 @@ import numpy as np import rasterio from rasterio.coords import disjoint_bounds from rasterio.enums import Resampling -from rasterio.errors import RasterioDeprecationWarning +from rasterio.errors import RasterioDeprecationWarning, RasterioError from rasterio import windows from rasterio.transform import Affine @@ -219,6 +219,7 @@ def merge( with dataset_opener(datasets[0]) as first: first_profile = first.profile + first_crs = first.crs first_res = first.res nodataval = first.nodatavals[0] dt = first.dtypes[0] @@ -324,10 +325,17 @@ def merge( # This approach uses the maximum amount of memory to solve the # problem. Making it more efficient is a TODO. + # 0. Precondition checks + # - Check that source is within destination bounds + # - Check that CRS is same + if disjoint_bounds((dst_w, dst_s, dst_e, dst_n), src.bounds): logger.debug("Skipping source: src=%r, window=%r", src) continue + if first_crs != src.crs: + raise RasterioError(f"CRS mismatch with source: {dataset}") + # 1. Compute spatial intersection of destination and source src_w, src_s, src_e, src_n = src.bounds diff --git a/tests/test_merge.py b/tests/test_merge.py index e8364376..01c0d3f6 100644 --- a/tests/test_merge.py +++ b/tests/test_merge.py @@ -9,6 +9,8 @@ import pytest import affine import rasterio from rasterio.merge import merge +from rasterio.crs import CRS +from rasterio.errors import RasterioError # Non-coincident datasets test fixture. # Three overlapping GeoTIFFs, two to the NW and one to the SE. @@ -42,6 +44,19 @@ def test_data_dir_overlapping(tmp_path): return tmp_path +def test_different_crs(test_data_dir_overlapping): + inputs = [x.name for x in test_data_dir_overlapping.iterdir()] + + # Create new raster with different crs + with rasterio.open(test_data_dir_overlapping.joinpath(inputs[-1])) as ds_src: + kwds = ds_src.profile + kwds['crs'] = CRS.from_epsg(3499) + with rasterio.open(test_data_dir_overlapping.joinpath("new.tif"), 'w', **kwds) as ds_out: + ds_out.write(ds_src.read()) + with pytest.raises(RasterioError): + result = merge(list(test_data_dir_overlapping.iterdir())) + + @pytest.mark.parametrize( "method,value", [("first", 1), ("last", 2), ("min", 1), ("max", 3), ("sum", 6), ("count", 3)],