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
```

`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]: