SciTools / cartopy

Cartopy - a cartographic python library with matplotlib support
https://scitools.org.uk/cartopy/docs/latest
BSD 3-Clause "New" or "Revised" License
1.43k stars 365 forks source link

'regrid_shape' produces incorrect wind barbs near the pole #1634

Open tomerburg opened 4 years ago

tomerburg commented 4 years ago

I'm generating numerical model output plots (e.g., GFS, ERA-5) for polar projections using the Nearside Perspective projection. When plotting wind barbs, the barb magnitudes are correct if plotted normally, such as:

ax.barbs(x, y, u, v, transform=ccrs.PlateCarree())

Once I add the regrid_shape argument, the barbs are inaccurate near the North Pole, with the magnitude decreasing towards zero within a few degrees latitude of the pole. I'm not entirely sure if this is user error or a bug in regrid_shape, but I am struggling to find any error on my end so I suspect it may be a bug.

I created a sample code & data to reproduce this, by assuming a constant global u=50 m/s and v=0 m/s, meaning the wind magnitude should be 50 m/s throughout the grid. Using regrid_shape, most of the barbs correctly show a 50 m/s flag, but right near the pole the barbs gradually decrease towards 0 m/s. If regrid_shape is omitted from ax.barbs(), then all barbs (including those near the pole) show a 50 m/s flag.

This is produced using Cartopy 0.17.0, in multiple versions of python (anywhere from 3.6.7 to 3.8.0) and on both Windows & Linux operating systems. Any help to identify if this is a user error or a bug, and if the latter then if it can be fixed, would be greatly appreciated.

import numpy as np
import matplotlib.pyplot as plt
from cartopy import crs as ccrs
import cartopy.feature as cfeature

#Retrieve instance of NearsidePerspective projection
proj = ccrs.NearsidePerspective(
    central_longitude = -100.0,
    central_latitude = 90.0,
    satellite_height = 4785831,
)

#Create figure
fig = plt.figure(figsize=(14,9))
ax = fig.add_axes([0.05,0.03,0.89,0.90],projection=proj)

#Plot geography boundaries
countries = ax.add_feature(cfeature.BORDERS.with_scale('50m'),linewidths=1.0,linestyle='solid')
coastlines = ax.add_feature(cfeature.COASTLINE.with_scale('50m'),linewidths=1.0,linestyle='solid')
continent_mask = ax.add_feature(cfeature.LAND.with_scale('50m'),facecolor='#eeeeee',edgecolor='face')

#Add sample data
lon = np.arange(-180,180,4) #lat/lon grid every 4 degrees
lat = np.arange(35,90,4) #lat/lon grid every 4 degrees
lons,lats = np.meshgrid(lon,lat)

u = np.zeros((lons.shape)) + 50.0 #set u-wind to 50 m/s throughout the grid
v = np.zeros((lons.shape)) + 0.0 #set v-wind to 0 m/s throughout the grid

#Plot barbs
ax.barbs(lons,lats,u,v,regrid_shape=(40,30),transform=ccrs.PlateCarree(),linewidth=0.5,length=6)

#Show plot and close
plt.show()
plt.close()
greglucas commented 4 years ago

Regridding is done in the projection frame and on either side of the pole you have a vector going opposite directions, so the vector will be interpolated to 0 at the pole, hence your decreasing values.

I think the intent you probably expected is for the interpolation to be done in the source frame, and I think this is what I would expect as well. I've put up a draft PR https://github.com/SciTools/cartopy/pull/1636 with that thought in mind. You should be able to adapt that methodology of transforming back and forth to do the interpolation yourself in the meantime.

dopplershift commented 4 years ago

I'll also add that this seems to depend on what you choose for your latitude range on input. I adjusted to start your lats at 37, thus giving me the highest lat point at 89 rather than 87, and the artifact went away.