Spatial Selection Operators¶

In this tutorial, you’ll learn:¶

Using UXarray to select specific regions from an unstructured grid

Prerequisites¶

Concepts	Importance	Notes
UXDataset & UxDataArray Notebook	Necessary

Time to learn: 10 minutes

Overview¶

When working with unstructured grids, a region or zone rather than the entire global grid might be of interest for data analysis and visualization. Such spatial selection of the grid/data might help with not only analyzing/plotting a specific region/zone, but also reducing the data size and increasing the performance as well as allowing the entire plots to be effectively displayed on the screen.

This notebook showcases how to spatially subset a grid and take a cross-section of the grid at a constant latitude and longitude.

Subsetting¶

UXarray provides functionality for subsetting the grid into a bounding box, bounding circle, or K-nearest neighbors.

Before demonstrating them, let us load some global data first:

Load Data¶

# Import
import cartopy.crs as ccrs
import geoviews.feature as gf
import uxarray as ux

grid_path = "../../meshfiles/x1.655362.grid.nc"
data_path = "../../meshfiles/x1.655362.data.nc"

# Open dataset and grab a data variable of interest
uxds = ux.open_dataset(grid_path, data_path)
uxda = uxds["relhum_200hPa"][0]

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/micromamba/envs/unstructured-grid-viz-cookbook-dev/lib/python3.14/site-packages/xarray/core/dataset.py:1303, in Dataset._construct_dataarray(self, name)
   1301             variable = self._variables[name]
   1302         except KeyError:
-> 1303             _, name, variable = _get_virtual_variable(self._variables, name, self.sizes)
   1304 

KeyError: 'verticesOnEdge'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
File ~/micromamba/envs/unstructured-grid-viz-cookbook-dev/lib/python3.14/site-packages/xarray/core/dataset.py:1421, in Dataset.__getitem__(self, key)
   1419                 if isinstance(key, tuple):
   1420                     message += f"\nHint: use a list to select multiple variables, for example `ds[{list(key)}]`"
-> 1421                 raise KeyError(message) from e
   1422 

File ~/micromamba/envs/unstructured-grid-viz-cookbook-dev/lib/python3.14/site-packages/xarray/core/dataset.py:1303, in Dataset._construct_dataarray(self, name)
   1301             variable = self._variables[name]
   1302         except KeyError:
-> 1303             _, name, variable = _get_virtual_variable(self._variables, name, self.sizes)
   1304 

File ~/micromamba/envs/unstructured-grid-viz-cookbook-dev/lib/python3.14/site-packages/xarray/core/dataset_utils.py:79, in _get_virtual_variable(variables, key, dim_sizes)
     78 if len(split_key) != 2:
---> 79     raise KeyError(key)
     81 ref_name, var_name = split_key

KeyError: 'verticesOnEdge'

The above exception was the direct cause of the following exception:

KeyError                                  Traceback (most recent call last)
Cell In[1], line 10
      6 grid_path = "../../meshfiles/x1.655362.grid.nc"
      7 data_path = "../../meshfiles/x1.655362.data.nc"
      8 
      9 # Open dataset and grab a data variable of interest
---> 10 uxds = ux.open_dataset(grid_path, data_path)
     11 uxda = uxds["relhum_200hPa"][0]

File ~/micromamba/envs/unstructured-grid-viz-cookbook-dev/lib/python3.14/site-packages/uxarray/core/api.py:151, in open_dataset(grid_filename_or_obj, filename_or_obj, latlon, use_dual, grid_kwargs, **kwargs)
    146     warn('source_grid is no longer a supported kwarg',
    147          DeprecationWarning,
    148          stacklevel=2)
    150 # Grid definition
--> 151 uxgrid = open_grid(grid_filename_or_obj,
    152                    latlon=latlon,
    153                    use_dual=use_dual,
    154                    **grid_kwargs)
    156 # UxDataset
    157 ds = xr.open_dataset(filename_or_obj, decode_times=False,
    158                      **kwargs)  # type: ignore

