Subsetting an Unstructured Grid: Analysis Over Chicago#
Authors: Philip Chmielowiec
This usage example showcases various ways of subsetting an unstructured grid using UXarray, focussing on analyzing a region around Chicago, Illinois.
import uxarray as ux
import geoviews.feature as gf
import cartopy.crs as ccrs
import holoviews as hv
import warnings
import geocat.datafiles as geodf
plot_opts = {"width" : 700, "height": 350}
hv.extension('bokeh')
warnings.filterwarnings("ignore")
Downloading file 'registry.txt' from 'https://github.com/NCAR/GeoCAT-datafiles/raw/main/registry.txt' to '/home/docs/.cache/geocat'.
Setup#
In this example, we will be using the geocat-datafiles package to obtain our grid and data files.
The dataset used in this example is a 30km global MPAS meshes. We will be investigating the relative humidity vertically interpolated to 200hPa (relhum200hPa) data variable.
datafiles = (
geodf.get(
"netcdf_files/MPAS/FalkoJudt/dyamond_1/30km/diag.2016-08-20_00.00.00_subset.nc"
),
geodf.get("netcdf_files/MPAS/FalkoJudt/dyamond_1/30km/x1.655362.grid_subset.nc"),
)
Downloading file 'netcdf_files/MPAS/FalkoJudt/dyamond_1/30km/diag.2016-08-20_00.00.00_subset.nc' from 'https://github.com/NCAR/GeoCAT-datafiles/raw/main/netcdf_files/MPAS/FalkoJudt/dyamond_1/30km/diag.2016-08-20_00.00.00_subset.nc' to '/home/docs/.cache/geocat'.
Downloading file 'netcdf_files/MPAS/FalkoJudt/dyamond_1/30km/x1.655362.grid_subset.nc' from 'https://github.com/NCAR/GeoCAT-datafiles/raw/main/netcdf_files/MPAS/FalkoJudt/dyamond_1/30km/x1.655362.grid_subset.nc' to '/home/docs/.cache/geocat'.
uxds = ux.open_dataset(datafiles[1], datafiles[0])
clim = (uxds['relhum_200hPa'][0].values.min(), uxds['relhum_200hPa'][0].values.max())
features = gf.coastline(projection=ccrs.PlateCarree(), line_width=1, scale='50m') * \
gf.states(projection=ccrs.PlateCarree(), line_width=1, scale='50m')
Global Grid#
Many unstructured grids, such as those from global climate models, span the entire surface of a sphere (both with or without masks, such as continents).
UXarray supports working with these global grids, handling cases that arise with the spherical geometry of the earth (wrap around at the antimeridian, pole points, etc.)
uxds['relhum_200hPa'][0].plot.rasterize(method='polygon',
exclude_antimeridian=True,
title = "Global Grid",
**plot_opts) * features
In addition to plotting global grids, we can perform analysis operations.
uxds['relhum_200hPa'][0].values.mean()
46.819023
Regional Subsets#
UXarray supports taking subsets of a grid, which allows us to select a region and perform analysis directly on that area, as opposed to the global grid.
There are currently three supported subsetting methods, both for the Grid and UxDataArray data structures.
uxds['relhum_200hPa'].subset
<uxarray.UxDataArray.subset>
Supported Methods:
* nearest_neighbor(center_coord, k, element, **kwargs)
* bounding_circle(center_coord, r, element, **kwargs)
* bounding_box(lon_bounds, lat_bounds, element, method, **kwargs)
uxds['relhum_200hPa'].uxgrid.subset
<uxarray.Grid.subset>
Supported Methods:
* nearest_neighbor(center_coord, k, element, **kwargs)
* bounding_circle(center_coord, r, element, **kwargs)
* bounding_box(lon_bounds, lat_bounds, element, method, **kwargs)
Bounding Box#
We can declare a bounding box centered about the Chicago area by specifying the minimum and maximum longitude and latitude bounds.
lon_bounds = (-87.6298 - 2, -87.6298 + 2)
lat_bounds = (41.8781 - 2, 41.8781 + 2)
Our bounding box ensures that the coordinates of our select element (nodes, edge_centers, or face_centers) are within the defined bounding box range.
Below is an example using the corner nodes for our subset.
bbox_subset_nodes = uxds['relhum_200hPa'][0].subset.bounding_box(lon_bounds, lat_bounds, element='nodes')
bbox_subset_nodes.plot.rasterize(method='polygon',
exclude_antimeridian=True,
clim=clim,
title="Bounding Box Subset",
**plot_opts) * features
And similarly using face centers.
bbox_subset_faces= uxds['relhum_200hPa'][0].subset.bounding_box(lon_bounds, lat_bounds, element='face centers')
bbox_subset_faces.plot.rasterize(method='polygon',
exclude_antimeridian=True,
clim=clim,
title="Bounding Box Subset (Face Center Query)",
**plot_opts) * features
While the bounding box is generally the same, you will notice differences along the border depending on which element is used to query.
Note
Specifying which element to query (i.e. nodes, edgecenters, or face centers) is supported by all subsetting methods.
Bounding Circle#
A bounding circle is defined using a center coordinate (lon, lat) and a radius (in degrees). The resulting subset will contain all elements within the radius of that circle.
center_coord = [-87.6298, 41.8781]
r = 2
bcircle_subset = uxds['relhum_200hPa'][0].subset.bounding_circle(center_coord, r)
bcircle_subset.plot.rasterize(method='polygon',
exclude_antimeridian=True,
clim=clim,
title="Bounding Circle Subset",
**plot_opts) * features
Nearest Neighbor#
Similar to the bounding circle, we can perform a nearest neighbor subset at some center coordinate (lon, lat) and query for some number of elements k
center_coord = [-87.6298, 41.8781]
nn_subset = uxds['relhum_200hPa'][0].subset.nearest_neighbor(center_coord, k = 30, element='nodes')
nn_subset.plot.rasterize(method='polygon', exclude_antimeridian=True,
clim=clim,
title= "Nearest Neighbor Subset",
**plot_opts) * features
nn_subset_120 = uxds['relhum_200hPa'][0].subset.nearest_neighbor(center_coord, k = 120, element="face centers")
nn_subset_120.plot.rasterize(method='polygon',
exclude_antimeridian=True,
clim=clim,
title= "Nearest Neighbor Subset (120 Faces)",
**plot_opts) * features
nn_subset_1 = uxds['relhum_200hPa'][0].subset.nearest_neighbor(center_coord, k = 1, element="face centers")
nn_subset_1.plot.rasterize(method='polygon', exclude_antimeridian=True,
clim=clim,
title= "Nearest Neighbor Subset (Closest Face)",
**plot_opts) * features
Analysis Operators#
Since each subset is a newly initialized UxDataArray, paired also with a newly initialized Grid, we can perform analysis operators directly on these new objects.
Below is a few examples of basic statical operations on the subset data arrays.
bbox_subset_nodes.values.mean(), bbox_subset_faces.values.mean(), bcircle_subset.values.mean()
(88.41318, 89.05569, 89.83152)
bbox_subset_nodes.values.std(), bbox_subset_faces.values.std(), bcircle_subset.values.std()
(21.11058, 20.208374, 19.313942)
bbox_subset_nodes.values.min(), bbox_subset_faces.values.min(), bcircle_subset.values.min()
(34.943027, 35.723145, 35.723145)
bbox_subset_nodes.values.max(), bbox_subset_faces.values.max(), bcircle_subset.values.max()
(121.48339, 121.48339, 121.48339)