# How to Create Contour Lines in QGIS

This tutorial will discuss how to create contour lines in QGIS. But before we are discussing about drawing contour lines in QGIS, let's start talking a little bit about the concept of contour line, what contour line definition is, contour lines characteristics and and how the method in drawing contour line . In the last section we will see how to create contour lines from points and also from DEM in QGIS. For this tutorial I am using QGIS 3.0.

### What is Contour Line?

Contour is a height isoline. Then what is an isoline? Isoline is a line on a map with a constant value. Therefore the contour definition is a line which has a constant elevation along the line. There are some characteristics of contour line such as contour interval, contour elevation and distance between contour lines. If we see a contour map we found many lines on it. Those lines are never crossing to each other and each line has a unique elevation value that define its elevation on the map. Contour interval defines height interval between the line. For example if a contour map has 5 m interval contour, it means the height different between the adjacent contour line is 5 m. If we know the distance between contour lines then we can calculate the slope. The rule of thumb in defining contour interval is 1/2000 x map scale.

### How to Draw Contour Lines

The basic method in drawing contour line is interpolation from a point dataset. The simplest interpolation method is linear interpolation which is frequently used in drawing contour lines manually.  In drawing a contour line with linear interpolation, at least three steps need to be done:
1. Connecting each point to form a TIN (Triangular Irregular Network) structure.
2. Do interpolation for each line in the TIN.
3. Connecting the interpolation point with the same elevation.
Figure 1 shows step by step approach in drawing contour line from connecting each point to form TIN, do interpolation along the TIN line and drawing contour line by connecting the same elevation of interpolation points.
 Figure 1. Contour line drawing
We already got some basic concepts about contour line. Now let's move to technical part how to create contour lines with QGIS. Firstly we will see how to create contour line from a set of height points.

### Create Contours from Points

To create contour lines from points, we are using Contour plugin. As usual we can get the plugin from Plugins menu >> Manage and install plugins. In the search option type "Contour". You will see the Contour plugin appears as in figure 2. Check it and Install the plugin if you don't have it in your QGIS 3.

 Figure 2. Installing Contour Plugin
Now, let's make contour lines from points dataset with the following steps:

1. Add the points into QGIS map canvas as in figure3.
 Figure 3. Points dataset
2. Open the Contour plugin by clicking the icon. The contour plugin window will appear as in figure 4.
3. Select the point's layer name and the elevation data field. In the contouring menu, we can choose  to create contour line, filled contour or both. Next, we have to define the number of contour lines. We can set any number, but I suggest the number of contour line follows the rule of thumb as mentioned above. For example if we want to make a contour map in 100.000 scale, then the contour interval will be 1/2000 x 100.000 = 50 m. Then the number of contour line will be: maximum elevation/contour interval. 3554/50=71 lines. When we set the calculated number, the contour interval can be seen in the right list. We can tweak the number if the interval is not correct. That's why I change the number to 72. Moreover the interval can be edited manually by double click the interval number.
4. Next we can specify the output name. The output will be written to memory, so we need to export to it to a GIS file to keep it permanently.
5. In the next setting we can specify the precision and unit for the contour label. Furthermore the color graduation for the contour lines can be changed by selecting a color ramp for the output contour lines.
6. When finish click the Add button. The result can be seen in the QGIS map canvas without closing the contour plugin window, therefore we can switch back to change a parameter if the result not looks as we want. Figure 5 shows the result of contour lines.

 Figure 4. Contour plugin window

 Figure 5. Contour lines result from points

### Create Contours from DEM

In this part we will create contour lines from DEM. To create contour lines from DEM we can use the Contour tool in the Raster menu >> Extraction >> Contour. See figure 6.

 Figure 6. Contour tool
To create the contour lines from DEM, can be done with several steps:
1. Add DEM data into QGIS map canvas, as in figure 7. If you don't have one. You can download DEM data directly from QGIS.

 Figure 7. DEM data in QGIS map canvas
