Creating An Elevation Profile Generator in Python

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.     
Figure 1. Elevation Profile

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).
Figure 2. Elevation Profile of Mountain Table, Cape Town
Figure 3. Table Montain, Cape Town. South Africa (Captured from Google Earth)
Find the complete code of elevation profile generator from online Open Elevation API using python 3 below.

"""
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.

Related Posts
Disqus Comments