I want to plot a velocity field (u, v) as a function of (longitude, latitude) as a streamplot, in plate-carree (lon-lat) and orthographic projections, using the matplotlib function plt.streamplot with the Cartopy package. I want the starting points for the stream plot fixed, as I would like to eventually create an animation of some features rotating over the sphere with the streamlines "anchored" somewhere. However I experience a number of issues with these projections, detailed further below.
I use the following package versions: Python 3.7.6, matplotlib 3.2.2, cartopy 0.18.0
A "minimal" example that creates streamplots with and without starting points for the mentioned projections would be:
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
# define longitude and latitude grids
nlon, nlat = 48, 24
dlon, dlat = 7.5, 7.5
lons = np.linspace(-180+dlon,180,nlon)
lats = np.linspace(-90,90-dlat,nlat)
# define random velocity fields
np.random.seed(999)
u = np.random.rand(nlat, nlon)
v = np.random.rand(nlat, nlon)
# get indices of latitudes north of the equator and longitudes within +/- 90 deg
ilat_north = np.where(lats >= 0)[0]
ilon = np.where(np.abs(lons) <= 90)[0]
# create a list of 5 starting points at various latitudes aligned at a longitude of 0 deg
start_points = list(zip(list(np.linspace(0,0,5)),list(np.linspace(0,80,5))))
# output plot directory and whether plots should be saved
plot_dir = '/home/proxauf/test_streamplot_cartopy'
save_plots = True
def plot_streamplot(mode, start_points=None, figtitle=None, outname=None, save_plot=False):
fig = plt.figure()
if mode == 'mpl-platecarree':
# a Plate-Carree projection without Cartopy (directly using matplotlib)
ax = plt.gca()
ax.streamplot(lons, lats[ilat_north], u[ilat_north,:], v[ilat_north,:], start_points=start_points)
elif 'cartopy' in mode:
if mode == 'cartopy-platecarree':
# a Plate-Carree projection with Cartopy
ax = plt.subplot(projection=ccrs.PlateCarree())
ax.streamplot(lons, lats[ilat_north], u[ilat_north,:], v[ilat_north,:], start_points=start_points, transform=ccrs.PlateCarree())
elif mode == 'cartopy-orthographic0':
# an orthographic projection with a latitudinal inclination of 0 deg with Cartopy
# restrict longitudes to +/- 90 deg, otherwise an error occurs (outside of projection area)
ax = plt.subplot(projection=ccrs.Orthographic(central_longitude=0,central_latitude=0))
ax.streamplot(lons[ilon], lats[ilat_north], u[ilat_north[:,None],ilon[None,:]], v[ilat_north[:,None], ilon[None,:]], start_points=start_points, transform=ccrs.PlateCarree())
if mode == 'cartopy-orthographic90':
# an orthographic projection with a latitudinal inclination of 90 deg with Cartopy
ax = plt.subplot(projection=ccrs.Orthographic(central_longitude=0,central_latitude=90))
ax.streamplot(lons, lats[ilat_north], u[ilat_north,:], v[ilat_north,:], start_points=start_points, transform=ccrs.PlateCarree())
if 'platecarree' in mode:
ax.set_xlim(-180,180)
ax.set_ylim(-90,90)
ax.set_aspect('equal')
if 'cartopy' in mode:
# draw gridlines with label for visual orientation
ax.gridlines(draw_labels=True)
if start_points is not None:
# plot starting points for the streamplot
for i in start_points:
if 'cartopy' in mode:
ax.plot(i[0],i[1],marker='o',transform=ccrs.PlateCarree())
else:
ax.plot(i[0],i[1],marker='o')
fig.suptitle(figtitle)
if save_plot:
# save the plot and close it
plt.savefig('%s/%s' % (plot_dir, outname))
plt.close()
# create and save streamplots for different projections, with and without starting points
plot_streamplot(mode='mpl-platecarree', start_points=None, figtitle='mpl-platecarree, without start points', outname='test_streamplot_mpl_pc_sp_no.png', save_plot=save_plots)
plot_streamplot(mode='cartopy-platecarree', start_points=None, figtitle='cartopy-platecarree, without start points', outname='test_streamplot_cartopy_pc_sp_no.png', save_plot=save_plots)
plot_streamplot(mode='cartopy-orthographic0', start_points=None, figtitle='cartopy-orthographic0, without start points', outname='test_streamplot_cartopy_ortho0_sp_no.png', save_plot=save_plots)
plot_streamplot(mode='cartopy-orthographic90', start_points=None, figtitle='cartopy-orthographic90, without start points', outname='test_streamplot_cartopy_ortho90_sp_no.png', save_plot=save_plots)
plot_streamplot(mode='mpl-platecarree', start_points=start_points, figtitle='mpl-platecarree, with start points', outname='test_streamplot_mpl_pc_sp_yes.png', save_plot=save_plots)
plot_streamplot(mode='cartopy-platecarree', start_points=start_points, figtitle='cartopy-platecarree, with start points', outname='test_streamplot_cartopy_pc_sp_yes.png', save_plot=save_plots)
plot_streamplot(mode='cartopy-orthographic0', start_points=start_points, figtitle='cartopy-orthographic0, with start points', outname='test_streamplot_cartopy_ortho0_sp_yes.png', save_plot=save_plots)
plot_streamplot(mode='cartopy-orthographic90', start_points=start_points, figtitle='cartopy-orthographic90, with start points', outname='test_streamplot_cartopy_ortho90_sp_yes.png', save_plot=save_plots)
Below is the output created by this minimal example.
Case 1: Mpl-Plate-Carree without starting points: works fine:
Case 2: Cartopy-Plate-Carree without starting points: works fine:
Case 3: Cartopy-Orthographic-0 without starting points: works fine (?), but high latitudes are cut off:
Case 4: Cartopy-Orthographic-90 without starting points: works fine:
Case 5: Mpl-Plate-Carree with starting points: works fine:
Case 6: Cartopy-Plate-Carree with starting points: two streamlines for two starting point are not shown:
Case 7: Cartopy-Orthographic-0 with starting points: no data are shown (no streamlines, starting points, etc.):
Case 8: Cartopy-Orthographic-90 with starting points: only a part of the northern hemisphere is shown and the streamline seems incorrect:
Summary: Cases 6, 7 and 8 (Cartopy Plate-Carree/Orthographic-0deg-inclination/Orthographic-90deg-inclination projections) give an unexpected, probably incorrect output. Can somebody explain to me why the plots look that way? Am I misinterpreting the input that I should give?
Small bonus question: If I want to set axis limits for the orthographic projections (e.g. latitude limits of [-90,+90] deg for 0 deg inclination, to [0, +90] deg for 90 deg inclination, and also longitude limits), e.g. assuming I only plot data above 30 deg, how does this work? ax.set_xlim or ax.set_ylim don't work, I assume (since they do not know the projection). If this needs to be done via ax.set_extent, which longitude and latitude limits would I need to set there?
