Rasterio – Fixing Rasterio Reprojection Issues

Question

I would like to reproject a raster using the coordinate system of a shapefile as my destination coordinate system. I am using one of the examples in the rasterio docs as inspiration:.

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
import geopandas

src_file = 'output/out.tif'
dst_file = 'output/out_reproj.tif'
dst_crs = geopandas.read_file('confidence.shp').crs['init']

with rasterio.open(src_file) as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
    })

    with rasterio.open(dst_file, 'w', **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=transform,
                dst_crs=dst_crs,
                resampling=Resampling.nearest)

The shapefile is in EPSG::3035 and the input raster is in EPSG::4326. The output raster is still in EPSG::4326.

Am I doing anything wrong in the script?

Answer 1

I just figured it out. You can't go passing ingdf.crs formatted dicts or their contents into rasterio metadata. It won't work. You actually have to build up a rasterio crs object so it will recognize it:

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
import geopandas

src_file = 'output/out.tif'
dst_file = 'output/out_reproj.tif'
gdf = geopandas.read_file('confidence.shp')
dst_crs = rasterio.crs.CRS.from_dict(gdf.crs)

with rasterio.open(src_file) as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
    })

    with rasterio.open(dst_file, 'w', **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=transform,
                dst_crs=dst_crs,
                resampling=Resampling.nearest)

The output will then be correctly projected.