06 Raster2 Demo

Reprojection, Clipping, Sampling, Zonal Stats

UW Geospatial Data Analysis
CEE467/CEWA567
David Shean, Eric Gagliano, Quinn Brencher

Objectives

• Demonstrate multiple approaches for “on-the-fly” raster data download

• Understand additional fundamental raster processing/analysis:

  • Reprojection

  • Clipping

  • Interpolation and sampling strategies

• Combine vector points and polygons with rasters for zonal statistics

• Understand processing strategies, derivative products, and common applications for a fundamental raster data product: DEMs

  • Slope and Aspect

  • Contour generation

  • Volume estimation (cut/fill analysis)

What is a DEM?

DEM = Digital Elevation Model

A generic term for a 2D raster grid with values representing surface elevation above some datum (e.g., WGS84 ellispoid or a geoid model representing mean sea level). Sometimes called 2.5D, as it’s not a true 3D dataset containing some value (e.g., temperature) at each (x,y,z) point.

There are subtypes:

• DSM = Digital Surface Model (“first-return” model includes top of canopy, buildings, etc.)

• DTM = Digital Terrain Model (bare ground model, with canopy, buildings, etc. removed)

Great resource on LiDAR and derivative products:

• https://www.earthdatascience.org/courses/use-data-open-source-python/data-stories/what-is-lidar-data/

• https://www.earthdatascience.org/courses/use-data-open-source-python/data-stories/what-is-lidar-data/lidar-chm-dem-dsm/

Airborne LiDAR example

• https://www.dnr.wa.gov/lidar

• https://wadnr.maps.arcgis.com/apps/Cascade/index.html?appid=b93c17aa1ef24669b656dbaea009b5ce

• https://wadnr.maps.arcgis.com/apps/Cascade/index.html?appid=36b4887370d141fcbb35392f996c82d9

• https://lidarportal.dnr.wa.gov/

SRTM

SRTM = Shuttle Radar Topography Mission

• https://en.wikipedia.org/wiki/Shuttle_Radar_Topography_Mission

• https://www2.jpl.nasa.gov/srtm/

This week, we’ll play with the landmark SRTM dataset. We briefly introduced this during Lab03, as I sampled the SRTM products for the GLAS point locations and included as the dem_z column in the csv.

• Collected February 11-22, 2000 (winter)

• Single-track InSAR (interferometric synthetic aperture radar) instrument

• Coverage: 56°S to 60°N (the shuttle orbit, plus radar look direction)

• Tiled in 1x1° raster data at different resolutions:

  • 1-arcsecond (~30 m)

  • 3-arcsecond (~90 m)

• Default elevation values in the SRTM tiles are relative to the EGM96 geoid (approximates mean sea level), not the WGS84 ellipsoid (as with the GLAS points)

For this lab, we’ll use the 3-arcsec (90 m) SRTM product for WA state to learn some new concepts. This is a relatively small dataset, with a limited number of 1x1° tiles required for WA state. However, the approaches we’ll learn (e.g., API subsetting, using vrt datasets), scale to larger datasets that are too big to fit in memory (like operations on the global SRTM dataset).

Copernicus DEM

• https://spacedata.copernicus.eu/web/cscda/dataset-details?articleId=394198

• https://spacedata.copernicus.eu/documents/20126/0/GEO1988-CopernicusDEM-SPE-002_ProductHandbook_I1.00.pdf

• https://spacedata.copernicus.eu/documents/20126/0/GEO1988-CopernicusDEM-RP-001_ValidationReport_V1.0.pdf

Interactive discussion topics

• Mixing command line utilities and Python API code

• Raster reprojection

  • https://support.esri.com/en/technical-article/000008915

  • https://en.wikipedia.org/wiki/Bilinear_interpolation

  • https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-r

  • Separate from visualization - another round of interpolation!

    • https://matplotlib.org/3.3.3/gallery/images_contours_and_fields/interpolation_methods.html

• Raster interpolation from unstructured points (e.g. Lidar point clouds)

  • https://www.earthdatascience.org/courses/use-data-open-source-python/data-stories/what-is-lidar-data/lidar-points-to-pixels-raster/

• Combining raster and vector

  • sampling at points

  • zonal statistics for polygons

    • https://pysal.org/scipy2019-intermediate-gds/deterministic/gds2-rasters.html

• Volume calculations

• DEM derivative products: shaded relief, slope, aspect

• Raster/array filtering

  • Moving window operations

  • https://docs.scipy.org/doc/scipy/reference/ndimage.html

• vrt

⚠️ Suggestion - shut down other unneeded kernels before running

import os
import requests
import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio as rio
from rasterio import plot, mask
from rasterio.warp import calculate_default_transform, reproject, Resampling
import matplotlib.pyplot as plt
import rasterstats
import rioxarray as rxr
from matplotlib_scalebar.scalebar import ScaleBar
from pathlib import Path

#%matplotlib widget

Download sample DEM data for WA state

demdir = f'{Path.home()}/gda_demo_data/dem_data'

if not os.path.exists(demdir):
    os.makedirs(demdir)

!ls -lh $demdir

