python4oceanographers

Turning ripples into waves

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:

• 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]:
"""Entre altitude and azimuth in degress."""
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))

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)
smoothing_iterations = 0
options = []
max_distance = 200
smoothing_iterations, options,
callback=None)
In [5]:

fig, ax = plt.subplots(figsize=(6, 6),
subplot_kw=dict(projection=ccrs.PlateCarree()))
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.