Creating a geodetic line that goes around the whole ellipsoid with pyproj

Question

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":

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

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

Answer 1

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.