I used gpxpy in a previous post to read and plot GPX files from smart-phones. Here I will show another module from the same author, Tomo Krajina, named srtm.py. This module reads SRTM files and extract the elevation data at the positions recorded in a GPX file. That is very useful! I never really trusted my cell phone GPS as an altimeter ;-)

First let's get an SRMT image centered at our GPS file.

In [10]:

import srtm
import gpxpy
import numpy as np
import numpy.ma as ma

fname = './data/2014_07_26_pai_ignacio.gpx'
gpx = gpxpy.parse(open(fname))

dx = dy = 3
s = gpx.tracks[0].segments[0]
lon, lat = np.array([(p.longitude, p.latitude) for p in s.points]).mean(axis=0)

data = srtm.get_data()
img = data.get_image((250, 250),
                     (lat - dy/2., lat + dy),
                     (lon - dx/2., lon + dx), 1000)

# Convert to numpy (and grayscale) array.
img = np.array(img.convert('L').getdata()).reshape(img.size)

I will do the same with cartopy just to be sure that srtm.py is getting the right image.

In [4]:

import cartopy.crs as ccrs
from cartopy.io.srtm import srtm_composite, add_shading, fill_gaps

elev, crs, extent = srtm_composite(lon - dx/2., lat - dy/2., dx, dy)
elev = ma.masked_equal(elev, -32768)

In [5]:

import matplotlib.pyplot as plt

subplot_kw = dict(projection=ccrs.PlateCarree())
fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(13, 13),
                               subplot_kw=subplot_kw)

kw = dict(extent=extent, transform=crs, cmap='Greys')
ax0.imshow(img, origin='upper', **kw)
ax0.plot(lon, lat, 'ro', alpha=0.5)
ax0.set_title('srtm.py')
gl = ax0.gridlines(draw_labels=True)
gl.xlabels_top = gl.ylabels_right = False

ax1.imshow(elev, origin='lower', **kw)
ax1.plot(lon, lat, 'ro', alpha=0.5)
ax1.set_title('cartopy')
gl = ax1.gridlines(draw_labels=True)
gl.xlabels_top = gl.ylabels_left = False

Weird, the same features are there, but I must be plotting it wrong or the data source is different... I must investigate this further another time. For now let's trust what srtm.py is doing and add the elevation data to the GPX file.

In [6]:

def get_line(gpx):
    previous_point = None
    length = 0
    lengths, elevations = [], []
    for point in gpx.walk(only_points=True):
        if previous_point:
            length += previous_point.distance_2d(point)
        previous_point = point
        lengths.append(length)
        elevations.append(point.elevation)
    return lengths, elevations

# Raw GPS.
gpx = gpxpy.parse(open(fname))
lengths_raw, elevs_raw = get_line(gpx)

# Raw SRTM.
gpx = gpxpy.parse(open(fname))
data = srtm.get_data()
data.add_elevations(gpx)
lengths_srtm, elevs_srtm = get_line(gpx)

# Smoothed SRTM.
gpx = gpxpy.parse(open(fname))
data = srtm.get_data()
data.add_elevations(gpx, smooth=True)
lengths_smoo, elevs_smoo = get_line(gpx)

In [7]:

import seaborn

fig, ax = plt.subplots(figsize=(7, 3.5))
ax.plot(lengths_raw, elevs_raw, label='Raw GPS')
ax.plot(lengths_srtm, elevs_srtm, label='SRTM')
ax.plot(lengths_smoo, elevs_smoo, label='Smoothed SRTM')
l = ax.legend()

The raw GPS altitude was not that bad (just one spike that could be easily removed). The SRTM data needs to be smoothed to avoid the step like plot, specially in small areas like this one. But at least they all agree, anyone that hiked Morro do Pai Inácio knows that this profile seems about right.

Looking forward to test this on my next hiking trip! I will end this post with this sunset photo taken that day (by a very special someone).

In [8]:

from IPython.display import Image
Image('../../images/sunset.jpg')

Out[8]:

In [9]:

HTML(html)

Out[9]:

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

Using SRTM data to correct GPS altimetry

Comments