-
Notifications
You must be signed in to change notification settings - Fork 15
Visualizing GEOS Cubed‐Sphere History Output
This page details ways to visualize GEOS Cubed-Sphere (CS) History netCDF output. Note that this does not include GEOS CS restarts as those use a slightly different file format. If needed, the GMAO SI Team has scripts which allow a user to convert from restart to history format.
The Panoply visualization software has native support for GEOS Cubed-Sphere output. Simply open the file and visualize the data as you would for a lat-lon output.
Below is an image of T2M from a GEOS CS C360 output file. The data is visualized using the viridis colormap.

For those desiring a Python workflow to visualize GEOS CS output, the uxarray package is the best method. As of version v2024.06.0, the package has the ability to read GEOS CS output directly. Install as usual with either pip or conda.
We highly recommend installing the latest version of UXarray as they have very good documentation at:
https://uxarray.readthedocs.io/
and as versions change, so can some of the code. All code below was run with version v2024.12.0
import uxarray as ux
import geoviews as gv
import geoviews.feature as gf
import cartopy
import cartopy.feature as cf
import cartopy.crs as ccrs
graticules = cf.NaturalEarthFeature(
category='physical',
name='graticules_15',
scale='50m')
features = gf.coastline(
projection=ccrs.PlateCarree(), line_width=1, scale="50m"
) * gf.borders(
projection=ccrs.PlateCarree(), line_width=1, scale="50m"
) * gv.Feature(graticules, group='Lines')
file_path = r"/path/to/file.nc4"
uxds = ux.open_dataset(file_path, file_path)
uxds['var'][0].plot.polygons(
rasterize=True,
cmap='viridis',
backend='bokeh',
height=500, width=1000) * featuresUsing the above code with a GEOS CS C360 T2M output (same as used for Panoply above) produces the following visualization.

NOTE: The above code uses the bokeh backend and the image was produced by saving from it. You can also use backend='matplotlib' to produce a matplotlib image.
As you can see from the above image, UXarray by default does not plot data that crosses the antimeridian. Thus, you get artifacts on the edges.
UXarray has different methods for handling this. The default is perodic_elements='exclude' which is the fastest method. If you desire these elements plotted, you can set periodic_elements='split'. For example:
uxds['var'][0].plot.polygons(
rasterize=True,
cmap='viridis',
backend='bokeh',
height=500, width=1000,
periodic_elements='split') * featureswhich leads to the following image:

Note that using split has a performance penalty. In testing the above code with %%time, the original image took:
CPU times: user 1.5 s, sys: 53.7 ms, total: 1.55 s
Wall time: 1.55 s
while with split:
CPU times: user 20.6 s, sys: 444 ms, total: 21 s
Wall time: 21.2 s
UXarray estimates a roughly 10-20x performance penalty for using split so use with caution.