Summary & code
Use the swisstopo API to get a profile of altitudes with the function get_swisstopo_elevation_profile by passing coordinates aka “Rechtswert” and “Hochwert” as array/list (any iterable should do). The gist on GitHub should always be the most up-to-date version
def get_swisstopo_elevation_profile(coords, # a path (2+ points in the form)
kind='csv',
# three heights are available for JSON
# COMB, DTM2, DTM25
which='COMB',
asnumpy=True,
quiet=True,
opts={"sr": None,
"nb_points": None,
"offset": None,
"distinct_points": True,
"callback": None,
}
):
"""
Call the swisstopo API for altitude data along a path.
Pass in a list or array of coordinates in EPSG 2056 (LV95) or 21781 (LV03)
to get altitude values along the provided path. For options see keywords or
parameters of call via swisstopo API.
Parameters
----------
coords : list or numpy array
The coordinates in either EPSG 2056 (LV95) or EPSG 21781 (LV03).
kind : str, optional
Which API backend should be queried. Available are json and csv.
If none of these are passed properly, fallback is csv.
The default is 'csv' (easier structure to parse).
which: str
If kind is json, three altitude values are available, DTM2, DTM25 and
COMB(INATION).
asnumpy : bool, optional
Whether to return a numpy array.
The default is True.
quiet : bool, optional
Whether to quiet the output (True) or not.
The default is True.
opts: dict
Further keywords than can be passed to the API call, see
https://api3.geo.admin.ch/services/sdiservices.html#profile
Returns
-------
list or numpy array
The returned points, in the form of distance along path, altitude,
coordinates.
"""
import requests
import numpy as np
if len(coords) > 5000:
print('Warning, number of coordinates exceeds swisstopo API'
'max number of points, reduce coords beforehand.')
return np.asarray([]) if asnumpy else []
payload = 'geom={"type":"LineString","coordinates":['
payload += ','.join([f"[{coord[0]},{coord[1]}]"
for coord in coords])
payload += ']'
opts = ','.join(['"'+str(key)+'"'+':'+'"'+str(opt)+'"'
for key, opt in opts.items()
if opt is not None])
if opts:
payload += ',' + opts
payload += '}'
kind = kind.lower()
if kind.lower() not in ['csv', 'json']:
if not quiet:
print('Only csv or json can be chosen, autoselecting csv')
kind = 'csv'
baseurl = "https://api3.geo.admin.ch/rest/services/profile." + kind
try:
profile = requests.get(baseurl, params=payload)
except ConnectionError:
if not quiet:
print('Connection timeout')
return np.asarray([]) if asnumpy else []
if profile.ok:
if not quiet:
print('Success')
else:
if not quiet:
print('Failed')
if kind == 'csv':
profile = profile.text.split('\r\n')
# Distance, Altitude, Easting, Northing -> not needed
# header = profile[0]
profile = list(map(lambda x: [float(j.strip('"'))
for j in x.split(';')],
profile[1:-1]))
elif kind == 'json':
profile = [[p['dist'], p['alts'][which],
p['easting'], p['northing']]
for p in profile.json()]
if asnumpy:
profile = np.asarray(profile)
return profile
if __name__ == '__main__':
# straightforward example usage
rw=[2611025.0, 2620975.0, 2633725.0]
hw=[1266400.0, 1256750.0, 1250000.0]
profile = get_swisstopo_elevation_profile(list(zip(rw, hw)))
import matplotlib.pyplot as plt
plt.fill_between(profile[:,0],
profile[:,1],
facecolor='grey',
edgecolor='k',
alpha=0.8)
plt.ylim(min(profile[:,1]), None)Background/motivation
For the CLOUDLAB project we regularly get forecasts that contain a profile at the bottom of the figure from MeteoSwiss that illustrates the terrain nicely (albeit, sometimes it is missing for unknown reasons):
Similarly, I produced a lot of figures of our scans and depending on azimuth and elevation they may be affected by ground clutter, e.g. at a distance of 2 km the increased reflectivity is caused by the terrain below (the hill at around 800 m above sea level). As you can see, I added the terrain already.
Initially, I thought I’d have to get the actual DEM for the area, find the path matching the scan direction and calculate the profile myself. While this might actually be a bit more accurate with a good DEM, it would be more work, not be transferable and mean I’d have to have extra data around. Instead, I realized that swisstopo offers the measure tool on map.geo.admin.ch and uses their own API with a matching request. So I choose to use this as I’m already adding the coordinates for the scan and it is relatively straightforward (min/max of lowest elevation scanline) to get the path we need. The main issue I faced was a misformatted quote that came from copy-pasting the swisstopo example. After finally figuring this out after a detailed URL decode and string comparison the function is relatively lean and can use either the JSON or the CSV backend of the API and by default gives you out a numpy array which can be visualized with ease and can look as follows:
I hope this helps someone to get a profile to enrich their Python (or other) graphs when measuring in Switzerland. In the likely case that you aren’t doing something in Switzerland it might be worthwhile to check out the Google Maps Elevation API (for which you need an API key and billing enabled)
One Reply to “geo utils – get altitude profile data from swisstopo”
Comments are closed.