GEOSChem-python-tutorial
GEOSChem-python-tutorial copied to clipboard
3D Globe Map showing a specie's concentration
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.