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 1.351 1.331 1.31 1.289 ... 0.6909 0.67 0.6495uxds.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.
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)
Edges#
# use matplotlib backend for rendering
hv.extension("matplotlib")
plot_kwargs = {"linewidth": 1.0, "xlabel":" Longitude", "ylabel": "Latitude", "coastline": True, "width": 1600 , "title": "Edge Plot (Matplotlib Backend)", "color": "black"}
import cartopy.crs as ccrs
gdf_grid.hvplot.paths(**plot_kwargs)
# use bokeh backend for rendering
hv.extension("bokeh")
plot_kwargs = {"line_width": 0.5, "xlabel": "Longitude", "ylabel": "Latitude", "coastline": True, "frame_width": 700, "title": "Edge Plot (Bokeh Backend)"}
gdf_grid.hvplot.paths(**plot_kwargs)
Visualizing Data Variables#
hv.extension("matplotlib")
plot_kwargs = {"c": "psi", "cmap": "inferno", "width": 400, "height": 200, "title": "Filled Polygon Plot (Matplotlib Backend, Rasterized)", "xlabel":" Longitude", "ylabel": "Latitude"}
gdf_data.hvplot.polygons(**plot_kwargs, rasterize=True)
Note
Visualizing filled polygons without rasterization using the matplotlib backend produces incorrect results, see hvplot/#1099
hv.extension("bokeh")
plot_kwargs = {"c": "psi", "cmap": "inferno", "line_width": 0.1, "frame_width": 500, "frame_height": 250, "xlabel":" Longitude", "ylabel": "Latitude"}
gdf_data.hvplot.polygons(**plot_kwargs, rasterize=True)
hv.Layout(gdf_data.hvplot.polygons(**plot_kwargs, rasterize=True).opts(title="Filled Polygon Plot (Bokeh Backend, Rasterized)") +
gdf_data.hvplot.polygons(**plot_kwargs, rasterize=False).opts(title="Filled Polygon Plot (Bokeh Backend, Vector)")).cols(1)
Geographic Projections#
hv.extension("matplotlib")
plot_kwargs = {"c": "psi", "cmap": "inferno", "width": 400, "height": 200, "title": "Filled Polygon Plot (Matplotlib Backend, Rasterized, Orthographic Projection)", "xlabel":" Longitude", "ylabel": "Latitude"}
gdf_data.hvplot.polygons(**plot_kwargs, rasterize=True, projection=ccrs.Orthographic())
hv.extension("matplotlib")
plot_kwargs = {"c": "psi", "cmap": "inferno", "width": 400, "height": 200, "title": "Filled Polygon Plot (Matplotlib Backend, Rasterized, Robinson Projection)", "xlabel":" Longitude", "ylabel": "Latitude"}
gdf_data.hvplot.polygons(**plot_kwargs, rasterize=True, projection=ccrs.Robinson())