2. Open the Contour tool. The contour tool window will appear as in figure 8.

 Figure 8. Raster contour tool
3. Select the DEM data. Then set the contour interval, in this case  I used 50 m. In the Contours option, we can specify the contour output path or save it to temporary file.
4. If all the settings already set. Click the Run in Background button. The result will shown in QGIS  map canvas as in figure 9.

 Figure 9. Contour lines from DEM

### Contour Smoothing

If we compare the result, the contour lines from DEM looks smoother than points. That's because the DEM contains more elevation information compare to points dataset. The DEM has information in each pixel, so it could consist of thousands pixel and elevation information. But don't worry. We can make the contour line from points smoother using Smooth tool as follow:

1.  In the processing toolbox type "smooth". The smooth tool will appear under  Vector Geometry, as in figure 10.
 Figure 10. Vector geometry smooth tool
2. Open the smooth tool. Select the input contour layer. Set the number of iteration, offset and maximum angle to smooth. The explanation of each parameter can be seen in the help window on the right. See figure 11.

 Figure 11. Contour smoothing
3. If all the parameters have been set. Click Run in Background button to start the smoothing process. Figure 12 shows the comparison of contour lines before and after smoothing.
 Figure 12. Contour lines before and after smoothing
That's all the tutorial how to create contour lines in QGIS both from points and DEM. I hope it's useful. Please give comment if you find another tool in QGIS or any improvement in creating contour lines using QGIS software.

# QGIS Tutorial: How to Create an Action

The previous tutorial discussed about how to perform spatial query in QGIS 3. In this tutorial, I will discuss about Actions in QGIS. What is Action? Action is a command or task to be executed based on a click trigger. For example if a user click a feature then an image will be opened. Here the click is the trigger and opening an image is the action.

We can define some types of actions in QGIS like opening an application or file in an operating system (Windows, Mac or UNIX), open a web page,  or create a specific action with a python script. In this tutorial I will cover how to create an open action and a python action. The action example is to open an URL or web link of traffic camera in Vancouver city. For that, we will use traffic camera location data that can be downloaded from Vancouver City Data Catalogue. For this tutorial I'm using QGIS 3.0.

### Opening an URL Action

Let's start with creating a simple action to view a camera traffic when a camera point is clicked on QGIS map canvas. The action can be done with the following steps:

1. Add a basemap. For this case I'm using Open Street Map.
2. Add traffic camera data into QGIS map canvas. The camera location over the basemap will look like as in figure 1.

 Figure 1. Vancouver City Traffic Camera Location
3. Open camera location layer properties and choose Actions icon. In this step, by clicking Create default actions, some actions will be listed in the Action List. To create a new one click the green + button as in figure 2.

 Figure 2. Add a new action
4. An Add action window will appear as in figure 3. In the Type option select Open. Give a description and a short name, for example I gave the name open traffic camera. In the Action Scope check Canvas and Feature Scope, cause we want the action will be executed for the feature in QGIS map canvas.
5. Next in the Action text, type [%URL%]. URL is a column name in the attribute table which contains a HTML link to each  traffic camera. You can see how the attribute looks like in figure 4.

 Figure 4. Traffic camera attribute table
6. Last step after clicking OK button both in the Add new action window and layer properties, let's run the action by selecting the action from Action Icons in the toolbar as in figure 5.
 Figure 5. Select an Action
The mouse pointer will change into cross hair. Click a camera point. Soon the URL of the camera will be opened if the default web browser as in figure 6.
 Figure 6. Traffic camera view in web browser

### Python URL Action

In this part we will make an action with Python script for the same example. The action will open a QT web browser instead of the default web browser. Most steps are similar like above, except:
1. Select Python in the type option.
2. Give a new name for example: open traffic camera python.
3. In the Action text, type the listing code below. Figure 7 shows the new action with python.

```from PyQt5 import QtCore
from PyQt5.QtCore import *
from PyQt5.QtWebKit import *
from PyQt5.QtWebKitWidgets import *
web = QWebView()
web.show()
```

 Figure 7. QGIS Action with Python
