Conversion to GeoDataFrame for Visualization with HoloViz Packages#

Authors: Philip Chmielowiec, Ian Franda

Overview#

This usage example notebook showcases how UXarray data structures (Grid, UxDataArray) can be converted in a Spatialpandas GeoDataFrame, which can be used directly by packages from the HoloViz stack, such as hvPlot, Datashader, Holoviews, and Geoviews. In addition to showing how to convert to a GeoDataFrame, a series of visualizations using hvPlot and GeoViews is showcased based around the converted data structure.

Note

While the visualizations in this notebook are not rendered using UXarray functionality, development of UXarray’s own visualization functions can be tracked here.

Imports#

This notebook requires the following visualization packages to be installed in order to run locally.

conda install -c conda-forge holoviews hvplot geoviews spatialpandas antimeridian
import uxarray as ux
import numpy as np
import warnings

import matplotlib
import matplotlib.pyplot as plt

import hvplot.pandas
import holoviews as hv

warnings.filterwarnings("ignore")

Data#

For this notebook, we will be using E3SM model output, which uses a cubed-sphere grid. However, this notebook can be adapted to datasets in any of our supported formats.

base_path = "../../test/meshfiles/ugrid/outCSne30/"
grid_path = base_path + "outCSne30.ug"
data_path = base_path + "outCSne30_vortex.nc"
uxds = ux.open_dataset(grid_path, data_path)
uxds
<xarray.UxDataset> Size: 43kB
Dimensions:  (n_face: 5400)
Dimensions without coordinates: n_face
Data variables:
    psi      (n_face) float64 43kB ...
uxds.uxgrid
<uxarray.Grid>
Original Grid Type: UGRID
Grid Dimensions:
  * n_node: 5402
  * n_face: 5400
  * n_max_face_nodes: 4
  * n_nodes_per_face: (5400,)
Grid Coordinates (Spherical):
  * node_lon: (5402,)
  * node_lat: (5402,)
Grid Coordinates (Cartesian):
Grid Connectivity Variables:
  * face_node_connectivity: (5400, 4)
Grid Descriptor Variables:

Conversion to spatialpandas.GeoDataFrame for Visualization#

In order to support visualization with the popular HoloViz stack of libraries (hvPlot, HoloViews, Datashader, etc.), UXarray provides methods for converting Grid and UxDataArray objects into a SpatialPandas GeoDataFrame, which can be used for visualizing the nodes, edges, and polygons that represent each grid, in addition to data variables.

Grid Conversion#

If you wish to plot only the grid’s geometrical elements such as nodes and edges (i.e. without any data variables mapped to them), you can directly convert a Grid instance into a GeoDataFrame.

Each UxDataset and UxDataArray has a Grid instance paired to it, which can be accessed using the .uxgrid attribute.

You can use the Grid.to_geodataframe() method to construct a GeoDataFrame.

