3

I have two rasters of the same geographic area, but using different CRS. My aim is to match them so I can stack them and perform per-pixel analysis between the two.

To do that, I'm trying to convert the "matching" raster's bounds to the input's CRS to calculate the default transformation to apply to the reprojection.

The match raster CRS is EPSG:3857, and my input's exact CRS string as returned by Rasterio is

PROJCS["ETRS89_ETRS_LAEA",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101004,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",52],PARAMETER["longitude_of_center",10],PARAMETER["false_easting",4321000],PARAMETER["false_northing",3210000],UNIT["metre",1],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

I'm using the following code to perform the conversion :

transformer = pyproj.Transformer.from_crs(match.crs, src.crs, always_xy=True)
new_bounds = transformer.transform_bounds(*match.bounds)

The result, however, is slightly off, as you can see here :

The bottom-left point of the matching raster

The bottom-left point after conversion in the input CRS

The difference is not massive, but big enough to prevent me from properly stacking the two rasters after reprojection.

I'm not sure what causes this issue, I've tried both PyProj (3.4.0) and Rasterio, thinking it could be a version issue, but no luck.

Is there a way to fix this conversion or, alternatively, another way of stacking the two rasters together (knowing that they don't have the same resolution either)

EDIT : I've added the full code of the method I'm using for more context

def match(
        infile: Union[str, pathlib.Path], 
        matchfile: Union[str, pathlib.Path], 
        outfile: Union[str, pathlib.Path],
        resampling_method: Resampling=Resampling.nearest
    ):
    """Reprojects the input raster so that it matches the 'matching' one.

    Parameters
    ----------
    infile : Union[str, pathlib.Path]
        The input raster
    matchfile : Union[str, pathlib.Path]
        The raster to match
    outfile : Union[str, pathlib.Path]
        The path where to save the resulting reprojection
    resampling_method : Resampling, optional
        The resampling method to use, by default Resampling.nearest
    """
    with rasterio.open(infile) as src:
        with rasterio.open(matchfile) as match:
            transformer = pyproj.Transformer.from_crs(match.crs, src.crs, always_xy=True)
            new_bounds = transformer.transform_bounds(*match.bounds)

            dst_transform, dst_width, dst_height = calculate_default_transform(
                src.crs,
                match.crs,
                src.width,
                src.height,
                *new_bounds,
                dst_width=match.width,
                dst_height=match.height
            )

        dst_kwargs = src.meta.copy()
        dst_kwargs.update({"crs": match.crs,
                           "transform": dst_transform,
                           "width": dst_width,
                           "height": dst_height})

        with rasterio.open(outfile, "w", **dst_kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=dst_transform,
                    dst_crs=match.crs,
                    resampling=resampling_method)

1 Answer 1

4

You do not need to pre-calculate the transform to do the resample, you can just give it the transform of the raster you want to match and reproject will do what is needed. That should ensure the resulting pixels match exactly.

Like so:

def match(
        infile: Union[str, pathlib.Path], 
        matchfile: Union[str, pathlib.Path], 
        outfile: Union[str, pathlib.Path],
        resampling_method: Resampling=Resampling.nearest
    ):
    """Reprojects the input raster so that it matches the 'matching' one.

    Parameters
    ----------
    infile : Union[str, pathlib.Path]
        The input raster
    matchfile : Union[str, pathlib.Path]
        The raster to match
    outfile : Union[str, pathlib.Path]
        The path where to save the resulting reprojection
    resampling_method : Resampling, optional
        The resampling method to use, by default Resampling.nearest
    """
    with rasterio.open(infile) as src:
        with rasterio.open(matchfile) as match_src:
            dst_kwargs = src.meta.copy()
            match_profile = match_src.profile
            dst_kwargs.update({"crs"       : match_profile['crs'],
                               "transform" : match_profile['transform'],
                               "width"     : match_profile['width'],,
                               "height"    : match_profile['height'],})

        with rasterio.open(outfile, "w", **dst_kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=match_profile['transform'],
                    dst_crs      =match_profile['crs'],
                    resampling=resampling_method)
1
  • There's a slight typo in the reprojection I believe, it should be dst_crs=match_profile['crs'] but this perfectly solves my issue, thank you Oct 20, 2022 at 14:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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