When the action is executed by clicking a traffic camera point. A QT web browser will be opened as in figure 8.
 Figure 8. Traffic camera view in QT web browser
Congratulation! Now you can observe the traffic condition around Vancouver City.
That's all the tutorial how to create an Action in QGIS. I think this is a great tool, cause gives more functionality to our map application. If you want to explore more, visit QGIS actions documentation.

# Creating GPX File Viewer in Python

GPX file is a GPS data that stored in XML format. There are some types of data that stored in GPS namely waypoint, route and track. Track is a type of data which is recorded regularly by  GPS in an interval. Therefore, with GPS tracker data we can visualize a trip that we did, which road we passed by, the length of track, time taken, etc. There are some tools or apps that available to view or visualize a GPS track file such as ArcGIS, QGIS, Google Earth and also some online GPX viewer app also available to visualize a GPS tracker data. Just do online searching and you will find many of them.

In this tutorial, I will discuss how to create a GPX track file viewer in Python. I'm using Jupyter Notebook with Python 3. So if you want to run the code make sure you have Jupyter Notebook installed in your machine. We will visualize the track itself, speed and elevation profile along the track. Moreover the visualization will be dynamic in an animation that plot each GPS tracking point as in the figure below.
To create the GPX file viewer we need some modules like matplotlib, time, IPython, xml dom and math. Import all the required modules as in the code below.

```import matplotlib.pyplot as plt
import time
from IPython import display
from xml.dom import minidom
import math  ```

Then open a GPX file and parse the data using minidom. Here we are taking some important elements: 'trkpt', 'ele' and 'time' for track point, elevation and time.

```#READ GPX FILE
data=open('F:/ride.gpx')
xmldoc = minidom.parse(data)
track = xmldoc.getElementsByTagName('trkpt')
elevation=xmldoc.getElementsByTagName('ele')
datetime=xmldoc.getElementsByTagName('time')
n_track=len(track)
```

After getting those elements, we need to get the value from each elements. Then we need to parse the elements to get latitude, longitude, elevation and time. At the end, the time is converted to second. It will be used later in speed calculation.

```#PARSING GPX ELEMENT
lon_list=[]
lat_list=[]
h_list=[]
time_list=[]
for s in range(n_track):
lon,lat=track[s].attributes['lon'].value,track[s].attributes['lat'].value
elev=elevation[s].firstChild.nodeValue
lon_list.append(float(lon))
lat_list.append(float(lat))
h_list.append(float(elev))
# PARSING TIME ELEMENT
dt=datetime[s].firstChild.nodeValue
time_split=dt.split('T')
hms_split=time_split[1].split(':')
time_hour=int(hms_split[0])
time_minute=int(hms_split[1])
time_second=int(hms_split[2].split('Z')[0])
total_second=time_hour*3600+time_minute*60+time_second
time_list.append(total_second)
```

Now we have all parameters that are needed to create the GPX file viewer. To do further calculation like distance and speed, we create some functions namely geo2cart, distance and speed. geo2cart is a function to convert Geodetic coordinate into Cartesian coordinate with WGS-84 ellipsoid reference datum. We need the Cartesian coordinate because it will be used to calculate a distance in meter between two GPS tracking points. After getting the distance, then the speed in m/s can be calculated by divided with time difference.

```#GEODETIC TO CARTERSIAN FUNCTION
def geo2cart(lon,lat,h):
a=6378137 #WGS 84 Major axis
b=6356752.3142 #WGS 84 Minor axis
e2=1-(b**2/a**2)
return X,Y

#DISTANCE FUNCTION
def distance(x1,y1,x2,y2):
d=math.sqrt((x1-x2)**2+(y1-y2)**2)
return d

#SPEED FUNCTION
def speed(x0,y0,x1,y1,t0,t1):
d=math.sqrt((x0-x1)**2+(y0-y1)**2)
delta_t=t1-t0
s=float(d/delta_t)
return s
```

