python4oceanographers

Turning ripples into waves

Plotting shaded relief

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.

Creative Commons License
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/.

Comments