I am using matplotlib and cartopy to draw lines overlaid on maps in python. As of right now, I am just identifying the lat/lon of two points and plotting a line between them. Since I am taking cross sections across these lines, I would like to find a way to make the line the same length (say 300km long) no matter where I place it on the map. Is this possible without just using trial and error and setting the points until they are the desired length?
lat1, lat2, lon1, lon2 = [34.5, 36, -100, -97] x, y = [lon1, lon2], [lat1, lat2]
ax1.plot(x, y, color="black", marker="o", zorder=3, transform = ccrs.PlateCarree(), linewidth = 2.5)
Here are the relevant parts of the code that I am using now. This works, but I am looking for a way to hold the line length constant rather than changing the values for the endpoints at "lat1, lat2, lon1, lon2." I envision setting a line length, a mid-point (lat/lon), and an angle that would pivot around that point. I don't know if that's even possible, but that's how I'd imagine it'd have to work!
The Geod
class from pyproj is very convenient for these types of operations. The example below moves in a line from a given lat/lon based on an azimuth and distance (in meters).
Since you mention starting with a mid-point, you can do this twice for each direction to get the full line.
In the example below, I just retrieve the end points and use Cartopy (ccrs.Geodetic()
) to do the interpolation to a great circle. But you can also do this yourself with the same Geod object (see the npts method), and sample a given amount of points a long the line. The latter might be convenient if you need those coordinates to extract data for example.
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from pyproj import Geod
import numpy as np
# start with a random point
np.random.seed(0)
lon = np.random.randint(-180,180)
lat = np.random.randint(-90,90)
# and a random direction and distance
half_width = (1 + np.random.rand()) * 5000000 # meters
azimuth = np.random.randint(0,360)
# calculate the end points
geod = Geod(ellps="WGS84")
lon_end1, lat_end1, azi_rev1 = geod.fwd(lon, lat, azimuth, half_width)
lon_end2, lat_end2, azi_rev2 = geod.fwd(lon, lat, azimuth-180, half_width)
# visualize the result
fig, ax = plt.subplots(
figsize=(8,4), dpi=86, layout="compressed", facecolor="w",
subplot_kw=dict(projection=ccrs.PlateCarree()),
)
ax.set_title(f"{lon=}, {lat=}, {azimuth=}, {half_width=:1.1f}m")
ax.plot(
[lon_end1, lon, lon_end2],
[lat_end1, lat, lat_end2],
"ro-",
transform=ccrs.Geodetic(),
)
ax.coastlines()
ax.set_global()