# This software is open source software available under the BSD-3 license.
#
# Copyright (c) 2022 Triad National Security, LLC. All rights reserved.
# Copyright (c) 2022 Lawrence Livermore National Security, LLC. All rights
# reserved.
# Copyright (c) 2022 UT-Battelle, LLC. All rights reserved.
#
# Additional copyright and license information can be found in the LICENSE file
# distributed with this code, or at
# https://raw.githubusercontent.com/MPAS-Dev/MPAS-Analysis/main/LICENSE
"""
Functions for plotting inset maps in plots (e.g. for transects)
"""
# Authors
# -------
# Xylar Asay-Davis
import matplotlib.path
import matplotlib.ticker as mticker
import cartopy
import cartopy.crs as ccrs
import numpy
import shapely.geometry
from geometric_features.plot import subdivide_geom
[docs]
def add_inset(fig, fc, latlonbuffer=45., polarbuffer=5., width=1.0,
height=1.0, lowerleft=None, xbuffer=None, ybuffer=None,
maxlength=1.):
"""
Plots an inset map showing the location of a transect or polygon. Shapes
are plotted on a polar grid if they are entirely poleward of +/-50 deg.
latitude and with a lat/lon grid if not.
Parameters
----------
fig : ``matplotlib.figure.Figure``
A matplotlib figure to add the inset to
fc : ``geometric_features.FeatureCollection``
A collection of regions, transects and/or points to plot in the inset
latlonbuffer : float, optional
The number of degrees lat/lon to use as a buffer around the shape(s)
to plot if a lat/lon plot is used.
polarbuffer : float, optional
The number of degrees latitude to use as a buffer equatorward of the
shape(s) in polar plots
width, height : float, optional
width and height in inches of the inset
lowerleft : pair of floats, optional
the location of the lower left corner of the axis in inches, default
puts the inset in the upper right corner of ``fig``.
xbuffer, ybuffer : float, optional
right and top buffers from the top-right corner (in inches) if
lowerleft is ``None``.
maxlength : float or ``None``, optional
Any segments longer than maxlength will be subdivided in the plot to
ensure curvature. If ``None``, no subdivision is performed.
Returns
-------
inset : ``matplotlib.axes.Axes``
The new inset axis
"""
# Authors
# -------
# Xylar Asay-Davis
minLon, minLat, maxLon, maxLat = _get_bounds(fc)
figsize = fig.get_size_inches()
width /= figsize[0]
height /= figsize[1]
if lowerleft is None:
if xbuffer is None:
xbuffer = 0.1*width
else:
xbuffer /= figsize[0]
if ybuffer is None:
ybuffer = xbuffer*figsize[0]/figsize[1]
else:
ybuffer /= figsize[1]
lowerleft = [1.0 - width - xbuffer, 1.0 - height - ybuffer]
else:
lowerleft = [lowerleft[0]/figsize[0], lowerleft[1]/figsize[1]]
bounds = [lowerleft[0], lowerleft[1], width, height]
if maxLat <= -50:
# an Antarctic-focused map makes the most sense
inset = fig.add_axes(bounds,
projection=ccrs.SouthPolarStereo())
extent = [-180., 180., -90., max(-65., maxLat+polarbuffer)]
_set_circular_boundary(inset)
xlocator = mticker.FixedLocator(numpy.linspace(-180., 180., 9))
ylocator = mticker.FixedLocator(numpy.linspace(-90., -50., 9))
elif minLat >= 50:
# an Arctic-focused map makes the most sense
inset = fig.add_axes(bounds,
projection=ccrs.NorthPolarStereo())
extent = [-180, 180, min(65., minLat-polarbuffer), 90]
_set_circular_boundary(inset)
xlocator = mticker.FixedLocator(numpy.linspace(-180., 180., 9))
ylocator = mticker.FixedLocator(numpy.linspace(50., 90., 9))
else:
inset = fig.add_axes(bounds,
projection=ccrs.PlateCarree())
extent = [max(-180., minLon-latlonbuffer),
min(180., maxLon+latlonbuffer),
max(-90., minLat-latlonbuffer),
min(90., maxLat+latlonbuffer)]
xlocator = None
ylocator = None
# kind of like "top" justified -- graphics are toward the "north" end of
# the subplot
inset.set_anchor('N')
inset.set_extent(extent, ccrs.PlateCarree())
inset.add_feature(cartopy.feature.LAND, zorder=1)
inset.add_feature(cartopy.feature.OCEAN, zorder=0)
gl = inset.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
linewidth=0.5, color='gray', alpha=0.5,
linestyle='--')
if xlocator is not None:
gl.xlocator = xlocator
if ylocator is not None:
gl.ylocator = ylocator
for feature in fc.features:
geomtype = feature['geometry']['type']
shape = shapely.geometry.shape(feature['geometry'])
if maxlength is not None:
shape = subdivide_geom(shape, shape.geom_type, maxlength)
if geomtype in ['Polygon', 'MultiPolygon']:
inset.add_geometries((shape,), crs=ccrs.PlateCarree(),
edgecolor='blue', facecolor='blue', alpha=0.4,
linewidth=1.)
elif geomtype in ['Point', 'MultiPoint']:
point_x, point_y = shape.xy
inset.scatter(point_x, point_y, s=9, color='k',
transform=ccrs.PlateCarree(), edgecolors='face')
else:
inset.add_geometries((shape,), crs=ccrs.PlateCarree(),
edgecolor='k', facecolor='none', alpha=1.,
linewidth=1.)
# put a red point at the beginning and a blue point at the end
# of the transect to help show the orientation
begin = shape.coords[0]
end = shape.coords[-1]
inset.plot(begin[0], begin[1], color='r', marker='o',
markersize=3., transform=ccrs.PlateCarree())
inset.plot(end[0], end[1], color='g', marker='o',
markersize=3., transform=ccrs.PlateCarree())
return inset
def _set_circular_boundary(ax):
"""Set the boundary of the given axis to be circular (for a polar plot)"""
# Compute a circle in axes coordinates, which we can use as a boundary
# for the map. We can pan/zoom as much as we like - the boundary will be
# permanently circular.
theta = numpy.linspace(0, 2*numpy.pi, 100)
center = numpy.array([0.5, 0.5])
radius = 0.5
verts = numpy.vstack([numpy.sin(theta), numpy.cos(theta)]).T
circle = matplotlib.path.Path(verts * radius + center)
ax.set_boundary(circle, transform=ax.transAxes)
def _get_bounds(fc):
"""Compute the lon/lat bounding box for all transects and regions"""
bounds = shapely.geometry.GeometryCollection()
for feature in fc.features:
shape = shapely.geometry.shape(feature['geometry'])
shape_bounds = shapely.geometry.box(*shape.bounds)
bounds = shapely.geometry.box(*bounds.union(shape_bounds).bounds)
return bounds.bounds