3D Globe Map showing a specie's concentration

[nicodgomez]

Hello GeosChem community. I have found this tutorial very useful to analyze GeosChem outputs. Also I have learned a lot about python programming and data types that can be treated. Recently I have been trying to illustrate species concentration in a 3D Interactive Globe Map. I have been trying the following code without success:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import xarray as xr # the major tool to work with NetCDF data!
import plotly.graph_objects as go
import plotly.express as px

ds = xr.open_dataset("initial_GEOSChem_rst.4x5_tropchem.nc")
dr = ds['TRC_O3']
dr_surf = dr.isel(time=0, lev=0)

latitudes = dr_surf['lat'].values
longitudes = dr_surf['lon'].values
concentracion = dr_surf.values

fig = go.Figure()
fig.add_trace(go.Contour(
z=concentracion,
x=longitudes,
y=latitudes,
colorscale='Viridis',
colorbar=dict(title='Concentración de Ozono')
))

fig.update_geos(
projection_type="orthographic",
showcountries=True, countrycolor="Black",
showcoastlines=True, coastlinecolor="Black"
)

fig.update_layout(
title="Concentración de Ozono en un Mapa Esférico Interactivo",
geo=dict(
showland=True,
landcolor="white",
showocean=True,
oceancolor="LightBlue"
)
)

The code runs without error, but it only draws a contour map without coastlines, or the orthographic projection. I have used the same test archive as the tutorial, "initial_GEOSChem_rst.4x5_tropchem.nc".

I was able to reproduce a 3D Globe Map, in static format, using the following code:

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

from gamap_colormap import WhGrYlRd

archivo_nc = 'initial_GEOSChem_rst.4x5_tropchem.nc'
datos = xr.open_dataset(archivo_nc)

ozono = datos['TRC_O3']

ozono_surf = ozono.isel(time=0, lev=46)

latitudes = ozono_surf['lat'].values
longitudes = ozono_surf['lon'].values
concentracion_ozono = ozono_surf.values

fig = plt.figure(figsize=(10, 5))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=50, central_latitude=40))

ax.coastlines()
ax.gridlines(linestyle='--')

lon_grid, lat_grid = np.meshgrid(longitudes, latitudes)
mesh = ax.pcolormesh(lon_grid, lat_grid, concentracion_ozono, transform=ccrs.PlateCarree(), cmap=WhGrYlRd)

cbar = plt.colorbar(mesh, orientation='horizontal', pad=0.05)
cbar.set_label('Concentración de Ozono')

plt.title('Concentración de Ozono en un Mapa Esférico')

plt.show()

I would like to do an equivalent thing, but with the capability of rotating the globe.

Thank you in advance.