File ~/micromamba/envs/unstructured-grid-viz-cookbook-dev/lib/python3.14/site-packages/uxarray/core/api.py:77, in open_grid(grid_filename_or_obj, latlon, use_dual, **kwargs)
     72 elif isinstance(grid_filename_or_obj, (str, Path, PurePath)):
     73     grid_ds = xr.open_dataset(grid_filename_or_obj,
     74                               decode_times=False,
     75                               **kwargs)
---> 77     uxgrid = Grid.from_dataset(grid_ds, use_dual=use_dual)
     79 elif isinstance(grid_filename_or_obj,
     80                 (list, tuple, np.ndarray, xr.DataArray)):
     81     uxgrid = Grid.from_face_vertices(grid_filename_or_obj, latlon=latlon)

File ~/micromamba/envs/unstructured-grid-viz-cookbook-dev/lib/python3.14/site-packages/uxarray/grid/grid.py:148, in Grid.from_dataset(cls, dataset, use_dual)
    146     grid_ds, source_dims_dict = _read_ugrid(dataset)
    147 elif source_grid_spec == "MPAS":
--> 148     grid_ds, source_dims_dict = _read_mpas(dataset, use_dual=use_dual)
    149 elif source_grid_spec == "Shapefile":
    150     raise ValueError("Shapefiles not yet supported")

File ~/micromamba/envs/unstructured-grid-viz-cookbook-dev/lib/python3.14/site-packages/uxarray/io/_mpas.py:401, in _read_mpas(ext_ds, use_dual)
    398     source_dim_map = _dual_to_ugrid(ext_ds, ds)
    399 # convert primal-mesh to UGRID
    400 else:
--> 401     source_dim_map = _primal_to_ugrid(ext_ds, ds)
    403 return ds, source_dim_map

File ~/micromamba/envs/unstructured-grid-viz-cookbook-dev/lib/python3.14/site-packages/uxarray/io/_mpas.py:102, in _primal_to_ugrid(in_ds, out_ds)
     92 out_ds["Mesh2_face_nodes"] = xr.DataArray(
     93     data=verticesOnCell,
     94     dims=["nMesh2_face", "nMaxMesh2_face_nodes"],
   (...)     98         "start_index": INT_DTYPE(0)
     99     })
    101 # vertex indices that saddle a given edge
--> 102 verticesOnEdge = np.array(in_ds['verticesOnEdge'].values, dtype=INT_DTYPE)
    104 # replace missing/zero values with fill value
    105 verticesOnEdge = _replace_zeros(verticesOnEdge)

File ~/micromamba/envs/unstructured-grid-viz-cookbook-dev/lib/python3.14/site-packages/xarray/core/dataset.py:1421, in Dataset.__getitem__(self, key)
   1417 
   1418                 # If someone attempts `ds['foo' , 'bar']` instead of `ds[['foo', 'bar']]`
   1419                 if isinstance(key, tuple):
   1420                     message += f"\nHint: use a list to select multiple variables, for example `ds[{list(key)}]`"
-> 1421                 raise KeyError(message) from e
   1422 
   1423         if utils.iterable_of_hashable(key):
   1424             return self._copy_listed(key)

KeyError: "No variable named 'verticesOnEdge'. Did you mean one of ('verticesOnCell',)?"

Plot The Global Data¶

Note!

The visualizations throughout this tutorial are only for demonstrating the results of the subsetting and cross-sections versus the global grid. Since the details of plotting with UXarray will be covered in the next chapter, we will not go over any details of the plots here.

plot_opts = {"width": 700, "height": 350}

features = gf.coastline(
    projection=ccrs.PlateCarree(), line_width=0.4, scale="50m"
) * gf.states(projection=ccrs.PlateCarree(), line_width=0.4, scale="50m")

clim = (uxda.values.min(), uxda.values.max())

uxda.plot(
    rasterize=True, periodic_elements="exclude", title="Global Grid", **plot_opts
) * features

Generate Various Subsets¶

