Limiting latitudinal extend of a cartopy orthographic projection
I am trying to plot a map of a sphere with an orthographic projection of the Northern (0-40N) and Southern (0-40S) hemispheres, and a Mollweide projection of the central latitudes (60N-60S). I get the following plot:
which shows a problem: there is a square bounding box with cut corners around the hemispherical plots. Note that the extent of the colours is the same for all three plots (-90 to 90).
When I plot a hemisphere without limiting its extent, however, I get a round bounding box, as expected from an orthographic projection:
Using plt.xlim(-90,-50) results in a vertical stripe, and plt.ylim(-90,-50) in a horizontal stripe, so that is no solution either.
How can I limit the latitudinal extent of my orthographic projection, whilst maintaining the circular bounding box?
The code to produce above graphs:
import numpy as np
from matplotlib import pyplot as plt
import cartopy.crs as ccrs
# Create dummy data, latitude from -90(S) to 90 (N), lon from -180 to 180
theta, phi = np.meshgrid(np.arange(0,180),np.arange(0,360));
theta = -1*(theta.ravel()-90)
phi = phi.ravel()-180
radii = theta
# Make masks for hemispheres and central
mask_central = np.abs(theta) < 60
mask_north = theta > 40
mask_south = theta < -40
data_crs= ccrs.PlateCarree() # Data CRS
# Grab map projections for various plots
map_proj = ccrs.Mollweide(central_longitude=0)
map_proj_N = ccrs.Orthographic(central_longitude=0, central_latitude=90)
map_proj_S = ccrs.Orthographic(central_longitude=0, central_latitude=-90)
fig = plt.figure()
ax1 = fig.add_subplot(2, 1, 2,projection=map_proj)
im1 = ax1.scatter(phi[mask_central],
theta[mask_central],
c = radii[mask_central],
transform=data_crs,
vmin = -90,
vmax = 90,
)
ax1.set_title('Central latitudes')
ax_N = fig.add_subplot(2, 2, 1, projection=map_proj_N)
ax_N.scatter(phi[mask_north],
theta[mask_north],
c = radii[mask_north],
transform=data_crs,
vmin = -90,
vmax = 90,
)
ax_N.set_title('Northern hemisphere')
ax_S = fig.add_subplot(2, 2, 2, projection=map_proj_S)
ax_S.scatter(phi[mask_south],
theta[mask_south],
c = radii[mask_south],
transform=data_crs,
vmin = -90,
vmax = 90,
)
ax_S.set_title('Southern hemisphere')
fig = plt.figure()
ax = fig.add_subplot(111,projection = map_proj_N)
ax.scatter(phi,
theta,
c = radii,
transform=data_crs,
vmin = -90,
vmax = 90,
)
ax.set_title('Northern hemisphere')
plt.show()
(1). In all of your plots, when you use scatter(), the size of the scatter points should be defined with proper s=value, otherwise the default value is used. I use s=0.2 and the resulting plots look better.
(2). For 'Central latitudes' case, you need to specify correct y-limits with set_ylim(). This involves the computation of them. The use of transform_point() is applied here.
(3). For the remaining plots that require elimination of unneeded features, proper circular clip paths can be used. Their perimeters are also used to plot as map boundaries in both cases. Their existence may cause trouble plotting other map features (such as coastlines) as I demonstrate with the code and its output.
# original is modified and extended
import numpy as np
from matplotlib import pyplot as plt
import cartopy.crs as ccrs
import matplotlib.patches as mpatches # need it to create clip-path
# Create dummy data, latitude from -90(S) to 90 (N), lon from -180 to 180
theta, phi = np.meshgrid(np.arange(0,180),np.arange(0,360));
theta = -1*(theta.ravel()-90)
phi = phi.ravel()-180
radii = theta
# Make masks for hemispheres and central
mask_central = np.abs(theta) < 60
mask_north = theta > 40
mask_south = theta < -40
data_crs= ccrs.PlateCarree() # Data CRS
# Grab map projections for various plots
map_proj = ccrs.Mollweide(central_longitude=0)
map_proj_N = ccrs.Orthographic(central_longitude=0, central_latitude=90)
map_proj_S = ccrs.Orthographic(central_longitude=0, central_latitude=-90)
# 'Central latitudes' plot
fig = plt.figure()
ax1 = fig.add_subplot(2, 1, 2, projection=map_proj)
# Note: Limits of plot depends on plotting data, but not exact!
im1 = ax1.scatter(phi[mask_central],
theta[mask_central],
s = 0.2,
c = radii[mask_central],
transform=data_crs,
vmin = -90,
vmax = 90,
)
# compute y limits
_, y_btm = map_proj.transform_point(0, -60, ccrs.Geodetic())
_, y_top = map_proj.transform_point(0, 60, ccrs.Geodetic())
# apply y limits
ax1.set_ylim(y_btm, y_top)
ax1.coastlines(color='k', lw=0.35)
ax1.set_title('Central latitudes')
ax_N = fig.add_subplot(2, 2, 1, projection=map_proj_N)
ax_N.scatter(phi[mask_north],
theta[mask_north],
s = 0.1, # not mandatory
c = radii[mask_north],
transform=data_crs,
vmin = -90,
vmax = 90,
)
# use a circular path as map boundary
clip_circle = mpatches.Circle(xy=[0,0], radius=4950000, facecolor='none', edgecolor='k')
ax_N.add_patch(clip_circle)
ax_N.set_boundary(clip_circle.get_path(), transform=None, use_as_clip_path=True)
# with `use_as_clip_path=True` the coastlines do not appear
ax_N.coastlines(color='k', lw=0.75, zorder=13) # not plotted!
ax_N.set_title('Northern hemisphere1')
# 'Southern hemisphere' plot
ax_S = fig.add_subplot(2, 2, 2, projection=map_proj_S)
ax_S.scatter(phi[mask_south],
theta[mask_south],
s = 0.02,
c = radii[mask_south],
transform=data_crs,
vmin = -90,
vmax = 90,
)
clip_circle = mpatches.Circle(xy=[0,0], radius=4950000, facecolor='none', edgecolor='k')
ax_S.add_patch(clip_circle)
# applying the clip-circle as boundary, but not use as clip-path
ax_S.set_boundary(clip_circle.get_path(), transform=None, use_as_clip_path=False)
# with `use_as_clip_path=False` the coastlines is plotted, but goes beyond clip-path
ax_S.coastlines(color='k', lw=0.75, zorder=13)
ax_S.set_title('Southern hemisphere')
# 'Northern hemisphere2' plot, has nice circular limit
fig = plt.figure()
ax = fig.add_subplot(111,projection = map_proj_N)
ax.scatter(phi,
theta,
s = 0.2,
c = radii,
transform=data_crs,
vmin = -90,
vmax = 90,
)
ax.coastlines(color='k', lw=0.5, zorder=13)
ax.set_title('Northern hemisphere2')
ax.set_global()
plt.show()
The output plot:
