1

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)

Figure with incorrect polygon:
enter image description here

Mr. T
  • 11,960
  • 10
  • 32
  • 54
Nabb96
  • 11
  • 2

0 Answers0