# Bounding box around Boulder, CO
ref_lon = -105.2705
ref_lat = 40.0150
ref_offset = 1

lon_bounds = (ref_lon - ref_offset, ref_lon + ref_offset)
lat_bounds = (ref_lat - ref_offset, ref_lat + ref_offset)

# Subset the global data variable and its grid using the bounding box
uxda_bbox = uxda.subset.bounding_box(lon_bounds, lat_bounds)

# Now use a bounding circle subsetting
ref_center = [ref_lon, ref_lat]

uxda_bcircle = uxda.subset.bounding_circle(ref_center, ref_offset)

# Now use a K-nearest neighbor subsetting
k_nn = 60

uxda_nn = uxda.subset.nearest_neighbor(ref_center, k=k_nn, element="nodes")

Plot The Subsets¶

plot_bbox = (
    uxda_bbox.plot(
        rasterize=True,
        periodic_elements="exclude",
        clim=clim,
        title="Bounding Box Subset around Boulder, CO (Corner Node Query)",
        **plot_opts,
    )
    * features
)

plot_bcircle = (
    uxda_bcircle.plot(
        rasterize=True,
        periodic_elements="exclude",
        clim=clim,
        title="Bounding Circle Subset around Boulder, CO (Corner Node Query)",
        **plot_opts,
    )
    * features
)

plot_nn = (
    uxda_nn.plot(
        rasterize=True,
        periodic_elements="exclude",
        clim=clim,
        title="K-Nearest Neighbor Subset around Boulder, CO (Corner Node Query)",
        **plot_opts,
    )
    * features
)

(plot_bbox + plot_bcircle + plot_nn).cols(1)

Cross-sections¶

Similarly, UXarray provides functionality for taking cross-sections of a grid/data at constant latitudes and longitudes.

Load Data¶

# Data paths
grid_path = "../../meshfiles/outCSne30.grid.ug"
data_path = "../../meshfiles/outCSne30.data.nc"

# Open dataset and grab a data variable of interest
uxds = ux.open_dataset(grid_path, data_path)
uxda = uxds["psi"]

Plot The Global Data¶

projection = ccrs.Robinson()

features_2 = gf.coastline(
    projection=projection, line_width=0.4, scale="50m"
) * gf.states(projection=projection, line_width=0.4, scale="50m")

uxda.plot(
    cmap="inferno",
    periodic_elements="split",
    projection=projection,
    title="Global Plot",
) * features_2

Generate Cross-Sections¶

# Cross-section at the latitude of Boulder, CO
cross_lat = 40.0150

uxda_cross_lat = uxda.cross_section.constant_latitude(cross_lat)

# Now cross-section at the longitude of Boulder, CO
cross_lon = -105.2705

uxda_cross_lon = uxda.cross_section.constant_longitude(cross_lon)

Plot The Cross-Sections¶

plot_cross_lat = (
    uxda_cross_lat.plot(
        rasterize=False,
        backend="bokeh",
        cmap="inferno",
        projection=projection,
        global_extent=True,
        coastline=True,
        title=f"Cross-section at ({cross_lat}) degrees latitude around Boulder, CO",
    )
    * features_2
)

plot_cross_lon = (
    uxda_cross_lon.plot(
        rasterize=False,
        backend="bokeh",
        cmap="inferno",
        projection=projection,
        global_extent=True,
        coastline=True,
        title=f"Cross-section at ({cross_lon}) degrees longitude around Boulder, CO",
    )
    * features_2
)

(plot_cross_lat + plot_cross_lon).cols(1)

What is next?¶

With this section, we have wrapped up this chapter, move on to the next chapter!

Spatial Selection Operators

Spatial Selection Operators¶

In this tutorial, you’ll learn:¶

Related Documentation¶

Prerequisites¶

Overview¶

Subsetting¶

Load Data¶

Plot The Global Data¶

Generate Various Subsets¶

Plot The Subsets¶

Cross-sections¶

Load Data¶

Plot The Global Data¶

Generate Cross-Sections¶

Plot The Cross-Sections¶

What is next?¶