total 1.1G
-rw-r--r-- 1 jovyan users 143M Feb 15  2014 egm08_25.gtx
-rw-r--r-- 1 jovyan users 4.0M Feb 15  2014 egm96_15.gtx
-rw-rw-r-- 1 jovyan users 3.8M Feb 10 23:27 SEA_COP30.tif
-rw-rw-r-- 1 jovyan users 780K Feb 10 23:27 SEA_SRTMGL1.tif
-rw-r--r-- 1 jovyan users 149M Feb 12 19:34 WA_COP90_ellipsoid.tif
-rw-rw-r-- 1 jovyan users 142M Feb 10 23:27 WA_COP90.tif
-rw-r--r-- 1 jovyan users  50M Feb 12 18:04 WA_COP90_utm_gdalwarp_hs.tif
-rw-r--r-- 1 jovyan users  50M Feb 12 19:31 WA_COP90_utm_gdalwarp_lzw_hs.tif
-rw-rw-r-- 1 jovyan users 177M Feb 10 23:41 WA_COP90_utm_gdalwarp_lzw.tif
-rw-rw-r-- 1 jovyan users 199M Feb 10 23:39 WA_COP90_utm_gdalwarp.tif
-rw-r--r-- 1 jovyan users  75M Feb 12 19:33 WA_SRTMGL3_ellipsoid.tif
-rw-rw-r-- 1 jovyan users  46M Feb 10 23:27 WA_SRTMGL3.tif

Define the Washington state bounds from previous notebook

• Use the lat/lon bounds here in decimal degrees

wa_bounds = (-124.733174, 45.543541, -116.915989, 49.002494)

seattle_bounds = (-122.4487793446, 47.4510910088, -122.1548950672, 47.7917227803)

Use OpenTopography GlobalDEM API to fetch DEM for WA state

• OpenTopgraphy is a fantastic organization that “facilitates community access to high-resolution, Earth science-oriented, topography data, and related tools and resources.”

  • https://opentopography.org/about

• One of the many services they provide is an API for several popular Global DEM datasets, with simple subsetting and delivery: https://opentopography.org/developers

• We’ll use this service to extract a small portion of the SRTM-GL3 DEM

• See also their nice notebooks on bulk access and processing of cloud-optimized geotiffs (many familiar packages and code snippets): https://github.com/OpenTopography/OT_BulkAccess_COGs/blob/main/OT_BulkAccessCOGs.ipynb

#We are using OpenTopography demo API key here - only valid for 25 requests from a single IP address
#please save local output and use moving forward
#API_key="demoapikeyot2022"
#GDA key
API_key="e93e356203e00da7a5532e19ae4c2012"

base_url="https://portal.opentopography.org/API/globaldem?demtype={}&API_Key={}&west={}&south={}&east={}&north={}&outputFormat=GTiff"

def get_OT_GlobalDEM(demtype, bounds, out_fn=None):
    if out_fn is None:
        out_fn = '{}.tif'.format(demtype)
    
    if not os.path.exists(out_fn):
        print(f'Querying API for {demtype}')
        #Prepare API request url
        #Bounds should be [minlon, minlat, maxlon, maxlat]
        url = base_url.format(demtype, API_key, *bounds)
        print(url)
        #Get
        print(f'Downloading and Saving {out_fn}')
        response = requests.get(url)
        if response.status_code == 200:
            #Write to disk
            open(out_fn, 'wb').write(response.content)
        else:
            print(response.text)
    else:
        print(f'Found existing file: {out_fn}')

#Get Seattle data
for demtype in ['SRTMGL1', 'COP30']:
    
    #Output filename
    out_fn = f"{demdir}/SEA_{demtype}.tif"
    #Call the function to query API and download if file doesn't exist
    get_OT_GlobalDEM(demtype, seattle_bounds, out_fn)

Found existing file: /home/jovyan/gda_demo_data/dem_data/SEA_SRTMGL1.tif
Found existing file: /home/jovyan/gda_demo_data/dem_data/SEA_COP30.tif

#Get WA state data
for demtype in ['SRTMGL3', 'COP90']:
    
    #Output filename
    out_fn = f"{demdir}/WA_{demtype}.tif"
    #Call the function to query API and download if file doesn't exist
    get_OT_GlobalDEM(demtype, wa_bounds, out_fn)

Found existing file: /home/jovyan/gda_demo_data/dem_data/WA_SRTMGL3.tif
Found existing file: /home/jovyan/gda_demo_data/dem_data/WA_COP90.tif

!ls -lh $demdir

total 1.1G
-rw-r--r-- 1 jovyan users 143M Feb 15  2014 egm08_25.gtx
-rw-r--r-- 1 jovyan users 4.0M Feb 15  2014 egm96_15.gtx
-rw-rw-r-- 1 jovyan users 3.8M Feb 10 23:27 SEA_COP30.tif
-rw-rw-r-- 1 jovyan users 780K Feb 10 23:27 SEA_SRTMGL1.tif
-rw-r--r-- 1 jovyan users 149M Feb 12 19:34 WA_COP90_ellipsoid.tif
-rw-rw-r-- 1 jovyan users 142M Feb 10 23:27 WA_COP90.tif
-rw-r--r-- 1 jovyan users  50M Feb 12 18:04 WA_COP90_utm_gdalwarp_hs.tif
-rw-r--r-- 1 jovyan users  50M Feb 12 19:31 WA_COP90_utm_gdalwarp_lzw_hs.tif
-rw-rw-r-- 1 jovyan users 177M Feb 10 23:41 WA_COP90_utm_gdalwarp_lzw.tif
-rw-rw-r-- 1 jovyan users 199M Feb 10 23:39 WA_COP90_utm_gdalwarp.tif
-rw-r--r-- 1 jovyan users  75M Feb 12 19:33 WA_SRTMGL3_ellipsoid.tif
-rw-rw-r-- 1 jovyan users  46M Feb 10 23:27 WA_SRTMGL3.tif

!gdalinfo $out_fn

