# great circles in QGIS and export in 3857 webmap

My plan is to visualize great circles in a leaflet webmap.

So coming from this great postgis article and R version I thought it would work like this in QGIS:

• create line from two points,
• store polyline in EPSG:4326
• add many vertexes using the densify algorithm
• store the densified line as 4326
• put it into leaflet and you should have nice "curves" on a map

what comes out is this:

And here a comparison of densified and normal lines in QGIS where you can see some small curvature of lines:

Is it just a problem with the EPSG 3857 or am I interpreting the post from @underdark wrong?

• Am I right, that I have to store all layers in EPSG:53027 instead of 4326 when using densify? – Riccardo Feb 4 '15 at 9:40

To create a great circle, you need either a custom azimuthal equidistant (aeqd) projection or a gnomonic (gnom) projection on one of your points.

In those projections, the great circle is a straight line, which you can densify to have a curved line in other projections. If the great circle crosses the 180° E/W line or the poles, it might be useful to cut the line on that line or point to avoid misplaced parts of the line in other projections.

updated workflow

1. create a text file with your start coordinates

E N

7 51

1. In QGIS, create a custom aeqd projection around your start point

+proj=aeqd +lat_0=51 +lon_0=7 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs

1. Load the text file as delimited text with WGS84 as CRS

2. Add the target points from any other backgrounds

3. switch the project CRS to your custom CRS

4. Create a new shapefile of type `line` and your custom CRS

5. Set snapping to your points layer with 10 pixels tolerance

6. connect your start point with all targets:

1. Densify the lines by 99 vertices

2. Switch the project CRS to EPSG:3857

Note that equidistant conic is a different projection, which does not show the great circles as straight lines:

A visualization of great circles in different projections can be found here:

http://www.progonos.com/furuti/MapProj/Normal/CartProp/Geodesic/geodesic.html

• EPSG:53027 - Sphere_Equidistant_Conic : so looks fine, right? – Riccardo Feb 4 '15 at 10:19
• No, that's a different projection: georeference.org/doc/equidistant_conic.htm – AndreJ Feb 4 '15 at 10:20
• thanks for clarification. greate example!... my problem is, that routes are starting at different points. defining only one CRS will not be enough ;-) – Riccardo Feb 4 '15 at 11:46
• A gnomonic projection can overcome that, but has problems with points on the other side of the world. So you have to repeat the process for other start points. You can collect all lines in a WGS84 layer at the end. – AndreJ Feb 4 '15 at 12:15
• A possible shortcut is that the source point will always have the coordinates 0,0 in a custom aeqd that uses its lat,lon for the center point, so you don't have to project the source points, just project the target points. – mkennedy Feb 4 '15 at 17:57

How about calculating linestrings of the true geodesics, see my answer here: How to make Great Circle Arcs which look good on a Web Mercator map? (with Python code example) and project it to whatever projection you need. The result on geojson.io looks like this:

in my special case I needed to create a custom crs for each point in the list of points in my point layer. it does now what it should: here is my code so far:

``````import processing
canvas = qgis.utils.iface.mapCanvas()
allLayers = canvas.layers()
point_layer = allLayers[0] #this is my pointlayer in my project
point_layer_dp=point_layer.dataProvider()
vlayer_crs=point_layer_dp.crs()
vertexes = []
features = point_layer.getFeatures()
for f in features:
geom = f.geometry()
point = geom.asPoint()
vertexes.append(point)

crs_src = QgsCoordinateReferenceSystem(vlayer_crs)
for i in range(0, len(vertexes)-1):
crs = QgsCoordinateReferenceSystem()
crs.createFromProj4("+proj=aeqd +lat_0=" + str(vertexes[i][1])+ " +lon_0=" + str(vertexes[i][0]) +" +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs")
#crs.saveAsUserCRS('azimuthal equidistant point' + str(i+1))
qgis.utils.iface.mapCanvas().mapRenderer().setDestinationCrs(crs)
layer = QgsVectorLayer('LineString?crs=' + crs.toWkt(), 'line'+str(i), "memory")
layer.setCrs(crs)
pr = layer.dataProvider()
line = QgsFeature()
xform = QgsCoordinateTransform(crs_src, crs)
start_point = QgsPoint(0,0)
end_point = xform.transform(QgsPoint(vertexes[i+1]))
seg = [start_point, end_point]
line.setGeometry(QgsGeometry.fromPolyline(seg))