Search code examples
pythonmatplotlibcartopyshapley

Cartopy Splitting Large Circle

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.

split

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()```
Solution

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

    tissot32

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

    tissot56