gdf_grid = uxds.uxgrid.to_geodataframe()
gdf_grid
geometry
0 MultiPolygon([[[-45.0, -35.26438903808594, -42...
1 MultiPolygon([[[-42.0, -36.61769485473633, -39...
2 MultiPolygon([[[-39.0, -37.852420806884766, -3...
3 MultiPolygon([[[-36.0, -38.973445892333984, -3...
4 MultiPolygon([[[-33.0, -39.98556900024414, -30...
... ...
5395 MultiPolygon([[[147.3314971923828, 43.07368087...
5396 MultiPolygon([[[144.1993408203125, 42.01158142...
5397 MultiPolygon([[[141.0996856689453, 40.83760452...
5398 MultiPolygon([[[138.03317260742188, 39.5490798...
5399 MultiPolygon([[[135.0, 38.143314361572266, 132...

5400 rows × 1 columns

UxDataArray & UxDataset Conversion#

If you are interested in mapping data to each face, you can index the UxDataset with the variable of interest (in this case “psi”) to return the same GeoDataFrame as above, but now with data mapped to each face.

Warning

UXarray currently only supports mapping data variables that are mapped to each face. Variables that reside on the corner nodes or edges are currently not supported for visualization.

gdf_data = uxds["psi"].to_geodataframe()
gdf_data
geometry psi
0 MultiPolygon([[[-45.0, -35.26438903808594, -42... 1.351317
1 MultiPolygon([[[-42.0, -36.61769485473633, -39... 1.330915
2 MultiPolygon([[[-39.0, -37.852420806884766, -3... 1.310140
3 MultiPolygon([[[-36.0, -38.973445892333984, -3... 1.289056
4 MultiPolygon([[[-33.0, -39.98556900024414, -30... 1.267717
... ... ...
5395 MultiPolygon([[[147.3314971923828, 43.07368087... 0.733539
5396 MultiPolygon([[[144.1993408203125, 42.01158142... 0.712085
5397 MultiPolygon([[[141.0996856689453, 40.83760452... 0.690883
5398 MultiPolygon([[[138.03317260742188, 39.5490798... 0.669989
5399 MultiPolygon([[[135.0, 38.143314361572266, 132... 0.649467

5400 rows × 2 columns

Challenges with Representing Geoscience Data as Geometries#

When we convert to a GeoDataFrame, we internally represent the surface of a sphere as a collection of polygons over a 2D projection. However, this leads to issues around the Antimeridian (180 degrees east or west) where polygons are incorrectly constructed and have incorrect geometries. When constructing the GeoDataFrame, UXarray detects and corrects any polygon that touches or crosses the antimeridian. An array of indices of these faces can be accessed as part of the Grid object.

Antimeridian Visual
uxds.uxgrid.antimeridian_face_indices
array([1814, 1844, 1874, 1904, 1934, 1964, 1994, 2024, 2054, 2084, 2114,
       2144, 2174, 2204, 2234, 2264, 2294, 2324, 2354, 2384, 2414, 2444,
       2474, 2504, 2534, 2564, 2594, 2624, 2654, 2684, 3615, 3645, 3675,
       3705, 3735, 3765, 3795, 3825, 3855, 3885, 3915, 3945, 3975, 4005,
       4034, 4035, 4964, 4965, 4995, 5025, 5055, 5085, 5115, 5145, 5175,
       5205, 5235, 5265, 5295, 5325, 5355, 5385])

Taking a look at one of these faces that crosses or touches the antimeridian, we can see that it’s split across the antimeridian and represented as a MultiPolygon, which allows us to properly render this two dimensional grid.

gdf_data.geometry[uxds.uxgrid.antimeridian_face_indices[0]]
MultiPolygon([[[180.0, -42.0, 177.0, -41.96092987060547, 177.0, -44.96071243286133, 180.0, -45.0, 180.0, -42.0]], [[-180.0, -45.0, -180.0, -45.0, -180.0, -42.0, -180.0, -42.0, -180.0, -45.0]]])

For more details about the algorithm used for splitting these polygons, see the Antimeridian Python Package, which we use internally for correcting polygons on the antimeridian.

Visualizing Geometries#

The next sections will show how the GeoDataFrame representing our unstructured grid can be used with hvPlot and GeoViews to generate visualizations. Examples using both the matplotlib and bokeh backends are shown, with bokeh examples allowing for interactive plots.

Nodes#

# use matplotlib backend for rendering
hv.extension("matplotlib")

plot_kwargs = {
    "size": 6.0,
    "xlabel": "Longitude",
    "ylabel": "Latitude",
    "coastline": True,
    "width": 1600,
    "title": "Node Plot (Matplotlib Backend)",
}


gdf_grid.hvplot.points(**plot_kwargs)
# use bokeh backend for rendering
hv.extension("bokeh")

plot_kwargs = {
    "s": 1.0,
    "xlabel": "Longitude",
    "ylabel": "Latitude",
    "coastline": True,
    "frame_width": 700,
    "title": "Node Plot (Bokeh Backend)",
}

gdf_grid.hvplot.points(**plot_kwargs)