Next step, using the function we do calculation of distance and speed. The result are stored in distance and speed list.

```#POPULATE DISTANCE AND SPEED LIST
d_list=[0.0]
speed_list=[0.0]
l=0
for k in range(n_track-1):
if k<(n_track-1):
l=k+1
else:
l=k
XY0=geo2cart(lon_list[k],lat_list[k],h_list[k])
XY1=geo2cart(lon_list[l],lat_list[l],h_list[l])

#DISTANCE
d=distance(XY0[0],XY0[1],XY1[0],XY1[1])
sum_d=d+d_list[-1]
d_list.append(sum_d)

#SPEED
s=speed(XY0[0],XY0[1],XY1[0],XY1[1],time_list[k],time_list[l])
speed_list.append(s)
```

Now we can visualize the GPX file track data with three plots: Track, speed and elevation profile.

```#PLOT TRACK
f,(track,speed,elevation)=plt.subplots(3,1)
f.set_figheight(8)
f.set_figwidth(5)
track.plot(lon_list,lat_list,'k')
track.set_ylabel("Latitude")
track.set_xlabel("Longitude")
track.set_title("Track Plot")

#PLOT SPEED
speed.bar(d_list,speed_list,30,color='w',edgecolor='w')
speed.set_title("Speed")
speed.set_xlabel("Distance(m)")
speed.set_ylabel("Speed(m/s)")

#PLOT ELEVATION PROFILE
base_reg=0
elevation.plot(d_list,h_list)
elevation.fill_between(d_list,h_list,base_reg,alpha=0.1)
elevation.set_title("Elevation Profile")
elevation.set_xlabel("Distance(m)")
elevation.set_ylabel("GPS Elevation(m)")
elevation.grid()
```

Lastly in a loop, dynamic plot is created in one second interval by plotting each tracking points over the track and elevation profile. The speed change is visualized in a bar chart.

```#ANIMATION/DYNAMIC PLOT
for i in range(n_track):
track.plot(lon_list[i],lat_list[i],'yo')
speed.bar(d_list[i],speed_list[i],30,color='g',edgecolor='g')
elevation.plot(d_list[i],h_list[i],'ro')
display.display(plt.gcf())
display.clear_output(wait=True)
time.sleep(1)

plt.show()
```

That's all the tutorial how to create a GPX file viewer to visualize track, speed and elevation profile with Python in Jupyter Notebook. You can try to run the code by downloading the GPX file viewer Jupyter Notebook and a GPX file sample that used in this tutorial. Finally use your own GPX track file and try to visualize it with the code.

# Spatial Query in QGIS 3.0

Spatial Query is selection of features that satisfies a certain condition which relates to other features in a space. Some example of spatial queries such as features that intersect with other features (i.e. parcel that intersect with road network), features within other features (i.e. hotspot within a forest area), features with a distance from another feature (i.e. buildings 20 meter from main road), etc. In QGIS 2.x we can perform spatial query with spatial query plugin, but for QGIS 3, I can't find the plugin in the plugin repository. So how to do spatial queries in QGIS 3.0?

In QGIS 3.0, the spatial queries operation can be performed with a tool called Select by Location. The tool can be found under the Vector Selection tool in the Processing Toolbox (See figure 1).Let see some spatial queries example using this tool.
 Figure 1. Select by Location and Select by Expression tool

For spatial query example, I used some Vancouver dataset such as: crime point 2016, parks polygon and public streets. Those data can be downloaded from Vancouver City data catalogue. Those data can be seen as in figure 2.
 Figure 2. Crimes, city park and street dataset of Vancouver
Based on the data we will do some spatial queries below.
1. Select crime events that happened in the city park
2. Which park that crimes took place?
3. Find crime event that took place on the street.
4. Find roads which cross city park
5. Find crime that occurred at radius 100 meter from main street.

