Saturday, April 14, 2018

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.

QGIS Action Vancouver City Traffic Camera Location
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.

Add a new action in QGIS
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.

Edit an action in QGIS
Figure 3. Add/Edit an action
Traffic camera attribute table QGIS Action
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.
Select an Action QGIS
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.
Traffic camera view in web browser
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.load(QUrl("[%URL%]"))
web.show()

QGIS Action with Python
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.
Traffic camera view in QT web browser
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.         

Monday, April 2, 2018

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)
    N=float(a/math.sqrt(1-e2*(math.sin(math.radians(abs(lat)))**2)))
    X=(N+h)*math.cos(math.radians(lat))*math.cos(math.radians(lon))
    Y=(N+h)*math.cos(math.radians(lat))*math.sin(math.radians(lon))
    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)
plt.subplots_adjust(hspace=0.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.
 

Wednesday, March 21, 2018

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.
Select by Location and Select by Expression 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.
Crimes, city park and street dataset of Vancouver
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.

Query option to select all crimes event in the city park
Figure 3. Query option to select all crimes event in the city park
 Spatial query result Crimes within 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.
Spatial query to find any park where crime took place in it
Figure 5. Spatial query to find any park where crime took place

The selected park where crimes took place in it
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.

Spatial query 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.
Make a definition query
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.

Arterial road buffer
Figure 9. Arterial road buffer
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.
Crime event 100 m from main 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.   

Sunday, March 11, 2018

This is How to Add Google Maps Layers in QGIS 3

Finally QGIS 3 was released and I am so pleased for that. QGIS 3.0 Girona brings many features and updates such as improved geometry editing, symbols, layout, 3D layer, etc. What about plugins? I found many plugins already ported to QGIS 3, but unfortunately not all, including one of my favorite OpenLayers plugin, which can be used to add Google Maps basemap into QGIS map canvas such as Google Satellite, terrain and road map. What could we do to add those Google Maps layers in QGIS 3.0? Is there any other way around so we can still use Google Map in QGIS? The answer is Yes and I will show it in this tutorial.

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:
  • Google Maps: https://mt1.google.com/vt/lyrs=r&x={x}&y={y}&z={z}
  • Google Satellite: http://www.google.cn/maps/vt?lyrs=s@189&gl=cn&x={x}&y={y}&z={z}
  • Google Satellite Hybrid: https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}
  • Google Terrain: https://mt1.google.com/vt/lyrs=t&x={x}&y={y}&z={z}
  • Google Roads: https://mt1.google.com/vt/lyrs=h&x={x}&y={y}&z={z}  

Adding Google Maps Layers in QGIS 3

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.

Make a new connection in XYZ Tiles
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.

XYZ Tiles connection details
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.

Google Maps Layers in QGIS
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.

Google Terrain Layer in QGIS 3
Figure 4. Google Terrain Layer in QGIS 3

Google Satellite 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 :)

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.

Thursday, February 22, 2018

Digital Elevation Model (DEM) 3D Visualization in QGIS

QGIS 3.0 comes with 3D layer view which enable us to visualize GIS data in 3D. So we can get a more vivid visualization of a data which contains elevation or height information such as Digital Elevation Model (DEM). DEM contains elevation information of earth's terrain surface. It will be really cool to visualize DEM data in 3D to see terrain elevation difference of earth's surface as in figure 1, rather than seeing it on a flat surface. This tutorial will guide you how to visualize DEM in 3D with QGIS 3.0. In this tutorial I'm using QGIS 3.0 master candidate which can be downloaded from QGIS weekly download.
DEM 3D Visualization
Figure 1. DEM 3D Visualization

Download DEM Data

SRTM is a very popular DEM data source which can be downloaded from USGS Earth Explorer or NASA Earth Data server. But in QGIS you can directly download SRTM DEM data using SRTM Downloader Plugin. You just need to add a data into QGIS map canvas then zoom in to an area that you want to get a DEM data. Read this tutorial How to Download SRTM DEM Data in QGIS.

DEM 3D Visualization

After downloaded a DEM Data. Let's start to visualize it in 3D.

1. Add DEM data into QGIS map canvas as in figure 2.
DEM data in QGIS Map Canvas
Figure 2. DEM data in QGIS Map Canvas
 2. The data contains elevation information in each pixel. In figure 2, the DEM layer rendered in gray singleband. Let's make it more beautiful by changing its render type. To do this, open the DEM layer properties. In Style Option change Render type to Singleband pseudocolor. Then select a color ramp. The DEM data now seen like figure 4.

DEM Layer Properties
Figure 3. DEM Layer Properties
DEM Data in Singleband Pseudocolor
Figure 4. DEM Data in Singleband Pseudocolor
3. To visualize the DEM data in 3D. We do it in 3D Map View. You can find it in menu View as in figure 5.

3D Map View Menu
Figure 5. 3D Map View Menu
 4. Suddenly might be you get an error like figure 6. That's because the 3D map view does not support unprojected coordinate system. In this case the DEM data has geographic coordinate system. We don't have to change the DEM data coordinate system, but just change the map view coordinate system. To do this click the right bottom QGIS window with current CRS. See figure 7

3D Map View Error
Figure 6. 3D Map View Error
QGIS Map Canvas with Geographic Coordinate System
Figure 7. QGIS Map Canvas with Geographic Coordinate System (EPSG:4326)
5. Soon will appear the properties window. Change the coordinate system into a projected coordinate system like Mercator or Universal Traverse Mercator (UTM). I chose World Mercator, EPSG:3395.

Change a Project Coordinate System
Figure 8. Change a Project Coordinate System
6. After changing coordinate system. Open again the New 3D Map View. The 3D map view window will appear as in figure 9. Hmm...It's flat, I can't see any 3D view. Don't worry try to rotate the DEM data with holding shift button on the keyboard + left mouse click. Move the mouse, you will see the DEM data rotate as in figure 10.
DEM data in 3D Map View
Figure 9. DEM data in 3D Map View
DEM Data in 3D View after rotation
Figure 10. DEM Data in 3D View after rotation
 7. Looks promising. But still I can't see a nice terrain surface. It's because the elevation terrain is not set.To do it, click Configure button. The configuration window will appear as in figure 11. Set Elevation to the DEM data, and change Vertical scale to some number, for example 2. Vertical scale is a scale to raise the terrain's elevation. Higher scale will raise the terrain higher. Then you will see the terrain as in figure 12. Play with it by zoom in/out and rotate. 

3D Map View Configuration
Figure 11. 3D Map View Configuration
DEM 3D Visualization
Figure 12. DEM 3D Visualization
8. Now let's add a satellite imagery layer. For exampe I added world satellite imagery layer from ESRI web map server. Read how to add satellite imagery in QGIS.  QGIS will automatically overlay the satellite imagery above the DEM data as in figure 13. Again play with it.
DEM 3D with satellite imagery overlay
Figure 13. DEM 3D with satellite imagery overlay

That's all the tutorial how to visualize DEM data in 3D with QGIS. In the video below you can see how to do it from the beginning.