8

My plan is to visualize great circles in a Leaflet web map.

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

  • create a 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: leaflet map in EPSG 3857

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

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

Improve this question
1
  • Am I right, that I have to store all layers in EPSG:53027 instead of 4326 when using densify? Commented Feb 4, 2015 at 9:40

3 Answers 3

Reset to default
13

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:

enter image description here

  1. Densify the lines by 99 vertices

  2. Switch the project CRS to EPSG:3857

enter image description here

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

enter image description here

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

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

Improve this answer
5
  • EPSG:53027 - Sphere_Equidistant_Conic : so looks fine, right? Commented Feb 4, 2015 at 10:19
  • No, that's a different projection: georeference.org/doc/equidistant_conic.htm Commented Feb 4, 2015 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 ;-) Commented Feb 4, 2015 at 11:46
  • 1
    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. Commented Feb 4, 2015 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. Commented Feb 4, 2015 at 17:57
1

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:

enter image description here

Improve this answer
0

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: enter image description here 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)) pr.addFeatures([line]) layer.updateExtents() QgsMapLayerRegistry.instance().addMapLayer(layer) dens_layer = processing.runalg("qgis:densifygeometriesgivenaninterval",'line'+str(i),1000,None) vlayer=QgsVectorLayer(dens_layer.get('OUTPUT'), "densified_layer" + str(i), "ogr") QgsMapLayerRegistry.instance().addMapLayer(vlayer) qgis.utils.iface.mapCanvas().mapRenderer().setDestinationCrs(crs_src)

Now it's only missing to get the densified layers into EPSG 4326 and merge them... But this is prob. another story to tell.

Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.