Elevation profile is a two dimensional cross section along a line or path. It is very helpful to visualize the elevation change along a line or path. In this post I will discuss how to create an elevation profile graph between two points using Open Elevation API in Python 3 as in figure 1.
Open Elevation is a free and open-source elevation API. It uses DEM data from SRTM, that is published and can be downloaded freely from USGS or using QGIS (Download DEM SRTM with QGIS). Current SRTM dataset has resolution 1 arc-second (around 30 m) with near global coverage from 56°S to 60°N.
Open Elevation is a free and open-source elevation API. It uses DEM data from SRTM, that is published and can be downloaded freely from USGS or using QGIS (Download DEM SRTM with QGIS). Current SRTM dataset has resolution 1 arc-second (around 30 m) with near global coverage from 56°S to 60°N.
Importing Libraries
We are starting to build the elevation profile application with importing some libraries such as urllib, json, math and matplotlib.
import urllib.request import json import math import matplotlib.pyplot as plt
Urllib library will be used to send a request to Open Elevation server. The request is made by sending a json data which contains point(s) location(latitude,longitude). The math library will be used to do some processing or calculation such as distance, average elevation, minimum and maximum elevation. Lastly the elevation profile will be plotted using matplotlib library.
Generating Elevation Points
To create an elevation profile a long a line we need two points as starting and end point. Those points are defined as P1 and P2. Then we need to generate a set of points between those two points. In this case I created a number of 100 points. You can change the number of points, but must be kept in mind, too fewer points will not give a good representation of a real terrain's profile on the other hand too many points will make the Open Elevation API server to response longer.#START-END POINT P1=[latitude,longitude] P2=[latitude,longitude] #NUMBER OF POINTS s=100 interval_lat=(P2[0]-P1[0])/s #interval for latitude interval_lon=(P2[1]-P1[1])/s #interval for longitude #SET A NEW VARIABLE FOR START POINT lat0=P1[0] lon0=P1[1] #LATITUDE AND LONGITUDE LIST lat_list=[lat0] lon_list=[lon0] #GENERATING POINTS for i in range(s): lat_step=lat0+interval_lat lon_step=lon0+interval_lon lon0=lon_step lat0=lat_step lat_list.append(lat_step) lon_list.append(lon_step)
Calculate Distance Between Two Elevation Points
The distance between each elevation points are calculated using Haversine Formula. A function called haversine is created with four parameters, there are: latitude and longitude of start-end point. The distances then are stored in a list which called d_list.
#HAVERSINE FUNCTION def haversine(lat1,lon1,lat2,lon2): lat1_rad=math.radians(lat1) lat2_rad=math.radians(lat2) lon1_rad=math.radians(lon1) lon2_rad=math.radians(lon2) delta_lat=lat2_rad-lat1_rad delta_lon=lon2_rad-lon1_rad a=math.sqrt((math.sin(delta_lat/2))**2+math.cos(lat1_rad)*math.cos(lat2_rad)*(math.sin(delta_lon/2))**2) d=2*6371000*math.asin(a) return d #DISTANCE CALCULATION d_list=[] for j in range(len(lat_list)): lat_p=lat_list[j] lon_p=lon_list[j] dp=haversine(lat0,lon0,lat_p,lon_p)/1000 #km d_list.append(dp) d_list_rev=d_list[::-1] #reverse list
Constructing JSON and Sending Request
The Open Elevation API receives request in JSON format. Therefore the points must be set into JSON as in specified form. Then the data in JSON format is sending to the server.#CONSTRUCT JSON d_ar=[{}]*len(lat_list) for i in range(len(lat_list)): d_ar[i]={"latitude":lat_list[i],"longitude":lon_list[i]} location={"locations":d_ar} json_data=json.dumps(location,skipkeys=int).encode('utf8') #SEND REQUEST url="https://api.open-elevation.com/api/v1/lookup" response = urllib.request.Request(url,json_data,headers={'Content-Type': 'application/json'}) fp=urllib.request.urlopen(response)
Processing The Elevation Response
If the Open Elevation server successfully response the request, then we have to process the result before using it. The processing steps include read and decoding the result. Then we extract the elevation of each point from the result to be plotted as elevation profile in the next step. Some basic stats are also computed such as average elevation, maximum and minimum elevation.
#RESPONSE PROCESSING res_byte=fp.read() res_str=res_byte.decode("utf8") js_str=json.loads(res_str) #print (js_mystr) fp.close() #GETTING ELEVATION response_len=len(js_str['results']) elev_list=[] for j in range(response_len): elev_list.append(js_str['results'][j]['elevation']) #BASIC STAT INFORMATION mean_elev=round((sum(elev_list)/len(elev_list)),3) min_elev=min(elev_list) max_elev=max(elev_list) distance=d_list_rev[-1]
Plotting Elevation Profile
Finally we plot elevation profile with distance as x-axis and elevation as y-axis. Including average elevation, minimum and maximum elevation.#PLOT ELEVATION PROFILE base_reg=0 plt.figure(figsize=(10,4)) plt.plot(d_list_rev,elev_list) plt.plot([0,distance],[min_elev,min_elev],'--g',label='min: '+str(min_elev)+' m') plt.plot([0,distance],[max_elev,max_elev],'--r',label='max: '+str(max_elev)+' m') plt.plot([0,distance],[mean_elev,mean_elev],'--y',label='ave: '+str(mean_elev)+' m') plt.fill_between(d_list_rev,elev_list,base_reg,alpha=0.1) plt.text(d_list_rev[0],elev_list[0],"P1") plt.text(d_list_rev[-1],elev_list[-1],"P2") plt.xlabel("Distance(km)") plt.ylabel("Elevation(m)") plt.grid() plt.legend(fontsize='small') plt.show()
Figure 2 is an elevation profile along a line across the Mountain Table in Cape Town, South Africa (figure 3), with P1(-33.965526,18.377388) and P2(-33.959881,18.478549).
"""
ELEVATION PROFILE APP GENERATOR
ideagora geomatics-2018
http://geodose.com
"""
import urllib.request import json import math import matplotlib.pyplot as plt #START-END POINT P1=[latitude,longitude] P2=[latitude,longitude] #NUMBER OF POINTS s=100 interval_lat=(P2[0]-P1[0])/s #interval for latitude interval_lon=(P2[1]-P1[1])/s #interval for longitude #SET A NEW VARIABLE FOR START POINT lat0=P1[0] lon0=P1[1] #LATITUDE AND LONGITUDE LIST lat_list=[lat0] lon_list=[lon0] #GENERATING POINTS for i in range(s): lat_step=lat0+interval_lat lon_step=lon0+interval_lon lon0=lon_step lat0=lat_step lat_list.append(lat_step) lon_list.append(lon_step) #HAVERSINE FUNCTION def haversine(lat1,lon1,lat2,lon2): lat1_rad=math.radians(lat1) lat2_rad=math.radians(lat2) lon1_rad=math.radians(lon1) lon2_rad=math.radians(lon2) delta_lat=lat2_rad-lat1_rad delta_lon=lon2_rad-lon1_rad a=math.sqrt((math.sin(delta_lat/2))**2+math.cos(lat1_rad)*math.cos(lat2_rad)*(math.sin(delta_lon/2))**2) d=2*6371000*math.asin(a) return d #DISTANCE CALCULATION d_list=[] for j in range(len(lat_list)): lat_p=lat_list[j] lon_p=lon_list[j] dp=haversine(lat0,lon0,lat_p,lon_p)/1000 #km d_list.append(dp) d_list_rev=d_list[::-1] #reverse list #CONSTRUCT JSON d_ar=[{}]*len(lat_list) for i in range(len(lat_list)): d_ar[i]={"latitude":lat_list[i],"longitude":lon_list[i]} location={"locations":d_ar} json_data=json.dumps(location,skipkeys=int).encode('utf8') #SEND REQUEST url="https://api.open-elevation.com/api/v1/lookup" response = urllib.request.Request(url,json_data,headers={'Content-Type': 'application/json'}) fp=urllib.request.urlopen(response) #RESPONSE PROCESSING res_byte=fp.read() res_str=res_byte.decode("utf8") js_str=json.loads(res_str) #print (js_mystr) fp.close() #GETTING ELEVATION response_len=len(js_str['results']) elev_list=[] for j in range(response_len): elev_list.append(js_str['results'][j]['elevation']) #BASIC STAT INFORMATION mean_elev=round((sum(elev_list)/len(elev_list)),3) min_elev=min(elev_list) max_elev=max(elev_list) distance=d_list_rev[-1] #PLOT ELEVATION PROFILE base_reg=0 plt.figure(figsize=(10,4)) plt.plot(d_list_rev,elev_list) plt.plot([0,distance],[min_elev,min_elev],'--g',label='min: '+str(min_elev)+' m') plt.plot([0,distance],[max_elev,max_elev],'--r',label='max: '+str(max_elev)+' m') plt.plot([0,distance],[mean_elev,mean_elev],'--y',label='ave: '+str(mean_elev)+' m') plt.fill_between(d_list_rev,elev_list,base_reg,alpha=0.1) plt.text(d_list_rev[0],elev_list[0],"P1") plt.text(d_list_rev[-1],elev_list[-1],"P2") plt.xlabel("Distance(km)") plt.ylabel("Elevation(m)") plt.grid() plt.legend(fontsize='small') plt.show()
I am sure, there are still many rooms to improve or to optimize the code, so we can make an elevation profile with more information such as the change slope or an elevation profile along a path. Trying a smoothing algorithm also a nice idea, so we can get a smooth elevation profile. For that, just use and modify the Python code and tell me how it goes.
