Monday, March 5, 2018

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.     
Elevation Profile
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).
Elevation Profile of Mountain  Table, Cape Town
Figure 2. Elevation Profile of Mountain  Table, Cape Town
Table Montain, Cape Town. South Africa (Captured from Google Earth)
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.

8 comments :

  1. This analysis uses spherical trigonometry. A more accurate method would use an ellipsoid, such as WGS84, but it wouldn't matter much at this scale.

    ReplyDelete
  2. Hi I'm triying to use your code, but I get an error when run the code, as follow:

    Traceback (most recent call last):
    File "C:/Users/HORUS10/PycharmProjects/PERFIL DE ELEVACION/PERFIL_ELEVACION.py", line 58, in
    location['locations']=d_ar
    NameError: name 'location' is not defined

    Could you help me.

    Thanks, best regards.

    Javier

    ReplyDelete
    Replies
    1. Hi Javier,

      Thanks for spotting the bug. Change the line to:location={"locations":d_ar}. It should work.

      Delete
  3. I simply wanted to write down a quick word to say thanks to you for those wonderful tips and hints you are showing on this site.
    It’s great to come across a blog every once in a while that isn’t the same out of date rehashed material. Fantastic read.

    Python Training in Chennai | Python Training Institutes in Chennai

    ReplyDelete
  4. Sir how can we load las file in 3.0

    ReplyDelete
  5. Needed to compose you a very little word to thank you yet again regarding the nice suggestions you’ve contributed here.
    Python Training in Bangalore

    ReplyDelete
  6. Oh, cool, you used my Open Elevation API! It's always nice to see it in use! Mind if I ask you a few questions?

    While you were using it, did you find any particular annoyances or features you'd like to see cleared up or implemented?

    Do you think you'll use it for other projects?

    How did you find it?

    If a payed option was provided, with faster response times, higher precision and a dashboard with your statistical historical usage (e.g. regions of the map you query the most) was available, would you consider it? If so, for what price?

    I am asking these questions because the service has become quite overloaded in recent months, and we're looking into ways of improving it in a self-sustaining way, which could improve the overall public offering (hardware upgrades, increase in dataset resolution, etc). Thanks in advance!

    ReplyDelete
    Replies
    1. Hi, Thanks for your Open Elevation API. It's great. I found it from online searching. I wanted to create an elevation profile with python and searching for some open elevation sources, and I find yours. So far, I just use it for the application in this post. It seems you're using SRTM 1 arc second resolution. So if you could use higher precision and deploy some useful tool or application using the data, I think it's worth for some pricing. How much? It depends on the service you're offering compare to others.

      Delete