GEOSChem-python-tutorial icon indicating copy to clipboard operation
GEOSChem-python-tutorial copied to clipboard

3D Globe Map showing a specie's concentration

Open nicodgomez opened this issue 1 year ago • 0 comments

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.

nicodgomez avatar Jun 24 '24 18:06 nicodgomez