1

https://jswhit.github.io/pyproj/pyproj.Geod-class.html
Pyproj's Geod class has a npts method to create points on the geodesic between two points:

npts(self, lon1, lat1, lon2, lat2, npts, radians=False)

Given a single initial point and terminus point (specified by python floats lon1,lat1 and lon2,lat2), returns a list of longitude/latitude pairs describing npts equally spaced intermediate points along the geodesic between the initial and terminus points.

How can I use this to create points on a geodesic that goes around the whole ellipsoid? The npts method uses the shortest of the two possible paths, so you can only use it for a line spans half the ellipsoid.

My attempt so far:

I thought I could use the fwd method to build the line in multiple steps but I do not know how to calculate/find the azimuth to continue on once I created the first segment with npts. If I use the same azimuth, the result is wrong (see image below). fwd returns the back azimuth, while I think I would need the forward azimuth at the terminus point.

fwd(self, lons, lats, az, dist, radians=False)

forward transformation - Returns longitudes, latitudes and back azimuths of terminus points given longitudes (lons) and latitudes (lats) of initial points, plus forward azimuths (az) and distances (dist).

Here is my failed attempt:

from pyproj import Geod
from shapely.geometry import Point, LineString

geod = Geod("+ellps=WGS84")
slon, slat = 42, -23
azimuth = 230
distance = 14600000

# first half
h1lon, h1lat, _backazi = geod.fwd(slon, slat, azimuth, distance/2)
h1points = geod.npts(slon, slat, h1lon, h1lat, npts=100)

# other half
h2lon, h2lat, _backazi = geod.fwd(h1lon, h1lat, azimuth, distance/2)
h2points = geod.npts(h1lon, h1lat, h2lon, h2lat, npts=100)

points = [Point(slon, slat)] + h1points + [Point(h1lon, h1lat)] + h2points + [Point(h2lon, h2lat)]
line = LineString(points)

The line resulting from this looks like this, you can see where the second half starts, it's where the direction "skips":

broken geodesic

This is how it would look like if the second azimuth was correct (this was created as a single line):

correct geodesic

If we could solve that azimuth issue, creating the geodesic around the ellipsoid could simply be done step by step.

2

pyproj hooks up with the geodesic routines (I am the author of these) in proj.

The geod_directline and geod_inverseline functions both initialize objects which allow you to compute points on a geodesic encircling the earth as many times as you like. So either implement calls in pyproj to invoke this functionality or request that the maintainer do so.

Alternative you can use my native implementation of the geodesic routines in python. The documentation includes an example of computing waypoints.

  • Thanks, I will try that! Had no idea that geographiclib had Python bindings :) – bugmenot123 Nov 18 '16 at 9:15
  • This worked great but I am not sure if my solution is the best way. To get an initial line I used Direct to create a point some kilometers infront of my starting point. Is there a way to put points on a line with just the starting point and the angle? Before I extend the answer with code I want to make it good. :) – bugmenot123 Nov 20 '16 at 20:41
  • You're making it too complicated for your application. Just use the Line method to define a geodesic in terms of initial position and azimuth. See geographiclib.sourceforge.net/html/python/… – cffk Nov 20 '16 at 21:03
  • Haha, I knew it would be easier. Should have read the docs. Thanks! – bugmenot123 Nov 21 '16 at 8:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.