Contour plot and Cartopy - Geostationary satellite view [Plots, Code,Attempt]

Question

I have data which is in X,Y co-ordinates that is transformed from spherical co-ordinates and I am trying to plot a contour of this data on top of a globe which I have taken from Cartopy. But its not working and I have been trying for a few hours googling and attempts at this! Any suggestions would be greatly appreciated.

I've amended the code to include dummy data, this is the data to be contour and overlay on the globe.

There are sample plots in below links too

import numpy as np
import matplotlib.pyplot as plt
# import cartopy.feature as cf
import cartopy.crs as ccrs
import cartopy.feature as cf

plt.close('all')


#%% Co-ordindate System
resolution = 201
theta = np.linspace(-90, 90, num=resolution) * np.pi/180
phi = np.linspace(0, 360, num=resolution) * np.pi/180
theta_grid, phi_grid = np.meshgrid(theta,phi)

u_grid = np.sin(theta_grid) * np.cos(phi_grid)
v_grid = np.sin(theta_grid) * np.sin(phi_grid)


#%% Dummy Data - End result called 'AF' is what is important, please ignore how its done. I cant get 'AF' to plot ontop of a globe.
x_num = 10
y_num = 10
num_elements = x_num * y_num
dx = 2
dy = 2
[array_x, array_y] = np.mgrid[0:x_num, 0:y_num]
array_x = (array_x - (x_num- 1) / 2)
array_y = np.flipud((array_y - (y_num - 1) / 2))
array_x = np.transpose(array_x) * dx
array_y = np.flip(np.transpose(array_y)) * dy  
array_x = np.reshape(array_x, (x_num * y_num, -1)).flatten()
array_y = np.reshape(array_y, (x_num * y_num, -1)).flatten()

AF = np.zeros(np.size(u_grid),dtype = np.complex128) #Make sure u_grid and v_grid are the same size
for n in range(num_elements):

    kx = array_x[n] * u_grid.ravel()
    ky = array_y[n] * v_grid.ravel()
    kr = kx+ky
    AF = AF  + np.exp(1j*kr)
    
AF = np.resize(AF, np.shape(u_grid))
AF = AF / np.max(abs(AF))
AF = 20*np.log10(abs(AF))
AF[AF< -50] =-50  


#%% Plots

c = np.linspace(np.min(AF),np.max(AF),50)
)

fig = plt.figure('UV') 
ax = fig.add_subplot(111)
plt.contourf(u_grid,v_grid,AF,c,extend="both", cmap = 'Spectral_r')  #plasma #gnuplot2
plt.colorbar(label ='Directivity (dBi)')
ax.set_xlabel('u = sin(theta)cos(phi)')
ax.set_ylabel('v = sin(theta)sin(phi)')
ax.legend(loc=1)


projection=ccrs.Geostationary(central_longitude=0.0, satellite_height=35785831)

fig = plt.figure('No Overlay') 
ax = plt.axes(projection=projection)
ax.add_feature(cf.COASTLINE)
plt.show()


fig = plt.figure('Attempt to Overlay') 
ax = plt.axes(projection=projection)
plt.contourf(u_grid,v_grid,AF,c,extend="both", cmap = 'Spectral_r', alpha = 0.2)  #plasma #gnuplot2
ax.add_feature(cf.COASTLINE)
plt.show()

Example plots of the data and the failing to overlay countour onto globe

Data Contour Plot in XY transform of spherical coordinates:

The Globe I want to plot a contour on top of:

I have tried many times googling the problem, Cartopy website, posts on stackexchange but I have not got this working yet... any suggestions greatly appreciated.

Answer 1

Since you want to show your contour plot together with a geographic map, you have to use appropriate coordinates, i.e. longitude and latitude in degrees, like in many of the Cartopy examples, e.g. in this one, where they convert the original values in radians into degrees. In your code it could be something like theta_deg = np.rad2deg(theta) and phi_deg = np.rad2deg(phi). Also, you have to set the transform parameter in matplotlib.pyplot.contourf(). Here is a simplified example using your data:

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

#co-ordindate system
resolution = 201
lats = np.linspace(-90, 90, num=resolution)
lons = np.linspace(0, 360, num=resolution)
lons, lats = np.meshgrid(lons, lats)

# load numpy array of data
AF = np.load('Plotted_data.npy')

projection=ccrs.Geostationary(central_longitude=0.0, satellite_height=35785831)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=projection)
# show coastlines
ax.coastlines()
# add color filled contours
filled_c = ax.contourf(lons, lats, AF, transform=ccrs.PlateCarree(), extend="both", cmap='Spectral_r', alpha=0.2)

plt.show()

The result:

Probably it doesn't look exactly as you want, esp. that ugly abrupt transition which happens to be right at the prime meridian. But this is how the data from Plotted_data.npy really maps onto the globe when using the geographical coordinates. You may need to verify the data or perform some kind of data transformation.