I am trying to create a "zoom" effect between two cartopy subplots (similar idea to the matplotlib mark_inset
function), where the main plot is Orthographic
and the zoomed in subplot is PlateCarree
. So far, I have figured out how to connect them with the patches.ConnectionPatch()
function and to mark the subregion using patches.Polygon()
.
I also want to draw a shaded polygon to cover the area of the plot between the two connection lines & two subplots as well. However, I can't seem to get the polygon to correctly cover the entire area between the two axes (see linked figure)... I suspect the issue might be in the transformation between the data coordinates and the projection coordinates? Can someone help me figure out where I've gone wrong? Thanks!
Function to add connection patches + polygons:
def custom_mark_zoom(axA, axB, direction='right', extent=None, fc=None, ec='k', alpha=1, transform=None):
# starting point:
# https://stackoverflow.com/questions/51268493/drawing-filled-shapes-between-different-axes-in-matplotlib
import matplotlib.patches as patches
import numpy as np
import matplotlib as mpl
xx = [extent[0], extent[1]]
yy = [extent[2], extent[3]]
xy = (xx[0], yy[0])
width = xx[1] - xx[0]
height = yy[1] - yy[0]
xyB1 = (0,1)
xyB2 = (0,0)
xyA1 = transform.transform_point(60,-40,ccrs.PlateCarree())
xyA2 = transform.transform_point(90,-40,ccrs.PlateCarree())
coordsA='data'
coordsB='axes fraction'
# First mark the patch in the main axes
pp = axA.add_patch(patches.Rectangle(xy, width, height, fc=fc, ec=ec, zorder=5, alpha=alpha, transform=ccrs.PlateCarree()))
# now draw an anchor line to the zoomed in axis
p1 = axA.add_patch(patches.ConnectionPatch(
xyA=xyA1, xyB=xyB1,
coordsA=coordsA, coordsB=coordsB,
axesA=axA, axesB=axB))
# draw a 2nd anchor line to the zoomed in axes
p2 = axA.add_patch(patches.ConnectionPatch(
xyA=xyA2, xyB=xyB2,
coordsA=coordsA, coordsB=coordsB,
axesA=axA, axesB=axB))
# create polygon between the lines to add shading (creates a zoomed effect)
verts = np.concatenate([p1.get_path().vertices,np.flipud(p2.get_path().vertices),])
p3 = axA.add_patch(mpl.patches.Polygon(verts, fc=fc, clip_on=False, alpha=alpha))
return pp, p1, p2, p3
Code to produce figure:
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
# set up figure
fig = plt.figure(figsize=(35,20))
gs = fig.add_gridspec(6, 6)
# Increase resolution of projection - needed to draw polygons accurately
map_proj = ccrs.Orthographic(central_latitude=-90.0, central_longitude=0)
map_proj._threshold /= 100
ax = fig.add_subplot(gs[0:4, 0:2], projection=map_proj)
ax2 = fig.add_subplot(gs[0:3, 2:6], projection=ccrs.PlateCarree())
# define extent
extent = [60, 90, -54, -40]
# format ax1
gl = ax.gridlines(draw_labels=True,
lw=3,
color="silver",
y_inline=True,
xlocs=range(-180,180,30),
ylocs=range(-80,91,10),
zorder=3,
)
ax.set_extent([-180,180,-90,-30], crs=ccrs.PlateCarree())
def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
see: https://stackoverflow.com/questions/4913349/haversine-formula-in-python-bearing-and-distance-between-two-gps-points
"""
from math import radians, cos, sin, asin, sqrt
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a))
r = 6371 # Radius of earth in kilometers. Use 3956 for miles
return c * r
# Create a circular map frame bounded by radius of map extent:
# see https://stackoverflow.com/questions/67877552/how-to-zoom-into-a-specific-latitude-in-cartopy-crs-orthographic
import matplotlib.path as mpath
r_limit = haversine(0,-90,0,-46)*1000 # convert km to m for CartoPy
circle_path = mpath.Path.unit_circle()
circle_path = mpath.Path(circle_path.vertices.copy() * r_limit,
circle_path.codes.copy())
# set circle boundary & extent
ax.set_boundary(circle_path)
# format ax2
gl = ax2.gridlines(draw_labels=True)
gl.top_labels = False
gl.right_labels = False
ax2.set_extent(extent)
# add the connection lines and shading
pp, p1, p2, p3 = custom_mark_zoom(ax, ax2, direction='right', extent=[60, 90, -54, -40], fc='lightgray', alpha=0.5, transform=map_proj)