A profile shows elevation variation along a line in two dimensional space. We use it for many applications such as: to know the elevation variation for a street, getting to know a landscape height variation between two points, pipe line or sewer planning, and so on.
|
|
Figure 1. Elevation Profile |
Import Libraries
For this tutorial we are using QGIS processing and matplotlib library. Therefore we need to import those two libraries as following.
#IMPORT MODULE/LIBRARY import processing import matplotlib.pyplot as plt
Define Input/Output Data
Next, we need to define input and output data. The input data which are required to create elevation profile are elevation in raster format and a line in vector format. Both input data must be in the same projection and preferred in a projected coordinate system. The output data are coming from processing result that involved in the process which are points along the path line and points with extracted elevation data from elevation raster data. In this tutorial I'm using elevation data from LiDAR (If you are interested how to make DEM from LiDAR check out this tutorial) and a line path as shown in figure 2.
|
| Figure 2. Elevation and Line Data |
The following code is used to define input and output data. For input data, I used DEM data from lidar and a line path in geopackage file type. The last two lines are the output data, which are actually temporary data, that can be removed after the process is done.
#DEFINE INPUT/OUTPUT DATA #CHANGE THIS WITH YOUR PATH AND INPUT DATA project_path='E:/Project/' elevation_data=project_path+'lidar_dem.tif' line=project_path+'line_path.gpkg' line_points=project_path+'out_points.gpkg' points_z=project_path+'point_z.gpkg'
Processing and Data Conversion
There are two processing steps that involved in this process. The first one is to convert the path line into a series of points along the line with a specified interval. This process is done using v.to.points tool from GRASS module. The second one is to extract elevation information from elevation data for each points. For this, the Add raster values to points tool from SAGA is used. The output for each process is stored in respected output file that already defined before.
The code below shows the processing steps using each tool. To generate points
along the path line with an interval, we need an interval value. Actually we
can define any positive value for interval as long as it's not exceed the
length of the line. But for this tutorial, I used the x resolution of the
elevation data. To get the x resolution can be done with
rasterUnitPerPixel method from a Qgs raster layer. Therefore firstly we
need to define a raster layer using QgsRasterLayer class.
#DEFINE RASTER LAYER AND GET X RESOLUTION raster_layer=QgsRasterLayer(elevation_data,'ogr') interval=raster_layer.rasterUnitsPerPixelX() #PROCESSING USING SAGA AND GRASS processing.run("grass7:v.to.points", {'input':line,'use':1,'dmax':interval,'-i':True,'-t':False,'output':out_points}) processing.run("saga:addrastervaluestopoints", {'SHAPES':line_points,'GRIDS':[dem],'RESAMPLING':0,'RESULT':point_z})
Extract Elevation Data and Calculate Distance
After we get the points with it's elevation data. Next we we will extract each elevation value and calculate the distance for each vertex from start node to end node. The distance and elevation data will be dumped into a list which is called z_list and list_distance. The process will be done in a loop for each point starting at line 22 in the following code. But before that, we need to set up some pre-processing steps such as define a vector layer and getting it's first feature to get start node coordinate and it's elevation. Then Initialize list and dictionary to store two points (start and to point) and initialize a class variable using QgsDistanceArea class and measureLine method with two points in the dictionary for distance calculation. Keep in mind to define a datum or ellipsoid in the distance calculation, with the datum as in your input data. For this case I'm using elevation data with North American Datum 1983(NAD83).
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 |
#CALCULATE DISTANCE AND GETTING ELEVATION #DEFINE VECTOR LAYER z_layer=QgsVectorLayer(points_z,'ogr') #GETTING VECTOR FEATURES features=z_layer.getFeatures() first_feat=next(features) first_point=first_feat.geometry().asPoint() first_z=first_feat.attributes()[-1] #INITIALIZE DICTIONARY AND LIST point_dict={'from':'','to':''} point_dict['from']=first_point list_distance=[0] z_list=[first_z] #INITIALIZE CLASS FOR DISTANCE CALCULATION d=QgsDistanceArea() d.setEllipsoid('NAD83') #GETTING Z FOR EACH POINT for f in features: z=f.attributes()[-1] z_list.append(z) xy_to=f.geometry().asPoint() point_dict['to']=xy_to distance=d.measureLine(point_dict['from'],point_dict['to']) segment_line=distance+list_distance[-1] list_distance.append(segment_line) point_dict['from']=point_dict['to'] |
Plotting The Profile
After getting a distance and elevation for each point along the line. Now we
are ready to plot it. For plotting we are using matplotlib library which
already imported in the first step. The following code is using for profile
plotting.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
#PLOTTING PROFILE min_z=round(min(z_list),3) max_z=round(max(z_list),3) mean_z=round(sum(z_list)/len(z_list),3) d_max=list_distance[-1] plt.figure(figsize=(10,4)) plt.plot(list_distance,z_list) plt.plot([0,d_max],[min_z,min_z],'g--',label='Min. : '+str(min_z)) plt.plot([0,d_max],[max_z,max_z],'r--',label='Max. : '+str(max_z)) plt.plot([0,d_max],[mean_z,mean_z],'y--',label='Mean : '+str(mean_z)) plt.grid() plt.legend(loc=1) plt.xlabel("Distance(ft)") plt.ylabel("Elevation(ft)") plt.fill_between(list_distance,z_list,min_z,alpha=0.1) plt.show() |
That's all the tutorial how to create elevation profile using python in QGIS
3. If you want to explore more tutorials about QGIS Python programming, please
visit
QGIS Python Programming Tutorial Series
for other interesting topics. Thanks for reading and happy coding!
