Cartopy Splitting Large Circle

Question

I am trying to plot the visible region of the sky from Denver, CO, and show this region as a shaded "circle" on a map projection.

I'm expecting a distorted circle (because of map projections) on the map, centered on Denver. When I try to plot the shapely circle with fill, it appears to be centered correctly but it splits along a latitude, as seen in the graphic. I cannot figure out why it does this and how to correct it. The desired result is that the region above the curve (all North America and ocean above Russia) are shaded red.

import numpy as np
import cartopy.geodesic as cgeod
import cartopy.crs as crs
import cartopy.feature as cfeature
import shapely
import matplotlib.pyplot as plt

# SET CONSTANTS
RE = 6371008.8 # meters, WGS84 Sphere
h = 35786000.0 # meters, GEO

def arc_dist_to_FOV(h, ah):
    a = (np.pi/180)*(90 + ah) # radians
    b = np.arcsin(RE * np.sin(a)/(RE + h)) # radians
    g = np.pi - a - b # radians
    d = (180/np.pi)*g*RE
    return d

cm = 0
proj = crs.PlateCarree(central_longitude=cm)
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(1,1,1, projection=proj)

ax.set_global()

ax.add_feature(cfeature.COASTLINE, edgecolor='black')
ax.add_feature(cfeature.BORDERS, edgecolor='black')
ax.add_feature(cfeature.STATES, edgecolor='black')

ax.stock_img()
# ax.background_img(name='BM', resolution='high')
ax.gridlines(draw_labels=True, crs=proj)
lat, lon = 39.7392, -104.985

plt.scatter(x=lon, y=lat, color='blue', s=10, transform=proj)

# COMPUTE FOV
ah = 20 # degrees above horizon
dFOV = arc_dist_to_FOV(h, ah)
site_lat = lat
site_lon = lon
    
# ADD SHAPES TO MAP
circle_points = cgeod.Geodesic().circle(lon=site_lon, lat=site_lat, radius=dFOV, n_samples=1000, endpoint=False)
geom = shapely.geometry.Polygon(circle_points)
ax.add_geometries((geom,), crs=proj, alpha=0.5, facecolor='red', edgecolor='black', linewidth=1)

# must save BEFORE show cmd
# plt.savefig('name.png', bbox_inches='tight', dpi=300)

plt.show()```

You'll have to insert the points (-180, 90) and (180, 90) at the right place of your polygon. Maybe there's a library function for that in cartopy or another plotting library that handles this. — Joooeey, Apr 23 '23 at 20:11
or perhaps it's simpler to skip shapely and use plt.fill_between with the coordinates and a "line" at the north pole https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.fill_between.html — Joooeey, Apr 23 '23 at 20:17
Thanks @Joooeey, I was going to pursue that until @swatchi pointed out the `.tissot()` method. — DoomedJupiter, Apr 25 '23 at 02:06
For anyone following, there is an error in my code above. The fourth line in the function `arc_dist_to_FOV` incudes the term `(180/np.pi)`, which converts `g` to degrees. To compute the arc length of a circle the angle value remain in radians, not degrees. This produces obviously incorrect results for values other than `ah = 20`. I did not correct the code since it produces the included graphic. — DoomedJupiter, Apr 25 '23 at 02:08

swatchai · Accepted Answer · 2023-04-24T01:05:57.263

For a (projection of) circle plotting on any projection, the method .tissot() is available. That simplifies drawing circles on the map greatly.

The (projection of) big circle in the plot below makes use of this line of code:

ax.tissot(rad_km=dFOV/1000, lons=[site_lon,], lats=[site_lat,], n_samples=128, fc="red", ec="black", alpha=0.75)

With a view of Orthographic projection, the plot will be similar to this:

Cartopy Splitting Large Circle

1 Answers1