cartopy rotate winds incorrectly?

Viewed 385

I am trying to plot some wind vectors on maps using cartopy. However, the winds look strange so I did a simple tests using wind vectors at a ring of points at 75N. The wind vectors should points 45 degree away from both lat/lon grids because u-winds and v-winds are of the same magnitude under the stereographic projection, which should preserve the angle.

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

pcproj = ccrs.PlateCarree()
lon0 = -150
mapproj = ccrs.Stereographic(
    central_longitude=lon0,central_latitude=75,
    true_scale_latitude=75,
    )
XLIM = 300e3; YLIM=300e3
dm =5; dp=1
fig = plt.figure(0,(7,7))
ax  = fig.add_axes([0.1,0.1,0.85,0.9],projection=mapproj)
      
ax.set_extent([-XLIM,XLIM,-YLIM,YLIM],crs=mapproj)
ax.coastlines(resolution='50m',color='.5',linewidth=1.5)
lon_grid = np.arange(-180,181,dm)
lat_grid = np.arange(-80,86,dp)
gl = ax.gridlines(draw_labels=True,
                  xlocs=lon_grid,ylocs=lat_grid,
                  x_inline=False,y_inline=False,
                  color='k',linestyle='dotted')

# --- draw 45 degree winds at 75N
lon = np.linspace(0,360,73)
lat = np.ones(len(lon))*75
uu  = np.ones(len(lon))*10
vv  = uu*1.
pts = mapproj.transform_points(pcproj,lon,lat)
xx = pts[...,0]; yy = pts[...,1]
ux,vx = mapproj.transform_vectors(pcproj,lon,lat,uu,vv)

ax.quiver(xx,yy,ux,vx,transform=mapproj)

It is clear that the wind vectors does not point 45 degrees from lat/lon grids.

winds using transform_vectors from cartopy

I check the ux, vx at 150W and 75N, the center of the projection. The wind here should be the same as in lat/lon grid but it is (3.54,13.7) rather than (10,10). The wind vectors looks the same as using this line

ax.quiver(lon,lat,uu,vv,transform=pcproj)

This may not be a surprise because I think quiver uses transform_vector and transform_points under the hood. I tried winds along lat/lon grid directions, they are transformed correctly. Is this a bug in transform_vectors or am I using it incorrectly?

UPDATE: As suggested by @swatchai, right now, u-wind should be divided by cos(lat) before being fed in to transform_vectors. I doubt this is the intended behavior but this normalization should be used until a future cartopy update changed the behavior. They may choose just to update the docstring. Still, something to be caution about in future cartopy updates.

1 Answers

When (lat, long) in degrees are used, a small angular distance ds on spherical surface has 2 components:

`R*dlat`  in NS direction
`R*dlon*cos(lat)`  in EW direction

where R is the earth radius,

the second component varies with latitude of the location considered,

with this in mind, the computation of ux, vx at latitude 75 degrees should be

rho = np.pi/180.
ux,vx = mapproj.transform_vectors(pcproj,lon,lat, uu/np.cos(75*rho), vv)
Related