Merge pull request SciTools#396 from pelson/srtm_shading

Elevation filling & shading by @ThomasLecocq.

docs/source/contributors.rst

@@ -24,6 +24,7 @@ the package wouldn't be as rich or diverse as it is today:
  * Andrew Dawson
  * Patrick Peglar
  * Christoph Gohlke
+ * Thomas Lecocq
 
 
 Thank you!

docs/source/whats_new.rst

@@ -1,3 +1,14 @@
+What's new in cartopy 0.11
+**************************
+
+Thomas Lecocq added functionality to :mod:`cartopy.io.srtm` allowing intelligent
+filling of missing elevation data, as well as a function to compute elevation
+shading for relief style mapping. An example has been added which uses both of
+these functions to produce a grayscale shaded relief map:
+
+.. plot:: examples/srtm_shading.py
+   :width: 300pt
+
 What's new in cartopy 0.10
 **************************
 
@@ -31,7 +42,6 @@ capabilities, including:
 .. plot:: examples/barbs.py
     :width: 300pt
 
-
 What's new in cartopy 0.9
 *************************
 

lib/cartopy/examples/srtm_shading.py

@@ -0,0 +1,42 @@
+__tags__ = ['Scalar data']
+"""
+This example illustrates the automatic download of
+STRM data, gap filling (using gdal) and adding shading
+to create a so-called "Shaded Relief SRTM"
+
+Contributed by: Thomas Lecocq (http://geophysique.be)
+"""
+
+import cartopy.crs as ccrs
+from cartopy.io import srtm
+import matplotlib.pyplot as plt
+
+
+def main():
+    ax = plt.axes(projection=ccrs.PlateCarree())
+
+    # Get the 1x1 degree SRTM tile for 12E, 47N
+    elev, crs, extent = srtm.srtm_composite(12, 47, 1, 1)
+
+    # Fill the gaps present in the elevation data
+    elev_filled = srtm.fill_gaps(elev, 15)
+
+    # Add shading simulating the Sun at 10am (South-East)
+    # and with a low angle (15 degrees above horizon)
+    shaded = srtm.add_shading(elev_filled, 135.0, 15.0)
+
+    # The plot the result :
+    plt.imshow(shaded, extent=extent, transform=crs,
+               cmap='Greys', origin='lower')
+
+    plt.title("SRTM Shaded Relief Map")
+
+    gl = ax.gridlines(draw_labels=True,)
+    gl.xlabels_top = False
+    gl.ylabels_left = False
+
+    plt.show()
+
+
+if __name__ == '__main__':
+    main()

lib/cartopy/io/srtm.py

@@ -38,13 +38,72 @@
 def srtm(lon, lat):
     """
     Return (elevation, crs, extent) for the given longitude latitude.
+    Elevation is in meters.
     """
     fname = SRTM3_retrieve(lon, lat)
     if fname is None:
         raise ValueError('No srtm tile found for those coordinates.')
     return read_SRTM3(fname)
 
 
+def add_shading(elevation, azimuth, altitude):
+    """Adds shading to SRTM elevation data, using azimuth and altitude
+    of the sun.
+
+    :type elevation: numpy.ndarray
+    :param elevation: SRTM elevation data (in meters)
+    :type azimuth: float
+    :param azimuth: azimuth of the Sun (in degrees)
+    :type altitude: float
+    :param altitude: altitude of the Sun (in degrees)
+
+    :rtype: numpy.ndarray
+    :return: shaded SRTM relief map.
+    """
+    azimuth = np.deg2rad(azimuth)
+    altitude = np.deg2rad(altitude)
+    x, y = np.gradient(elevation)
+    slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))
+    # -x here 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
+
+
+def fill_gaps(elevation, max_distance=10):
+    """Fills gaps in SRTM elevation data for which the distance from
+    missing pixel to nearest existing one is smaller than `max_distance`.
+
+    This function requires osgeo/gdal to work.
+
+    :type elevation: numpy.ndarray
+    :param elevation: SRTM elevation data (in meters)
+    :type max_distance: int
+    :param max_distance: maximal distance (in pixels) between a missing point
+    and the nearest valid one.
+
+    :rtype: numpy.ndarray
+    :return: SRTM elevation data with filled gaps..
+    """
+    # Lazily import osgeo - it is only an optional dependency for cartopy.
+    from osgeo import gdal
+    from osgeo import gdal_array
+
+    src_ds = gdal_array.OpenArray(elevation)
+    srcband = src_ds.GetRasterBand(1)
+    dstband = srcband
+    maskband = srcband
+    smoothing_iterations = 0
+    options = []
+    gdal.FillNodata(dstband, maskband,
+                    max_distance, smoothing_iterations, options,
+                    callback=None)
+    elevation = dstband.ReadAsArray()
+    return elevation
+
+
 def srtm_composite(lon_min, lat_min, nx, ny):
 
     # XXX nx and ny have got confused in the code (numpy array ordering?).