Driver: GTiff/GeoTIFF
Files: /home/jovyan/gda_demo_data/dem_data/WA_COP90.tif
Size is 9381, 4151
Coordinate System is:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        MEMBER["World Geodetic System 1984 (G2296)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-124.733750033333308,49.002916700000000)
Pixel Size = (0.000833333333333,-0.000833333333333)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  LAYOUT=COG
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-124.7337500,  49.0029167) (124d44' 1.50"W, 49d 0'10.50"N)
Lower Left  (-124.7337500,  45.5437500) (124d44' 1.50"W, 45d32'37.50"N)
Upper Right (-116.9162500,  49.0029167) (116d54'58.50"W, 49d 0'10.50"N)
Lower Right (-116.9162500,  45.5437500) (116d54'58.50"W, 45d32'37.50"N)
Center      (-120.8250000,  47.2733334) (120d49'30.00"W, 47d16'24.00"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray

Review the metadata

• Note the input data type - are these values unsigned or signed (meaning values can be positive or negative)?

Nice! WA state really does have interesting topography.

• Cascade mountains, Olympic Mountains, Active stratovolcanoes, Puget Sound, Columbia River

• Channeled Scablands: https://en.wikipedia.org/wiki/Channeled_Scablands

Much more interesting than Kansas:

• https://www.usu.edu/geo/geomorph/kansas.html

Compare SRTM And Copernicus DEM

Read the data and plot with matplotlib

srtm_wa_fn = f'{demdir}/WA_SRTMGL3.tif'
cop_wa_fn = f'{demdir}/WA_COP90.tif'

srtm_wa_da = rxr.open_rasterio(srtm_wa_fn, masked=True).squeeze()
cop_wa_da = rxr.open_rasterio(cop_wa_fn, masked=True).squeeze()
srtm_wa_da

<xarray.DataArray (y: 4151, x: 9381)> Size: 156MB
[38940531 values with dtype=float32]
Coordinates:
    band         int64 8B 1
  * x            (x) float64 75kB -124.7 -124.7 -124.7 ... -116.9 -116.9 -116.9
  * y            (y) float64 33kB 49.0 49.0 49.0 49.0 ... 45.55 45.55 45.54
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0

xarray.DataArray

• y: 4151
• x: 9381

• ...

  [38940531 values with dtype=float32]

• Coordinates: (4)
  • band
    ()
    int64
    1

    array(1)

  • x
    (x)
    float64
    -124.7 -124.7 ... -116.9 -116.9

    array([-124.733333, -124.7325  , -124.731667, ..., -116.918333, -116.9175  ,
           -116.916667])

  • y
    (y)
    float64
    49.0 49.0 49.0 ... 45.55 45.54

    array([49.0025  , 49.001667, 49.000833, ..., 45.545833, 45.545   , 45.544167])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    -124.73375000001255 0.000833333333333144 0.0 49.002916666669165 0.0 -0.000833333333333144

    array(0)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([-124.73333333334588, -124.73250000001255, -124.73166666667922,
           -124.73083333334588, -124.73000000001255, -124.72916666667922,
           -124.72833333334589, -124.72750000001255, -124.72666666667922,
           -124.72583333334589,
           ...
           -116.92416666668099, -116.92333333334766, -116.92250000001432,
           -116.92166666668099, -116.92083333334766, -116.92000000001433,
             -116.919166666681, -116.91833333334766, -116.91750000001433,
             -116.916666666681],
          dtype='float64', name='x', length=9381))

  • y
    PandasIndex

    PandasIndex(Index([  49.0025000000025, 49.001666666669166, 49.000833333335834,
             49.0000000000025,  48.99916666666917, 48.998333333335836,
             48.9975000000025, 48.996666666669164,  48.99583333333583,
             48.9950000000025,
           ...
            45.55166666666995,  45.55083333333662,  45.55000000000328,
            45.54916666666995, 45.548333333336615,  45.54750000000328,
            45.54666666666995,  45.54583333333662, 45.545000000003284,
            45.54416666666995],
          dtype='float64', name='y', length=4151))

• Attributes: (3)

  AREA_OR_POINT :

  Area

  scale_factor :

  1.0

  add_offset :

  0.0

f, axs = plt.subplots(1,2, sharex=True, sharey=True, figsize=(10,3))
srtm_wa_da.plot.imshow(ax=axs[0])
axs[0].set_title('SRTMGL3')
cop_wa_da.plot.imshow(ax=axs[1])
axs[1].set_title('COP90')

Text(0.5, 1.0, 'COP90')

Xarray defaults to the RdBu colormap when there are negative values! Let’s use a vmin and vmax to fix this.

f, axs = plt.subplots(1,2, sharex=True, sharey=True, figsize=(10,3))
srtm_wa_da.plot.imshow(ax=axs[0],vmin=0,vmax=4000)
axs[0].set_title('SRTMGL3')
cop_wa_da.plot.imshow(ax=axs[1],vmin=0,vmax=4000)
axs[1].set_title('COP90')

Text(0.5, 1.0, 'COP90')

But the COP90 DEM is showing lots of 0 values around water…

Set nodata value for Copernicus DEM to 0

• Note: The Copernicus DEM creators filled water with values of 0.0 above the EGM2008 geoid (approximates mean sea level)

• We can set these to nodata to only consider land. This is a bit of a hack, and it is possible that there are real elevations at 0 m that will be masked here.

print(f"nodata value for SRTMGL3: {srtm_wa_da.rio.nodata}")
print(f"nodata value for COP90: {cop_wa_da.rio.nodata}")

nodata value for SRTMGL3: nan
nodata value for COP90: None

cop_da = cop_wa_da.where(cop_wa_da != 0) # mask 0 values
cop_da = cop_da.rio.write_nodata(0) # update the attribute

print(f"nodata value for SRTMGL3: {srtm_wa_da.rio.nodata}")
print(f"nodata value for COP90: {cop_da.rio.nodata}")

nodata value for SRTMGL3: nan
nodata value for COP90: 0.0

f, axs = plt.subplots(1,2, sharex=True, sharey=True, figsize=(10,3))
srtm_wa_da.plot.imshow(ax=axs[0],vmin=0,vmax=4000)
axs[0].set_title('SRTMGL3')
cop_da.plot.imshow(ax=axs[1],vmin=0,vmax=4000)
axs[1].set_title('COP90')

Text(0.5, 1.0, 'COP90')

Raster reprojection

So easy to do in rioxarray!

Reproject the Copernicus DEM to UTM 10N

• The COP grids are distributed with crs of EPSG:4326

• Let’s reproject to a more appropriate coordinate system for WA state

  • Use EPSG:32610

dst_crs = 'EPSG:32610'

cop_wa_utm_da = cop_wa_da.rio.reproject(dst_crs)
cop_wa_utm_da

<xarray.DataArray (y: 5851, x: 8877)> Size: 208MB
array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)
Coordinates:
  * x            (x) float64 71kB 3.647e+05 3.648e+05 ... 9.748e+05 9.749e+05
  * y            (y) float64 47kB 5.446e+06 5.446e+06 ... 5.043e+06 5.043e+06
    band         int64 8B 1
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:  Point
    scale_factor:   1.0
    add_offset:     0.0
    _FillValue:     nan

xarray.DataArray

• y: 5851
• x: 8877

• nan nan nan nan nan nan nan nan ... nan nan nan nan nan nan nan nan

  array([[nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan],
         ...,
         [nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)

• Coordinates: (4)
  • x
    (x)
    float64
    3.647e+05 3.648e+05 ... 9.749e+05

    axis :

    X

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    units :

    metre

    array([364687.335257, 364756.083719, 364824.832181, ..., 974761.185449,
           974829.933911, 974898.682373])

  • y
    (y)
    float64
    5.446e+06 5.446e+06 ... 5.043e+06

    axis :

    Y

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    units :

    metre

    array([5445601.595528, 5445532.847066, 5445464.098604, ..., 5043560.590804,
           5043491.842342, 5043423.09388 ])

  • band
    ()
    int64
    1

    array(1)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 10N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -123.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    GeoTransform :

    364652.9610264343 68.74846182011875 0.0 5445635.969758815 0.0 -68.74846182011875

    array(0)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([364687.33525734436,  364756.0837191645,  364824.8321809846,
           364893.58064280474,  364962.3291046248, 365031.07756644493,
           365099.82602826506,  365168.5744900852,  365237.3229519053,
           365306.07141372544,
           ...
            974279.9462163374,  974348.6946781573,  974417.4431399775,
            974486.1916017977,  974554.9400636179,  974623.6885254378,
             974692.436987258,  974761.1854490782,  974829.9339108982,
            974898.6823727184],
          dtype='float64', name='x', length=8877))

  • y
    PandasIndex

    PandasIndex(Index([ 5445601.595527905,  5445532.847066085,  5445464.098604265,
           5445395.3501424445,  5445326.601680624,  5445257.853218804,
            5445189.104756984,  5445120.356295164, 5445051.6078333445,
            5444982.859371524,
           ...
            5044041.830036592,  5043973.081574772,  5043904.333112951,
            5043835.584651131,  5043766.836189311,  5043698.087727491,
            5043629.339265671,   5043560.59080385,   5043491.84234203,
             5043423.09388021],
          dtype='float64', name='y', length=5851))

• Attributes: (4)

  AREA_OR_POINT :

  Point

  scale_factor :

  1.0

  add_offset :

  0.0

  _FillValue :

  nan

Let’s also try gdalwarp command-line utility

• A very simple, efficient way to accomplish this - let GDAL worry about all of the underlying math

• Review the documentation and options here: https://gdal.org/programs/gdalwarp.html

  • Resampling options: https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-r

• Save the projected file as a GeoTiff on disk

  • Use cubic resampling algorithm

• Note the filesize using ls -lh

proj_fn = os.path.splitext(out_fn)[0]+'_utm_gdalwarp.tif'

%%time
if not os.path.exists(proj_fn):
    !gdalwarp -overwrite -r cubic -t_srs $dst_crs $out_fn $proj_fn

CPU times: user 25 µs, sys: 29 µs, total: 54 µs
Wall time: 48.2 µs

!ls -lh $proj_fn

-rw-rw-r-- 1 jovyan users 199M Feb 10 23:39 /home/jovyan/gda_demo_data/dem_data/WA_COP90_utm_gdalwarp.tif

!gdalinfo $proj_fn

Driver: GTiff/GeoTIFF
Files: /home/jovyan/gda_demo_data/dem_data/WA_COP90_utm_gdalwarp.tif
Size is 8877, 5851
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 10N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 10N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-123,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK)."],
        BBOX[0,-126,84,-120]],
    ID["EPSG",32610]]
Data axis to CRS axis mapping: 1,2
Origin = (364652.961026434320956,5445635.969758815132082)
Pixel Size = (68.748461820118749,-68.748461820118749)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  364652.961, 5445635.970) (124d51'21.71"W, 49d 8'55.04"N)
Lower Left  (  364652.961, 5043388.720) (124d44' 0.08"W, 45d31'51.18"N)
Upper Right (  974933.057, 5445635.970) (116d30'26.42"W, 48d58'49.95"N)
Lower Right (  974933.057, 5043388.720) (116d56' 0.37"W, 45d22'57.49"N)
Center      (  669793.009, 5244512.345) (120d45' 9.46"W, 47d19'55.33"N)
Band 1 Block=8877x1 Type=Float32, ColorInterp=Gray

A note about “creation options” when writing a new file to disk

• Review this great reference on GeoTiff file format: https://www.gdal.org/frmt_gtiff.html

  • See Creation Options section toward bottom of page: https://gdal.org/drivers/raster/gtiff.html#creation-options

• I almost always use -co COMPRESS=LZW -co TILED=YES -co BIGTIFF=IF_SAFER

  • Uses lossless LZW compression (we explored in LS8 example comparing filesize on disk to your calculated array filesize)

  • Writes the tif image in “tiles” of 256x256 px instead of one large block - makes it much more efficient to read and extract a subwindow from the tif

  • If necessary, use the BIGTIFF format for files that are >4 GB

Write out the same file, but this time use tiling and LZW compression

• Compare the new filesize on disk

proj_fn_lzw = os.path.splitext(out_fn)[0]+'_utm_gdalwarp_lzw.tif'

%%time
if not os.path.exists(proj_fn_lzw):
    !gdalwarp -r cubic -t_srs $dst_crs -co COMPRESS=LZW -co TILED=YES -co BIGTIFF=IF_SAFER $out_fn $proj_fn_lzw

CPU times: user 29 µs, sys: 33 µs, total: 62 µs
Wall time: 55.1 µs

!ls -lh $proj_fn_lzw

-rw-rw-r-- 1 jovyan users 177M Feb 10 23:41 /home/jovyan/gda_demo_data/dem_data/WA_COP90_utm_gdalwarp_lzw.tif

#Let's use the compressed file from here on out
#Should compare runtimes when using compressed vs. uncompressed
proj_fn = proj_fn_lzw

!gdalinfo $proj_fn

Driver: GTiff/GeoTIFF
Files: /home/jovyan/gda_demo_data/dem_data/WA_COP90_utm_gdalwarp_lzw.tif
Size is 8877, 5851
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 10N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 10N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-123,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK)."],
        BBOX[0,-126,84,-120]],
    ID["EPSG",32610]]
Data axis to CRS axis mapping: 1,2
Origin = (364652.961026434320956,5445635.969758815132082)
Pixel Size = (68.748461820118749,-68.748461820118749)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  364652.961, 5445635.970) (124d51'21.71"W, 49d 8'55.04"N)
Lower Left  (  364652.961, 5043388.720) (124d44' 0.08"W, 45d31'51.18"N)
Upper Right (  974933.057, 5445635.970) (116d30'26.42"W, 48d58'49.95"N)
Lower Right (  974933.057, 5043388.720) (116d56' 0.37"W, 45d22'57.49"N)
Center      (  669793.009, 5244512.345) (120d45' 9.46"W, 47d19'55.33"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray

Discussion Questions

• What is the original x and y pixel size and units?

  • Remember this is a 3-arcsecond grid

  • Info on an arcsecond: https://www.esri.com/news/arcuser/0400/wdside.html

• What is the projected x and y pixel size and units?

  • Is this consistent with your expectation for the 90-m Copernicus DEM or 3-arcsec SRTM data?

    • Hint: Think about the dimensions of a pixel in meters near the top and bottom of the unprojected raster.

      • Roughly 45.5° and 49°

    • Then think back to Lab03 when we calculated length of a degree of longitude and a degree of latitude at different locations on the planet.

print(f"the original x and y pixel size: {cop_wa_da.rio.resolution()}")
print(f"the reprojected x and y pixel size: {cop_wa_utm_da.rio.resolution()}")

the original x and y pixel size: (0.0008333333333333334, -0.0008333333333333334)
the reprojected x and y pixel size: (68.74846182011875, -68.74846182011875)

• Note that if the output resolution is unspecified, GDAL and rasterio will calcluate for you.

  • More info on the default “Suggested Warp Output” resolution: https://gdal.org/api/gdal_alg.html#_CPPv423GDALSuggestedWarpOutput12GDALDatasetH19GDALTransformerFuncPvPdPiPi

…resolution is computed with the intent that the length of the distance from the top left corner of the output imagery to the bottom right corner would represent the same number of pixels as in the source image. Note that if the image is somewhat rotated the diagonal taken isn’t of the whole output bounding rectangle, but instead of the locations where the top/left and bottom/right corners transform. The output pixel size is always square. This is intended to approximately preserve the resolution of the input data in the output file.

• But you can also specify this output resolution (-tr argument in gdalwarp), if, for example, you wanted to resample your output to 180 m resolution.

Part 2: Create a shaded relief map (hillshade) from your projected DEM

• See background info on hillshades (and other DEM derivative products like slope and aspect) here: http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-hillshade-works.htm

• Can do this easily with gdaldem command line utility: https://gdal.org/programs/gdaldem.html

  • Try this and output to a new file

  • Note that if you were to run on a DEM with different horizontal and vertical units, (like say, degrees and meters), you will have scaling issues. There are ways around this, but always generate hillshade using a projected DEM.

hs_fn = os.path.splitext(proj_fn)[0]+'_hs.tif'

%%time
if not os.path.exists(hs_fn):
    !gdaldem hillshade $proj_fn $hs_fn

CPU times: user 19 µs, sys: 21 µs, total: 40 µs
Wall time: 47.9 µs

Use rasterio to load your GDAL created hillshade

• Note that the datatype is 8-bit (Byte), with values from 0-255 (it’s a grayscale image, not elevation data)

• Make sure you mask nodata values

hs_wa_da = rxr.open_rasterio(hs_fn, masked=True).squeeze()
hs_wa_da

<xarray.DataArray (y: 5851, x: 8877)> Size: 208MB
[51939327 values with dtype=float32]
Coordinates:
    band         int64 8B 1
  * x            (x) float64 71kB 3.647e+05 3.648e+05 ... 9.748e+05 9.749e+05
  * y            (y) float64 47kB 5.446e+06 5.446e+06 ... 5.043e+06 5.043e+06
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0

xarray.DataArray

• y: 5851
• x: 8877

• ...

  [51939327 values with dtype=float32]

• Coordinates: (4)
  • band
    ()
    int64
    1

    array(1)

  • x
    (x)
    float64
    3.647e+05 3.648e+05 ... 9.749e+05

    array([364687.335257, 364756.083719, 364824.832181, ..., 974761.185449,
           974829.933911, 974898.682373])

  • y
    (y)
    float64
    5.446e+06 5.446e+06 ... 5.043e+06

    array([5445601.595528, 5445532.847066, 5445464.098604, ..., 5043560.590804,
           5043491.842342, 5043423.09388 ])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 10N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -123.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    GeoTransform :

    364652.9610264343 68.74846182011875 0.0 5445635.969758815 0.0 -68.74846182011875

    array(0)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([364687.33525734436,  364756.0837191645,  364824.8321809846,
           364893.58064280474,  364962.3291046248, 365031.07756644493,
           365099.82602826506,  365168.5744900852,  365237.3229519053,
           365306.07141372544,
           ...
            974279.9462163374,  974348.6946781573,  974417.4431399775,
            974486.1916017977,  974554.9400636179,  974623.6885254378,
             974692.436987258,  974761.1854490782,  974829.9339108982,
            974898.6823727184],
          dtype='float64', name='x', length=8877))

  • y
    PandasIndex

    PandasIndex(Index([ 5445601.595527905,  5445532.847066085,  5445464.098604265,
           5445395.3501424445,  5445326.601680624,  5445257.853218804,
            5445189.104756984,  5445120.356295164, 5445051.6078333445,
            5444982.859371524,
           ...
            5044041.830036592,  5043973.081574772,  5043904.333112951,
            5043835.584651131,  5043766.836189311,  5043698.087727491,
            5043629.339265671,   5043560.59080385,   5043491.84234203,
             5043423.09388021],
          dtype='float64', name='y', length=5851))

• Attributes: (3)

  AREA_OR_POINT :

  Area

  scale_factor :

  1.0

  add_offset :

  0.0

f,ax=plt.subplots(figsize=(10, 5))
hs_wa_da.plot.imshow(ax=ax,vmin=0,vmax=360,cmap='gray', interpolation='bilinear')
ax.set_aspect('equal')

Part 3: Clipping Raster data

Load the states GeoDataFrame

#states_url = 'http://eric.clst.org/assets/wiki/uploads/Stuff/gz_2010_us_040_00_5m.json'
states_url = 'http://eric.clst.org/assets/wiki/uploads/Stuff/gz_2010_us_040_00_500k.json'
states_gdf = gpd.read_file(states_url)
states_gdf.head()

	GEO_ID	STATE	NAME	LSAD	CENSUSAREA	geometry
0	0400000US23	23	Maine		30842.923	MULTIPOLYGON (((-67.61976 44.51975, -67.61541 ...
1	0400000US25	25	Massachusetts		7800.058	MULTIPOLYGON (((-70.83204 41.6065, -70.82374 4...
2	0400000US26	26	Michigan		56538.901	MULTIPOLYGON (((-88.68443 48.11578, -88.67563 ...
3	0400000US30	30	Montana		145545.801	POLYGON ((-104.0577 44.99743, -104.25014 44.99...
4	0400000US32	32	Nevada		109781.180	POLYGON ((-114.0506 37.0004, -114.05 36.95777,...

Reproject the GeoDataFrame to match your DEM

states_gdf_proj = states_gdf.to_crs(dst_crs)

Isolate the WA state geometry object

• We want the geometry, not a GeoDataFrame or GeoSeries

wa_state = states_gdf_proj.loc[states_gdf_proj['NAME'] == 'Washington']
wa_geom = wa_state.iloc[0].geometry

wa_geom

f,axs = plt.subplots(1,2, figsize=(12,4))
cop_wa_utm_da.plot.imshow(ax=axs[0],vmin=0,vmax=4000)
axs[0].set_title('COP')
axs[0].set_aspect('equal')
cop_wa_utm_da.rio.clip(wa_state.geometry).plot.imshow(ax=axs[1],vmin=0,vmax=4000)
axs[1].set_title('cop clipped to WA')
axs[1].set_aspect('equal')

Part 4: Comparing and differencing rasters

• Often have two rasters with different resolution, extent, projection

• Want to combine, or use for some kind of raster math (e.g., subtract one from the other to compute difference)

• Relatively new package, rioxarray makes this very easy

  • https://corteva.github.io/rioxarray/stable/index.html

  • More on this in Week09 and Week10

  • Built for NetCDF model, but wraps rasterio for working with rasters

  • Lazy evaluation, doesn’t read from disk util necessary

  • Built-in support for parallel operations using Dask

#Original SRTM DEM, still in EPSG:4326
print(f"Original SRTM DEM: {srtm_wa_da.rio.crs}")
print(f"Original SRTM dimensions: {srtm_wa_da.rio.shape}")
print(f"Reprojected COP DEM: {cop_wa_utm_da.rio.crs}")
print(f"Reprojected COP dimensions: {cop_wa_utm_da.rio.shape}")

Original SRTM DEM: EPSG:4326
Original SRTM dimensions: (4151, 9381)
Reprojected COP DEM: EPSG:32610
Reprojected COP dimensions: (5851, 8877)

#Projected COP DEM, reprojected in EPSG:32610
srtm_wa_utm_da = srtm_wa_da.rio.reproject_match(cop_wa_utm_da)

print(f"Reprojected SRTM DEM: {srtm_wa_utm_da.rio.crs}")
print(f"Reprojected SRTM dimensions: {srtm_wa_utm_da.rio.shape}")
print(f"Reprojected COP DEM: {cop_wa_utm_da.rio.crs}")
print(f"Reprojected COP dimensions: {cop_wa_utm_da.rio.shape}")

Reprojected SRTM DEM: EPSG:32610
Reprojected SRTM dimensions: (5851, 8877)
Reprojected COP DEM: EPSG:32610
Reprojected COP dimensions: (5851, 8877)

Create a difference map

• Now that both rasters have identical projection, extent, and resolution, this is a simple subtraction

diff_da = cop_wa_utm_da - srtm_wa_utm_da

f, ax = plt.subplots(dpi=150)
diff_da.plot.imshow(cmap='RdBu', robust=True, ax=ax)
ax.set_aspect('equal')
ax.set_title('COP90 minus SRTMGL3')

Text(0.5, 1.0, 'COP90 minus SRTMGL3')

f, ax = plt.subplots(dpi=150)
diff_da.plot.imshow(cmap='RdBu', robust=True, vmin=-5, vmax=5, ax=ax)
ax.set_aspect('equal')
ax.set_title('COP90 minus SRTMGL3')

Text(0.5, 1.0, 'COP90 minus SRTMGL3')

diff_da.median()

<xarray.DataArray ()> Size: 4B
array(1.5637817, dtype=float32)
Coordinates:
    band         int64 8B 1
    spatial_ref  int64 8B 0

xarray.DataArray

• 1.564

  array(1.5637817, dtype=float32)

• Coordinates: (2)
  • band
    ()
    int64
    1

    array(1)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 10N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -123.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    GeoTransform :

    364652.9610264343 68.74846182011875 0.0 5445635.969758815 0.0 -68.74846182011875

    array(0)

• Indexes: (0)
• Attributes: (0)

diff_da.mean()

<xarray.DataArray ()> Size: 4B
array(1.1240803, dtype=float32)
Coordinates:
    band         int64 8B 1
    spatial_ref  int64 8B 0

xarray.DataArray

• 1.124

  array(1.1240803, dtype=float32)

• Coordinates: (2)
  • band
    ()
    int64
    1

    array(1)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 10N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -123.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    GeoTransform :

    364652.9610264343 68.74846182011875 0.0 5445635.969758815 0.0 -68.74846182011875

    array(0)

• Indexes: (0)
• Attributes: (0)

diff_da.std()

<xarray.DataArray ()> Size: 4B
array(6.0349627, dtype=float32)
Coordinates:
    band         int64 8B 1
    spatial_ref  int64 8B 0

xarray.DataArray

• 6.035

  array(6.0349627, dtype=float32)

• Coordinates: (2)
  • band
    ()
    int64
    1

    array(1)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 10N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -123.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    GeoTransform :

    364652.9610264343 68.74846182011875 0.0 5445635.969758815 0.0 -68.74846182011875

    array(0)

• Indexes: (0)
• Attributes: (0)

What about the vertical datum?

• Remember that these DEMs contain elevation values that are relative to some reference surface, the vertical datum!

• The SRTM DEM uses the EGM96 geoid, while the Copernicus DEM uses the EGM2008 geoid. These are slightly different. https://en.wikipedia.org/wiki/Earth_Gravitational_Model

• So before we difference these DEMs, let’s change the elevation values to be relative to the WGS84 ellipsoid.

• We can use gdalwarp for this, but we need to download files that contain the offset between the geoid and the ellipsoid. https://proj.org/en/stable/usage/transformation.html

# files containing vertical datum shifts
egm96 = '/home/jovyan/gda_demo_data/dem_data/egm96_15.gtx'
egm08 = '/home/jovyan/gda_demo_data/dem_data/egm08_25.gtx'

# download the files
if not os.path.exists(egm96):
    !wget -P $demdir https://download.osgeo.org/proj/vdatum/egm96_15/egm96_15.gtx
if not os.path.exists(egm08):
    !wget -P $demdir https://download.osgeo.org/proj/vdatum/egm08_25/egm08_25.gtx

# in files
cop90_wa_EGM08_fn = '/home/jovyan/gda_demo_data/dem_data/WA_COP90.tif'
srtm_wa_EGM96_fn = '/home/jovyan/gda_demo_data/dem_data/WA_SRTMGL3.tif'

# out files
cop90_wa_wgs84_ellipsoid_fn = '/home/jovyan/gda_demo_data/dem_data/WA_COP90_ellipsoid.tif'
srtm_wa_wgs84_ellipsoid_fn = '/home/jovyan/gda_demo_data/dem_data/WA_SRTMGL3_ellipsoid.tif'

# transform the vertical datum of our SRTM DEM
if not os.path.exists(srtm_wa_wgs84_ellipsoid_fn):
    !gdalwarp -s_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=/home/jovyan/gda_demo_data/dem_data/egm96_15.gtx" -t_srs "+proj=longlat +datum=WGS84 +no_def" $srtm_wa_EGM96_fn $srtm_wa_wgs84_ellipsoid_fn

# transform the vertical datum of our SRTM DEM
if not os.path.exists(cop90_wa_wgs84_ellipsoid_fn):
    !gdalwarp -s_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=/home/jovyan/gda_demo_data/dem_data/egm08_25.gtx" -t_srs "+proj=longlat +datum=WGS84 +no_def" $cop90_wa_EGM08_fn $cop90_wa_wgs84_ellipsoid_fn

# read DEMs back in
srtm_wa_da = rxr.open_rasterio(srtm_wa_wgs84_ellipsoid_fn, masked=True).squeeze()
cop_wa_da = rxr.open_rasterio(cop90_wa_wgs84_ellipsoid_fn, masked=True).squeeze()

srtm_wa_utm_da = srtm_wa_da.rio.reproject(dst_crs)

cop_wa_utm_da = cop_wa_da.rio.reproject_match(srtm_wa_utm_da)

diff_da = cop_wa_utm_da - srtm_wa_utm_da

f, ax = plt.subplots(dpi=150)
diff_da.plot.imshow(cmap='RdBu', robust=True, vmin=-5, vmax=5, ax=ax)
ax.set_aspect('equal')
ax.set_title('COP90 minus SRTMGL3')

Text(0.5, 1.0, 'COP90 minus SRTMGL3')

diff_da.median()

<xarray.DataArray ()> Size: 4B
array(1.2817993, dtype=float32)
Coordinates:
    band         int64 8B 1
    spatial_ref  int64 8B 0

xarray.DataArray

• 1.282

  array(1.2817993, dtype=float32)

• Coordinates: (2)
  • band
    ()
    int64
    1

    array(1)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 10N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -123.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    GeoTransform :

    364652.9635476386 68.74846185120099 0.0 5445635.966252405 0.0 -68.74846185120093

    array(0)

• Indexes: (0)
• Attributes: (0)

diff_da.mean()

<xarray.DataArray ()> Size: 4B
array(1.0107714, dtype=float32)
Coordinates:
    band         int64 8B 1
    spatial_ref  int64 8B 0

xarray.DataArray

• 1.011

  array(1.0107714, dtype=float32)

• Coordinates: (2)
  • band
    ()
    int64
    1

    array(1)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 10N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -123.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    GeoTransform :

    364652.9635476386 68.74846185120099 0.0 5445635.966252405 0.0 -68.74846185120093

    array(0)

• Indexes: (0)
• Attributes: (0)

diff_da.std()

<xarray.DataArray ()> Size: 4B
array(6.018347, dtype=float32)
Coordinates:
    band         int64 8B 1
    spatial_ref  int64 8B 0

xarray.DataArray

• 6.018

  array(6.018347, dtype=float32)

• Coordinates: (2)
  • band
    ()
    int64
    1

    array(1)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 10N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -123.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    GeoTransform :

    364652.9635476386 68.74846185120099 0.0 5445635.966252405 0.0 -68.74846185120093

    array(0)

• Indexes: (0)
• Attributes: (0)

# plot change at Mt. Baker
f, ax = plt.subplots(1, 3, figsize=(14, 4), sharex=True, sharey=True)
cop_wa_utm_da.isel(x=slice(3120, 3320), y=slice(510, 710)).plot(ax=ax[0])
hs_wa_da.isel(x=slice(3120, 3320), y=slice(510, 710)).plot(ax=ax[1], cmap='Grays_r')
diff_da.isel(x=slice(3120, 3320), y=slice(510, 710)).plot(ax=ax[2], vmin=-30, vmax=30, cmap='RdBu')
ax[0].set_aspect('equal')
ax[1].set_aspect('equal')
ax[2].set_aspect('equal')
plt.tight_layout()

Save output raster to disk

diff_fn = 'WA_COP90_SRTMGL3_diff.tif'

diff_da.rio.to_raster(diff_fn)

!gdalinfo $diff_fn

Driver: GTiff/GeoTIFF
Files: WA_COP90_SRTMGL3_diff.tif
Size is 8877, 5851
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 10N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 10N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-123,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK)."],
        BBOX[0,-126,84,-120]],
    ID["EPSG",32610]]
Data axis to CRS axis mapping: 1,2
Origin = (364652.963547638617456,5445635.966252405196428)
Pixel Size = (68.748461851200986,-68.748461851200929)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  364652.964, 5445635.966) (124d51'21.71"W, 49d 8'55.04"N)
Lower Left  (  364652.964, 5043388.716) (124d44' 0.08"W, 45d31'51.18"N)
Upper Right (  974933.059, 5445635.966) (116d30'26.42"W, 48d58'49.94"N)
Lower Right (  974933.059, 5043388.716) (116d56' 0.37"W, 45d22'57.49"N)
Center      (  669793.011, 5244512.341) (120d45' 9.46"W, 47d19'55.33"N)
Band 1 Block=8877x1 Type=Float32, ColorInterp=Gray

Profile plot

#y=5300000
y=5190000
f, ax = plt.subplots()
cop_wa_utm_da.plot.imshow(cmap='inferno', robust=True, ax=ax)
ax.set_aspect('equal')
ax.axhline(y)

f, axs = plt.subplots(2, 1)
cop_wa_utm_da.sel(y=y, method='nearest').plot(ax=axs[0], label='COP90')
srtm_wa_utm_da.sel(y=y, method='nearest').plot(ax=axs[0], label='SRTM_GL1')
axs[0].legend()
axs[0].set_ylabel('Elevation (m)')
diff_da.sel(y=y, method='nearest').plot(ax=axs[1])
axs[1].axhline(0, color='k')
axs[1].set_ylabel('Difference (m)')
plt.tight_layout()

Some other newer packages worth exploring

Work in Progress: xarray-spatial

• Newer project implementing many raster/DEM analysis functions in Python with Dask/CUDA support

• Uses xarray data model (more complicated than rasterio dataset model)

• https://xarray-spatial.org/

• Surface functions: https://github.com/makepath/xarray-spatial#surface

  • Hillshade: https://xarray-spatial.org/reference/_autosummary/xrspatial.hillshade.hillshade.html#xrspatial.hillshade.hillshade

#import xrspatial

#Run the hillshade operation
#hs = xrspatial.hillshade(cop_xds, azimuth=315, angle_altitude=45)
#Note: there are some coordinate scaling and nodata issues here
#hs.plot.imshow(cmap='gray');

xdem

• https://xdem.readthedocs.io/en/latest/quick_start.html