This is a quick post to test Cartopy's shaded relief plots. I never had to plot something like this before, but recently a friend asked me how to this with python and I got interested.

After a quick Google search I found several blogs post about this. The best one IMHO was from my favorite geo-blog: http://www.geophysique.be/. Here is a summary of what I learner:

Cartopy can download srtm composite for you.
GDAL has all the tools needed to read, smooth, and fill data gaps.
Before plotting play with the altitude and the azimuth to find the desired shading.

So let's start,

In [2]:

import numpy as np
import matplotlib.pyplot as plt

import cartopy.crs as ccrs
from cartopy.io.srtm import srtm_composite

from osgeo import gdal
from osgeo import gdal_array

The hill_shade() function is just a compilation of the steps I found here.

In [3]:

def hill_shade(elev, altitude=45, azimuth=315):
    """Entre altitude and azimuth in degress."""
    altitude = np.deg2rad(altitude)
    azimuth = np.deg2rad(azimuth)
    x, y = np.gradient(elev)
    slope = np.pi/2. - np.arctan(np.sqrt(x**2 + y**2))
    # Negative x because of pixel orders in the SRTM tile.
    aspect = np.arctan2(-x, y)

    shaded = (np.sin(altitude) * np.sin(slope) +
              np.cos(altitude) * np.cos(slope) *
              np.cos((azimuth - np.pi/2.) - aspect))
    return shaded

Cartopy's srtm_composite() downloads the data from http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/ auto-magically. The functions takes a starting lon, lat position and a integer for the step.

In [4]:

lon, dx = -42.8, 1
lat, dy = -10.8, 1

elev, crs, extent = srtm_composite(lon, lat, dx, dy)
src_ds = gdal_array.OpenArray(elev)
srcband = src_ds.GetRasterBand(1)
dstband = maskband = srcband
smoothing_iterations = 0
options = []
max_distance = 200
result = gdal.FillNodata(dstband, maskband, max_distance,
                         smoothing_iterations, options,
                         callback=None)
elev = dstband.ReadAsArray()

In [5]:

shaded = hill_shade(elev, altitude=45, azimuth=90)

fig, ax = plt.subplots(figsize=(6, 6),
                       subplot_kw=dict(projection=ccrs.PlateCarree()))
ax.imshow(shaded, extent=extent, transform=crs,
          cmap='Greys', origin='lower')

ax.set_title(u"São Francisco River")
gl = ax.gridlines(draw_labels=True,)
gl.xlabels_top = gl.ylabels_left = False

That plot show my next vacation hiking trip!

In [6]:

HTML(html)

Out[6]:

This post was written as an IPython notebook. It is available for download or as a static html.

python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.

python4oceanographers

Turning ripples into waves

Plotting shaded relief

Comments