For first query we will select all crimes that took place in the city park. We can do the query with the following steps:
1. Select crime_shp as the selected features.
2. Check are within option below geometric predicate.
3. Select park_polygon as comparing features.
4. Select creating new selection, because we make a new selection (not to select from selected features).
5. Run the the query with the selected option as in figure 3. All crimes that happened in the city park will be selected. I captured a part of the query result in figure 4.

 Figure 3. Query option to select all crimes event in the city park
 Figure 4. Spatial query result: Crimes within park
For the second question, we use contain. The logic is we select any park that contains crime point. The query is set as in figure 5. Figure 6 shows the result.
 Figure 5. Spatial query to find any park where crime took place

 Figure 6. The selected park where crimes took place
The third query is to find any crime that happened on the street. For this we can use touch or overlap. But I think touch is more correct than overlap, because overlap is more suits with polygon features.   While street is a line/polyline features, so in this query we just need to find a point that touch the line of a street. If there is any crime point that touch a street then it will be selected. But I got none.

In the fourth query we want to select any road in a park. For this case we use intersect option. There is also a cross option but it will not select a street segment that completely within a park. The figure 7 gives the difference between cross and intersect.

 Figure 7. Cross vs Intersect
The last spatial query is a little bit complex, we want to find any crime that happened within a radius of 100 m from main/arterial road. In this query there is a dynamic variable, radius length. For this case we can't use Select by location tool right away because there is no option to set define a radius. Because of that, firstly we have to make a buffer for arterial road with a distance of 100 m, and then we select any crime points within the buffer area. The query can be done with the following steps:
1. Make a definition query to select arterial road only.
2. Buffer 100 m the arterial road.
3. Use Select by location tool to find any crime within the buffer area.

A definition query can be made through street's layer properties in the source option. In the option select Query Builder and make a query expression to select feature by its attribute. For this case the expression is USE='Arterial'. See figure 8.
 Figure 8. Make a definition query
Then using Buffer tool, buffer the arterial road 100 m. The result can be seen as in figure 9.

Finally using Select by location tool select any crime point within the buffer area. Figure 10 shows the result crime that took place 100 m from main/arterial road.
 Figure 10. Crime event 100 m from main road
That's all the tutorial how to perform spatial query in QGIS 3.0. With the example above hopefully you can get an idea how to do spatial query for you own case. If you want to explore more, QGIS 3.0 documentation provides more explanation about vector query.

# This is How to Add Google Maps Layers in QGIS 3

Google Maps Layer is hosting somewhere in Google server and sends the tile to the user who request it. Technically we call this as Tile Map Service (TMS). So, we just need to find the TMS which Google uses to serve Google Maps layers. After searching for a while, finally I found the Google Maps Layer TMS from NextGIS that are listed below:

Then how to add the Google Maps tile service in QGIS 3? Simply, we uses XYZ Tiles. For example we want to add a Google Maps layer as a basemap. We can do it with the following steps:
1. Make a new connection, with right click the XYZ Tiles as in figure 1.

 Figure 1. Make a new connection in XYZ Tiles
2. Give the name, for example: Google Maps.  In the URL copy and paste the Google Maps TMS above. Set the max zoom level to 19. See figure 2.

 Figure 2. XYZ Tiles connection details
3. Add the Google Maps XYZ tiles connection that already created to the QGIS layer. The Google Maps layers basemap will be added to the QGIS map canvas direclty as in figure 3.

 Figure 3. Google Maps Layers in QGIS

With the same steps, the other Google basemap layers like Google Terrain, Google Satellite, etc, can be added to QGIS 3. I captured some of the results in the images below.

 Figure 4. Google Terrain Layer in QGIS 3

 Figure 5. Google Satellite Layer in QGIS 3
Great, now we can add a basemap from Google Maps layer in QGIS 3, even without any plugin. Hope this tutorial helps :)

# 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):
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"
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_str=res_byte.decode("utf8")
#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):
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"
fp=urllib.request.urlopen(response)

#RESPONSE PROCESSING
res_str=res_byte.decode("utf8")
#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.