polar pcolormesh plot projected onto cartopy map

Question

To simplify, as much as possible, a question I already asked, how would you OVERLAY or PROJECT a polar plot onto a cartopy map.

phis = np.linspace(1e-5,10,10) # SV half cone ang, measured up from nadir
thetas = np.linspace(0,2*np.pi,361)# SV azimuth, 0 coincides with the vel vector

X,Y = np.meshgrid(thetas,phis)
Z   = np.sin(X)**10 + np.cos(10 + Y*X) * np.cos(X)

fig, ax = plt.subplots(figsize=(4,4),subplot_kw=dict(projection='polar'))
im = ax.pcolormesh(X,Y,Z, cmap=mpl.cm.jet_r,shading='auto') 

ax.set_theta_direction(-1)
ax.set_theta_offset(np.pi / 2.0)
ax.grid(True)

that results in

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?

Answer 1

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.

Edit1

Update of the code and plots.

Part 1 (Data)

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

#
# Part 1
#
phis = np.linspace(1e-5,10,10) # SV half cone ang, measured up from nadir
thetas = np.linspace(0,2*np.pi,361)# SV azimuth, 0 coincides with the vel vector

X,Y = np.meshgrid(thetas,phis)
Z   = np.sin(X)**10 + np.cos(10 + Y*X) * np.cos(X)

fig, ax = plt.subplots(figsize=(4,4),subplot_kw=dict(projection='polar'))
im = ax.pcolormesh(X,Y,Z, cmap=mpl.cm.jet_r,shading='auto') 

ax.set_theta_direction(-1)
ax.set_theta_offset(np.pi / 2.0)
ax.grid(True)

Part 2 The required code and output.

#
# Part 2
#
flatMap = ccrs.PlateCarree()
resolution = '110m'
fig = plt.figure(figsize=(12,6), dpi=96)

ax = fig.add_subplot(111, projection=flatMap)

ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', 
            resolution, edgecolor='black', alpha=0.7, 
            facecolor=cfeature.COLORS['land']))
ax.set_extent([-180, 180, -90, 90], crs=ccrs.PlateCarree())


def scale_position(lat_deg, lon_deg, rad_deg):
    # Two operations:
    # 1. manipulates X,Y data and get (XX,YY)
    # 2. create proper projection of (XX,YY), `rotpole_proj`
    # Returns: XX,YY,rotpole_proj

    # For X data
    XX = X/np.pi*180   #always wrap around EW direction

    # For Y data
    # The cone data: min=0, max=10  -->  (90-rad),90
    # rad_deg = radius of the display area
    top = 90
    btm = top-rad_deg
    YY = btm + (Y/Y.max())*rad_deg

    # The proper coordinate system
    rotpole_proj = ccrs.RotatedPole(pole_latitude=lat_deg, pole_longitude=lon_deg)
    # Finally,
    return XX,YY,rotpole_proj

# Location 1 (Asia)
XX1, YY1, rotpole_proj1 = scale_position(20, 100, 20)
ax.pcolormesh(XX1, YY1, Z, cmap=mpl.cm.jet_r, 
              transform=rotpole_proj1) 

# Location 2 (Europe)
XX2, YY2, rotpole_proj2 = scale_position(62, -6, 8)
ax.pcolormesh(XX2, YY2, Z, cmap=mpl.cm.jet_r,
              transform=rotpole_proj2) 

# Location 3 (N America)
XX3, YY3, rotpole_proj3 = scale_position(29, -75, 30)
ax.pcolormesh(XX3, YY3, Z, cmap=mpl.cm.jet_r,shading='auto', 
              transform=rotpole_proj3) 

#gc.collect()
plt.show()