Slight error after converting coordinates between CRS

Question

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)

Answer 1

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)