Lesson 8. Reproject Raster Data Python


Learning objectives

After completing this tutorial, you will be able to:

  • Reproject a raster using rasterio

What You Need

You will need a computer with internet access to complete this lesson.

Download Colorado Flood Teaching Data Subset data

Reprojecting

Sometimes you will work with multiple rasters and they might not always be in the same projections. You will need to reproject the raster so they are in the same coordinate reference system.

Reproject using Rasterio

import rasterio as rio
import os
import earthpy as et
import numpy as np
from rasterio.warp import calculate_default_transform, reproject, Resampling
os.chdir(os.path.join(et.io.HOME, 'earth-analytics'))
lidar_dem = rio.open('data/colorado-flood/spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM.tif')
print(lidar_dem.meta)
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.4028234663852886e+38, 'width': 4000, 'height': 2000, 'count': 1, 'crs': CRS({'init': 'epsg:32613'}), 'transform': Affine(1.0, 0.0, 472000.0,
       0.0, -1.0, 4436000.0)}
dst_crs = 'EPSG:3857' # CRS for web meractor 

with rio.open('data/colorado-flood/spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM.tif') 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 rio.open('data/colorado-flood/spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM_repoject.tif', 'w', **kwargs) as dst:
        for i in range(1, src.count + 1):
            reproject(
                source=rio.band(src, i),
                destination=rio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=Resampling.nearest)
lidar_dem2 = rio.open('data/colorado-flood/spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM_repoject.tif')
print(lidar_dem2.meta)
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.4028234663852886e+38, 'width': 4004, 'height': 2020, 'count': 1, 'crs': CRS({'init': 'epsg:3857'}), 'transform': Affine(1.3063652820086313, 0.0, -11725101.307458913,
       0.0, -1.3063652820086313, 4876690.453258085)}
lidar_dem2.close()

If you have many raster files to re-project the rasterio method has several lines of code that could get repetative to type. Therefore your instructions have wrapped the rasterio reproject code into a function called reproject_et notice that this function contains all of the same code, but allows the input path, the output path, and the new CRS components to change every time the function is called.

def reproject_et(inpath, outpath, new_crs):
    dst_crs = new_crs # CRS for web meractor 

    with rio.open(inpath) 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 rio.open(outpath, 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rio.band(src, i),
                    destination=rio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.nearest)
reproject_et(inpath = 'data/colorado-flood/spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM.tif', 
              outpath = 'data/colorado-flood/spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM_repoject2.tif', 
              new_crs = 'EPSG:3857')
# Check to make sure function worked, then close raster
lidar_dem3 = rio.open('data/colorado-flood/spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM_repoject2.tif')
print(lidar_dem3.meta)
lidar_dem3.close()
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.4028234663852886e+38, 'width': 4004, 'height': 2020, 'count': 1, 'crs': CRS({'init': 'epsg:3857'}), 'transform': Affine(1.3063652820086313, 0.0, -11725101.307458913,
       0.0, -1.3063652820086313, 4876690.453258085)}