Drawing Circles with cartopy in orthographic projection

Question

I am trying to draw circles at a given geographical coordinate with a certain radius using cartopy. I want to draw using an orthographic projection, which is centred at the centre of the circle.

I use the following python code for testing:

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

# example: draw circle with 45 degree radius around the North pole
lon = 0
lat = 90
r = 45

# find map ranges (with 5 degree margin)
minLon = lon - r - 5
maxLon = lon + r + 5
minLat = lat - r - 5
maxLat = lat + r + 5

# define image properties
width = 800
height = 800
dpi = 96
resolution = '50m'

# create figure
fig = plt.figure(figsize=(width / dpi, height / dpi), dpi=dpi)
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Orthographic(central_longitude=lon, central_latitude=lat))
ax.set_extent([minLon, maxLon, minLat, maxLat])
ax.imshow(np.tile(np.array([[cfeature.COLORS['water'] * 255]], dtype=np.uint8), [2, 2, 1]), origin='upper', transform=ccrs.PlateCarree(), extent=[-180, 180, -180, 180])
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', resolution, edgecolor='black', facecolor=cfeature.COLORS['land']))
ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_0_countries', resolution, edgecolor='black', facecolor='none'))
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'lakes', resolution, edgecolor='none', facecolor=cfeature.COLORS['water']), alpha=0.5)
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', resolution, edgecolor=cfeature.COLORS['water'], facecolor='none'))
ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_1_states_provinces_lines', resolution, edgecolor='gray', facecolor='none'))
ax.add_patch(mpatches.Circle(xy=[lon, lat], radius=r, color='red', alpha=0.3, transform=ccrs.PlateCarree(), zorder=30))
fig.tight_layout()
plt.savefig('CircleTest.png', dpi=dpi)

plt.show()

I get a correct result at the equator (set lat to 0 in example above):

But when I move towards a pole the shape is distorted (lat = 45):

At the pole I only see one quarter of the circle:

I would expect to always see a perfect circle in orthographic projection, if the view is centred correctly. I also tried to use a different transform in the add_patch method, but then the circle completely vanishes!

Answer 1

You approach of defining the circle in PlateCarree coordinates is not going to work, because this is a cartesian projection and a circle drawn on it is not necessarily circular in spherical geometry (unless the circle is at (0, 0) as you saw).

Since you want the result to be circular in the Orthographic projection, you could draw the circle in native coordinates. This requires first defining your Orthographic projection centred on the centre of your circle, then computing what the radius of the circle (which you specify in degrees) would be in projection coordinates (distance from the centre of the projection). Doing it this way is convenient because it also gives you a neat way of determining the correct map extents. The example below computes the radius in orthographic coordinates by transforming a point 45 degrees north (or south if more convenient) away from the centre of the projection and gives the following:

The full code is below:

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

# example: draw circle with 45 degree radius around the North pole
lat = 51.4198101
lon = -0.950854653584
r = 45

# Define the projection used to display the circle:
proj = ccrs.Orthographic(central_longitude=lon, central_latitude=lat)


def compute_radius(ortho, radius_degrees):
    phi1 = lat + radius_degrees if lat <= 0 else lat - radius_degrees
    _, y1 = ortho.transform_point(lon, phi1, ccrs.PlateCarree())
    return abs(y1)

# Compute the required radius in projection native coordinates:
r_ortho = compute_radius(proj, r)

# We can now compute the correct plot extents to have padding in degrees:
pad_radius = compute_radius(proj, r + 5)

# define image properties
width = 800
height = 800
dpi = 96
resolution = '50m'

# create figure
fig = plt.figure(figsize=(width / dpi, height / dpi), dpi=dpi)
ax = fig.add_subplot(1, 1, 1, projection=proj)
# Deliberately avoiding set_extent because it has some odd behaviour that causes
# errors for this case. However, since we already know our extents in native
# coordinates we can just use the lower-level set_xlim/set_ylim safely.
ax.set_xlim([-pad_radius, pad_radius])
ax.set_ylim([-pad_radius, pad_radius])
ax.imshow(np.tile(np.array([[cfeature.COLORS['water'] * 255]], dtype=np.uint8), [2, 2, 1]), origin='upper', transform=ccrs.PlateCarree(), extent=[-180, 180, -180, 180])
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', resolution, edgecolor='black', facecolor=cfeature.COLORS['land']))
ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_0_countries', resolution, edgecolor='black', facecolor='none'))
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'lakes', resolution, edgecolor='none', facecolor=cfeature.COLORS['water']), alpha=0.5)
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', resolution, edgecolor=cfeature.COLORS['water'], facecolor='none'))
ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_1_states_provinces_lines', resolution, edgecolor='gray', facecolor='none'))
ax.add_patch(mpatches.Circle(xy=[lon, lat], radius=r_ortho, color='red', alpha=0.3, transform=proj, zorder=30))
fig.tight_layout()
plt.savefig('CircleTest.png', dpi=dpi)
plt.show()

Answer 2

This might be a little late, but there is a convient function in Cartopy for this.

We can use Cartopy's .circle function (documentation) to generate a ring of points with a specified radius from a particular (longitude & latitude) in the Geodesic coordinate frame and then plot a polygon with those points using Shapely.

This would look something like the following

circle_points = cartopy.geodesic.Geodesic().circle(lon=lon, lat=lat, radius=radius_in_meters, n_samples=n_points, endpoint=False)
geom = shapely.geometry.Polygon(circle_points)
ax.add_geometries((geom,), crs=cartopy.crs.PlateCarree(), facecolor='red', edgecolor='none', linewidth=0)

Specifying the crs as PlateCarree does not matter and merely avoids a warning with Shapely. You will keep your desired projection. However, if you are plotting directly with the circle center on the pole, you might still have an issue and may need to do some fancy transformations (Haven't tested it recently, but recall from a few months ago it being a little wonky).

You could also manually compute these points using the pyproj library Cartopy makes use of, specifically the Geod class. Pick a point with a radius and loop through the azmoths for however fine you want your circle to be with the .inv or .fwd function similar to the suggestion in https://stackoverflow.com/a/57002776/2430454. I don't recommend this method, but used it a long while back to accomplish the same thing.