Digital Geography

10. February 2015

Creating ARCs in QGIS: The Python Way

So I came across this nice little project which focussed on trip planning and route over large distances… And there was this nice little post from Nathan Yau at flowingdata.com where he describes the making of great circles from one point to different other points in R and the other example from Anita Graser where she shows how to deal with an Arc in QGIS but using postgis functionality. So what about QGIS itself and a programmatic way? See yourself…Andre Joost pointed me into the right direction by describing the process with a nice example after I asked on gis.stackexchange.com . I extracted the information and would like to show you the process in Python/PyQGIS

Creating a Line with Python

The connection and creation of an Arc or great circles starts with a starting and endpoint. Let’s take it for first try as these are given in a pointshapefile. Lets extract the CRS and get the coordinates of those two points:
point_layer = qgis.utils.iface.mapCanvas().currentLayer()
point_layer_dp=point_layer.dataProvider() #we need this to extract the coordinate system
source_crs=point_layer_dp.crs() #so lets extract it
features = point_layer.getFeatures() #now determine each feature
vertexes = [] #here we will store the coordinates as I like arrays
for f in features:
	geom = f.geometry() #access to the geometry of a point
	point = geom.asPoint() #here we extract the coordinates as one point
	vertexes.append(point) #and store it in the array vertexes by appending the vertexes array
Now we have the point coordinates in a nice little array which is iterable by an index. As I just have two points in this file we can easily create a line connecting these two points in the crs defined by the points:
layer = QgsVectorLayer('LineString', 'connecting line', "memory") #this will create an in memory layer which is not stored on a HDD
layer.setCrs(source_crs) #we set the CRS to the one provided by the point layer
pr = layer.dataProvider() #once again we need the provider (should be shapefile/ESRI)
line = QgsFeature() #let's create an empty Feature which will be the box to put our segment(s) into
start_point = vertexes[0] #this is the start of the first segment
end_point = vertexes[1] #and the end of course
seg = [start_point, end_point] #so this is our whole segment
line.setGeometry(QgsGeometry.fromPolyline(seg)) #we will set the geometry of our so far empty feature from this array of points
pr.addFeatures([line]) #and add it to the layer
layer.updateExtents() #update it
QgsMapLayerRegistry.instance().addMapLayer(layer) #and put it on the map. Store it manually if you need to.
The result might look something like this:

connection with line from one point to another

Custom CRS for a Real Great Circle/Arc

What we will need for real great circles in a GIS is something that preserves distances and angles. Andre pointed out to use a azimuthal equidistant projection which needs to be defined for a custom starting point. The CRS for our starting point of the line would look something like this:
crs_new = QgsCoordinateReferenceSystem()
crs_new.createFromProj4("+proj=aeqd +lat_0=" + str(vertexes[0][1])+ " +lon_0=" + str(vertexes[0][0]) +" +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs") #which uses the Lat and Lon values of our starting point
But before we can do this, we need to get the coordinates of the line segments. According to this new CRS the start point will always get the coordinate (0,0) but we need to find the target coordinates for the endpoint. Therefore we will transform coordinates in Python:
crs_src = QgsCoordinateReferenceSystem(source_crs)
xform = QgsCoordinateTransform(crs_src, crs_new)
start_point_new = QgsPoint(0,0)
end_point_new = xform.transform(QgsPoint(vertexes[1]))
seg_new = [start_point_new, end_point_new]
So now we have the start and end points stored in a new segment but with transformed coordinates in a very special CRS. As we want a smooth line we will not come across to insert vertexes in this new line. The main problem now is to set the current project CRS in QGIS to the crs from the new created line. If we have done this, we can call an algroithm defined in QGIS which densifies the line with new vertexes in a given distance:
crs_old = qgis.utils.iface.mapCanvas().mapRenderer().destinationCrs() #just store it for later 😉
qgis.utils.iface.mapCanvas().mapRenderer().setDestinationCrs(crs_new) #change current crs to fit with the desired one.
layer_new = QgsVectorLayer('LineString', 'line_new CRS', "memory") #once more...
layer_new.setCrs(crs_new) #and again
pr_new = layer_new.dataProvider() #yap, 2nd time
line_new = QgsFeature() #and...
line_new.setGeometry(QgsGeometry.fromPolyline(seg_new)) #and...
pr_new.addFeatures([line_new]) #okay now we get serious
layer_new.updateExtents() #whoha
QgsMapLayerRegistry.instance().addMapLayer(layer_new) #now we have a line with new coordinates
import processing #we need this library to use the QGIS processing tools
dens_layer = processing.runalg("qgis:densifygeometriesgivenaninterval",layer_new,'100000',None) #which densifies the line by new vertexes each 100 km
dens_layer=QgsVectorLayer(dens_layer.get('OUTPUT'), "densified_layer", "ogr") #and we would like to store the result in a new layer
QgsMapLayerRegistry.instance().addMapLayer(dens_layer) #which we put on the map
QgsMapLayerRegistry.instance().removeMapLayer(layer_new.id()) #we don't need the old one anymore
qgis.utils.iface.mapCanvas().mapRenderer().setDestinationCrs(crs_old) #and change back the CRS to the old one.

what a nice curve

If you want to use this in another project zou need to store this separately.

Difficulties

With the above situation it is not possible to create nice lines when you cross -180/180 longitude:

crossing the date line…

So how to overcome this? Therefore we retransform the densified layer into a layer based on EPSG 4326 and extract the vertices of this densified layer. We will identify the vertex which is near the border and split the line into two features at this vertex (sorry, no comments from now on!):
crs_4326 = QgsCoordinateReferenceSystem(4326)
xform2 = QgsCoordinateTransform(crs_new, crs_4326)
layer_4326 = QgsVectorLayer('LineString', 'line_4326', "memory")
layer_4326.setCrs(crs_4326)
pr_4326 = layer_4326.dataProvider()
outFeat = QgsFeature()
for f in dens_layer.getFeatures():
	geom = f.geometry()
	geom.transform(xform2)
	outFeat.setGeometry(geom)
	outFeat.setAttributes(f.attributes())

pr_4326.addFeatures([outFeat])
layer_4326.updateExtents()
QgsMapLayerRegistry.instance().addMapLayer(layer_4326)
vertexes_merge = []
for f in layer_4326.getFeatures():
	line = f.geometry()
	for i in line.asPolyline():
		vertexes_merge.append(QgsPoint(i))

import math
breaks = []
for i in range(0,len(vertexes_merge)-1):
	if math.sqrt((vertexes_merge[i][0]-vertexes_merge[i+1][0])**2) > 2:
		breaks.append(i)

breaks.append(len(vertexes_merge)-1)
layer_merge = QgsVectorLayer('LineString', 'line', "memory")
layer_merge.setCrs(crs_4326)
pr_merge = layer_merge.dataProvider()
start = 0
for index in range(0,len(breaks)):
	end = breaks[index]
	line_merge = QgsFeature()
	line_merge.setGeometry(QgsGeometry.fromPolyline(vertexes_merge[start:end]))
	start = breaks[index] + 1
	pr_merge.addFeatures([line_merge])
	layer_merge.updateExtents()

QgsMapLayerRegistry.instance().addMapLayer(layer_merge)
Now we will get a nicely splitted line with two features which can be exported using the qgis2leaf plugin and be shown in a webmap…

splitted line to avoid bad rendering

ps: this all very rough but fit the needs. You can download the sourcecode as python file here: two_points. Yet it is computational intensive if you work with a list of several points. See the current implementation status on github and the video here: