Monday, 9 July 2018

coordinate system - How can I make a polyline wrap around the world?



I am using leaflet maps to create a representation of a round-the-world challenge. I would like to add a polyline which heads east from tokyo and then appears west of south america on the map - but instead I get a line which crosses the map in the opposite direction (see yellow line).


enter image description here


I think this is probably related to datelines and /or coordinate systems but am a bit sketchy on the detail. Can anyone explain the theory behind what I need to do to get this to work? I am using the bluemarble projection from Nasa:


var bluemarble = new L.TileLayer.WMS("http://demo.opengeo.org/geoserver/wms", {
layers: 'bluemarble',
attribution: "Data © NASA Blue Marble, image service by OpenGeo",
minZoom: 2,
maxZoom: 5
});

Answer




You need to break the polyline at the +-180 degree meridian. This requires finding the latitude at which the polyline crosses that meridian. Your GIS probably has methods to do the breaking. If not, a simple solution can be derived from code shown in a related thread. Here are some details.




  • A polyline is represented as a sequence of vertices, each given in (lat, lon) form, with -180 <= lon <= 180. You need to check each successive pair to see whether it crosses the +-180 meridian. There's a quick test: if the absolute value of the difference of longitudes is 180 or greater, there is a crossing.




  • Within each segment (lat0,lon0) --> (lat1,lon1) that crosses the +-180 meridian, you need to break the polyline into two pieces where it crosses.




The key is finding the latitude of the break point with reasonable accuracy. This is most easily done with a spherical earth model: the error (compared to a more accurate ellipsoidal model) will be too small to notice.



Let the segment in question go from point 0 at (lat0,lon0) to point 1 at (lat1,lon1). The break point can be found by running a straight line segment in 3D between the two points as represented in Cartesian coordinates and finding where the y coordinate is zero. The Cartesian coordinates are


(x0, y0, z0) = (cos(lon0)*sin(lat0), sin(lon0)*sin(lat0), cos(lat0))

and a similar expression giving (x1, y1, z1) for point 1. Solve the equation


t * y0 + (1-t) * y1 = 0

for t; that is,


t = y1 / (y1 - y0).

The coordinates of the intersection are therefore



(x, y, z) = (t * x0 + (1-t) * x1, 0, t * z0 + (1-t) * z1)

This point (which lies beneath the earth's surface somewhere beneath the +-180 meridian) has latitude equal to


lat2 = ATan(z/x).

The break point needs to be represented in two ways. When attaching it after (lat0, lon0) to terminate the first part of the broken polyline, use (lat2, -180) if lon0 is negative and otherwise use (lat2, 180). When attaching it before (lat1, lon1) to start the second part of the broken polyline, follow a similar rule.


In exceptional cases, one or both of point 0 and point 1 may be on the +-180 meridian. Following this procedure will cause you to place a zero-length segment on one of the polyline pieces you create. If this might cause a problem with the GIS, test for this condition.


Note that a polyline can cross this meridian more than once. Therefore, after finding the first break and breaking the polyline into two parts, you need to process the second part in the same way.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...