Creating a Line with Python

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

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.

Custom CRS for a Real Great Circle/Arc

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

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]

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.

Difficulties

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)

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/PyQGISThe 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: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:The result might look something like this: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: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: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:If you want to use this in another project zou need to store this separately.With the above situation it is not possible to create nice lines when you cross -180/180 longitude: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!):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…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: