Search code examples
pythonmatplotliblinecartopy

Is there a way to plot a line that is a consistent length using matplotlib/cartopy?


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!

Example of a line that the cross section would be through


Solution

  • 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()
    

    enter image description here