How to reproject single and multiple rasters with Python and Rasterio - Tutorial

Raster reprojection is a common task on GIS analysis however to do it with only Python commands has some challenges. We have done an applied example of raster reprojection for single and multiple rasters from WGS 84 UTM to WGS 84 Geographic. The codes work on monoband and multiband rasters and can reproject from to any projection by specifiyin it EPSG code.

The tutorial is done on a conda enviroment that has the geospatial libraries installed. You can set up the environment by following this tutorial:

https://hatarilabs.com/ih-en/how-to-install-python-geospatial-libraries-gdal-fiona-rasterio-etc-under-a-conda-env-in-windows

Tutorial

Code

This is the code for the single raster reprojection:

import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

#open source raster
srcRst = rasterio.open('../singleRaster/landsatImage.tif')

print("source raster crs:")
print(srcRst.crs)

dstCrs = {'init': 'EPSG:4326'}
print("destination raster crs:")
print(dstCrs)

source raster crs:
EPSG:32614
destination raster crs:
{'init': 'EPSG:4326'}

#calculate transform array and shape of reprojected raster
transform, width, height = calculate_default_transform(
        srcRst.crs, dstCrs, srcRst.width, srcRst.height, *srcRst.bounds)

print("transform array of source raster")
print(srcRst.transform)

print("transform array of destination raster")
print(transform)

transform array of source raster
| 10.00, 0.00, 567810.00|
| 0.00,-10.00, 2127440.00|
| 0.00, 0.00, 1.00|
transform array of destination raster
| 0.00, 0.00,-98.36|
| 0.00,-0.00, 19.24|
| 0.00, 0.00, 1.00|

#working of the meta for the destination raster
kwargs = srcRst.meta.copy()
kwargs.update({
        'crs': dstCrs,
        'transform': transform,
        'width': width,
        'height': height
    })

#open destination raster
dstRst = rasterio.open('../outputRaster/landsatImageWgs84.tif', 'w', **kwargs)

#reproject and save raster band data
for i in range(1, srcRst.count + 1):
    reproject(
        source=rasterio.band(srcRst, i),
        destination=rasterio.band(dstRst, i),
        #src_transform=srcRst.transform,
        src_crs=srcRst.crs,
        #dst_transform=transform,
        dst_crs=dstCrs,
        resampling=Resampling.nearest)

#close destination raster
dstRst.close()

Input data

You can download the input data from this link.