polar pcolormesh plot projected onto cartopy map-CodePudding

To simplify, as much as possible, a

Over a cartopy map like this...

flatMap = ccrs.PlateCarree()
resolution = '110m'
fig = plt.figure(figsize=(12,6), dpi=96)
ax = fig.add_subplot(111, projection=flatMap)

ax.imshow(np.tile(np.array([[cfeature.COLORS['water'] * 255]], dtype=np.uint8), [2, 2, 1]), origin='upper', transform=ccrs.PlateCarree(), extent=[-180, 180, -180, 180])
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', resolution, edgecolor='black', facecolor=cfeature.COLORS['land']))
ax.pcolormesh(X,Y,Z, cmap=mpl.cm.jet_r,shading='auto') 
gc.collect()

I'd like to project this polar plot over an arbitrary lon/lat... I can convert the polar theta/phi into lon/lat, but lon/lat coords (used on the map) are more 'cartesian like' than polar, hence you cannot just substitute lon/lat for theta/phi ... This is a conceptual problem. How would you tackle it?

CodePudding user response：

Firstly, the data must be prepared/transformed into certain projection coordinates for use as input. And the instruction/option of the data's CRS must be specified correctly when used in the plot statement.

In your specific case, you need to transform your data into (long,lat) values.

XX = X/np.pi*180  # wrap around data in EW direction
YY = Y*9          # spread across N hemisphere

And plot it with an instruction transform=ccrs.PlateCarree().

ax.pcolormesh(XX,YY,Z, cmap=mpl.cm.jet_r,shading='auto', 
          transform=ccrs.PlateCarree())

The same (XX,YY,Z) data set can be plotted on orthographic projection.

CodePudding user response：

This solution does NOT account for the projection point being at some altitude above the globe... I can do that part, so I really have trouble mapping the meshgrid to lon/lat so the work with the PREVIOUSLY GENERATES values of Z.

Here's a simple mapping directly from polar to cart:

X_cart = np.array([[p*np.sin(t) for p in phis] for t in thetas]).T
Y_cart = np.array([[p*np.cos(t) for p in phis] for t in thetas]).T
# Need to map cartesian XY to Z that is compatbile with above... 

Z_cart = np.sin(X)**10   np.cos(10   Y*X) * np.cos(X) # This Z does NOT map to cartesian X,Y
print(X_cart.shape,Y_cart.shape,Z_cart.shape) 

flatMap = ccrs.PlateCarree()
resolution = '110m'
fig = plt.figure(figsize=(12,6), dpi=96)
ax = fig.add_subplot(111, projection=flatMap)

ax.imshow(np.tile(np.array([[cfeature.COLORS['water'] * 255]], dtype=np.uint8), [2, 2, 1]), origin='upper', transform=ccrs.PlateCarree(), extent=[-180, 180, -180, 180])
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', resolution, edgecolor='black', facecolor=cfeature.COLORS['land']))

im = ax.pcolormesh(X_cart*2,Y_cart*2, Z_cart, cmap=mpl.cm.jet_r, shading='auto') # c=mapper.to_rgba(Z_cart), cmap=mpl.cm.jet_r) 


gc.collect()

Which maps the polar plot center to lon/lat (0,0):

I'm close... I somehow need to move my cartesian coords to the proper lon/lat (the satellite sub-point) and then scale it appropriately. Have the set of lon/lat but I'm screwing up the meshgrid somehow... ???

The sphere_intersect() routine returns lon/lat for projection of theta/phi on the globe (that works)... The bit that doesn't work is building the meshgrid that replaces X,Y:

lons = np.array([orbits.sphere_intersect(SV_pos_vec, SV_vel_vec, az << u.deg, el << u.deg, 
                                        lonlat=True)[0] for az in thetas for el in phis], dtype='object')
lats = np.array([orbits.sphere_intersect(SV_pos_vec, SV_vel_vec, az << u.deg, el << u.deg, 
                                        lonlat=True)[1] for az in thetas for el in phis], dtype='object')

long, latg = np.meshgrid(lons,lats) # THIS IS A PROBLEM I BELIEVE...

and the pcolormesh makes a mess...