Select unstructured data#

Xarray has flexible tools for label based selection, in the form of .sel and .isel for index selection. This works well for structured data since the orthogonality of the x and y axes is reflected in the axes of the underlying arrays. This orthogonality does not exist for unstructured grids, as the data for all faces cannot be stored in a two-dimensional array and is stored in a one-dimensional array instead.

Xugrid provides tools for convenient spatial selection, primarily via the .ugrid.sel method; its behavior is comparable to xarray’s .sel method. The .ugrid.sel method should only be used for selection in the x or y dimension. Selections along other dimension (such as time) should be performed by xarray’s .sel instead (without the ugrid accessor).

The examples below demonstrate the various ways to select data.

Imports#

The following imports suffice for the examples.

import matplotlib.pyplot as plt
import numpy as np
import shapely

import xugrid as xu

We will take a look at a sample dataset: a triangular grid with the surface elevation of the Netherlands.

uda = xu.data.elevation_nl()
uda.ugrid.plot(vmin=-20, vmax=90, cmap="terrain")
selection
<matplotlib.collections.PolyCollection object at 0x7f60488c9010>

We will start by demonstrating the behavior of .ugrid.sel. This method takes several types of arguments, like its xarray equivalent. The return type and shape of the selection operation depends on the argument given.

Selection

Result type

Subset

xugrid

Point

xarray

Line

xarray

Grid subset selection#

A subset of the unstructured grid is returned by using slices without a step:

subset = uda.ugrid.sel(x=slice(100_000.0, 200_000.0), y=slice(450_000.0, 550_000.0))
subset.ugrid.plot(vmin=-20, vmax=90, cmap="terrain")
selection
<matplotlib.collections.PolyCollection object at 0x7f6047833410>

The default arguments of x and y are: slice(None, None). In such a case the entire grid is returned.

subset = uda.ugrid.sel()
subset.ugrid.plot(vmin=-20, vmax=90, cmap="terrain")
selection
<matplotlib.collections.PolyCollection object at 0x7f60488c6390>

Note

None in a Python slice can be interpreted as “from the start” or “up to and including the end”.

This means we can easily select along a single dimension:

subset = uda.ugrid.sel(x=slice(100_000.0, 200_000.0))
subset.ugrid.plot(vmin=-20, vmax=90, cmap="terrain", aspect=1, size=5)
selection
<matplotlib.collections.PolyCollection object at 0x7f6052108c20>

Or, using None if we only care about the start:

subset = uda.ugrid.sel(x=slice(100_000.0, None))
subset.ugrid.plot(vmin=-20, vmax=90, cmap="terrain", aspect=1, size=5)
selection
<matplotlib.collections.PolyCollection object at 0x7f604791bd40>

Point selection#

Since point data can be represented as an ordinary xarray DataArray with x and y coordinates, all point selection result in xarray DataArrays rather than UgridDataArrays with an associated unstructured grid topology.

We will use a utility function to show what is selected on the map:

def show_point_selection(uda, da):
    _, ax = plt.subplots()
    uda.ugrid.plot(ax=ax, vmin=-20, vmax=90, cmap="terrain")
    ax.scatter(da["mesh2d_x"], da["mesh2d_y"], color="red")
    ax.set_aspect(1.0)

Two values will select a point:

da = uda.ugrid.sel(x=150_000.0, y=463_000.0)
show_point_selection(uda, da)
da
selection
<xarray.DataArray 'elevation' (mesh2d_nFaces: 1)> Size: 4B
[1 values with dtype=float32]
Coordinates:
    mesh2d_face_x  (mesh2d_nFaces) float64 8B ...
    mesh2d_face_y  (mesh2d_nFaces) float64 8B ...
  * mesh2d_nFaces  (mesh2d_nFaces) int64 8B 4221
    mesh2d_index   (mesh2d_nFaces) int64 8B 0
    mesh2d_x       (mesh2d_nFaces) float64 8B 1.5e+05
    mesh2d_y       (mesh2d_nFaces) float64 8B 4.63e+05
Attributes:
    unit:     m NAP


Multiple values are broadcasted against each other (“outer indexing”). If we select by three x values and two y values, the result is a collection of six points:

da = uda.ugrid.sel(x=[125_000.0, 150_000.0, 175_000.0], y=[400_000.0, 465_000.0])
show_point_selection(uda, da)
da
selection
<xarray.DataArray 'elevation' (mesh2d_nFaces: 6)> Size: 24B
[6 values with dtype=float32]
Coordinates:
    mesh2d_face_x  (mesh2d_nFaces) float64 48B ...
    mesh2d_face_y  (mesh2d_nFaces) float64 48B ...
  * mesh2d_nFaces  (mesh2d_nFaces) int64 48B 1905 1356 372 3057 4198 4113
    mesh2d_index   (mesh2d_nFaces) int64 48B 0 1 2 3 4 5
    mesh2d_x       (mesh2d_nFaces) float64 48B 1.25e+05 1.5e+05 ... 1.75e+05
    mesh2d_y       (mesh2d_nFaces) float64 48B 4e+05 4e+05 ... 4.65e+05 4.65e+05
Attributes:
    unit:     m NAP


To select points without broadcasting, use .ugrid.sel_points instead:

da = uda.ugrid.sel_points(
    x=[125_000.0, 150_000.0, 175_000.0], y=[400_000.0, 430_000.0, 465_000.0]
)
show_point_selection(uda, da)
da
selection
<xarray.DataArray 'elevation' (mesh2d_nFaces: 3)> Size: 12B
[3 values with dtype=float32]
Coordinates:
    mesh2d_face_x  (mesh2d_nFaces) float64 24B ...
    mesh2d_face_y  (mesh2d_nFaces) float64 24B ...
  * mesh2d_nFaces  (mesh2d_nFaces) int64 24B 1905 2675 4113
    mesh2d_index   (mesh2d_nFaces) int64 24B 0 1 2
    mesh2d_x       (mesh2d_nFaces) float64 24B 1.25e+05 1.5e+05 1.75e+05
    mesh2d_y       (mesh2d_nFaces) float64 24B 4e+05 4.3e+05 4.65e+05
Attributes:
    unit:     m NAP


We can sample points along a line as well by providing slices with a step:

da = uda.ugrid.sel(x=slice(100_000.0, 200_000.0, 10_000.0), y=465_000.0)
show_point_selection(uda, da)
da
selection
<xarray.DataArray 'elevation' (mesh2d_nFaces: 10)> Size: 40B
[10 values with dtype=float32]
Coordinates:
    mesh2d_face_x  (mesh2d_nFaces) float64 80B ...
    mesh2d_face_y  (mesh2d_nFaces) float64 80B ...
  * mesh2d_nFaces  (mesh2d_nFaces) int64 80B 3263 3225 3144 ... 2941 3135 4256
    mesh2d_index   (mesh2d_nFaces) int64 80B 0 1 2 3 4 5 6 7 8 9
    mesh2d_x       (mesh2d_nFaces) float64 80B 1e+05 1.1e+05 ... 1.8e+05 1.9e+05
    mesh2d_y       (mesh2d_nFaces) float64 80B 4.65e+05 4.65e+05 ... 4.65e+05
Attributes:
    unit:     m NAP


Two slices with a step results in broadcasting:

da = uda.ugrid.sel(
    x=slice(100_000.0, 200_000.0, 10_000.0), y=slice(400_000.0, 500_000.0, 10_000.0)
)
show_point_selection(uda, da)
da
selection
<xarray.DataArray 'elevation' (mesh2d_nFaces: 100)> Size: 400B
[100 values with dtype=float32]
Coordinates:
    mesh2d_face_x  (mesh2d_nFaces) float64 800B ...
    mesh2d_face_y  (mesh2d_nFaces) float64 800B ...
  * mesh2d_nFaces  (mesh2d_nFaces) int64 800B 1867 1882 1890 ... 4121 4077 1581
    mesh2d_index   (mesh2d_nFaces) int64 800B 0 1 2 3 4 5 ... 94 95 96 97 98 99
    mesh2d_x       (mesh2d_nFaces) float64 800B 1e+05 1.1e+05 ... 1.9e+05
    mesh2d_y       (mesh2d_nFaces) float64 800B 4e+05 4e+05 ... 4.9e+05 4.9e+05
Attributes:
    unit:     m NAP


As well as a slice with a step and multiple values:

da = uda.ugrid.sel(x=slice(100_000.0, 200_000.0, 10_000.0), y=[400_000.0, 430_000.0])
show_point_selection(uda, da)
da
selection
<xarray.DataArray 'elevation' (mesh2d_nFaces: 20)> Size: 80B
[20 values with dtype=float32]
Coordinates:
    mesh2d_face_x  (mesh2d_nFaces) float64 160B ...
    mesh2d_face_y  (mesh2d_nFaces) float64 160B ...
  * mesh2d_nFaces  (mesh2d_nFaces) int64 160B 1867 1882 1890 ... 2889 226 3035
    mesh2d_index   (mesh2d_nFaces) int64 160B 0 1 2 3 4 5 ... 14 15 16 17 18 19
    mesh2d_x       (mesh2d_nFaces) float64 160B 1e+05 1.1e+05 ... 1.9e+05
    mesh2d_y       (mesh2d_nFaces) float64 160B 4e+05 4e+05 ... 4.3e+05 4.3e+05
Attributes:
    unit:     m NAP


Line selection#

Since line data can be represented as an ordinary xarray DataArray with x and y coordinates, all line selection result in xarray DataArrays rather than UgridDataArrays with an associated unstructured grid topology.

Line selection is performed by finding all faces that are intersected by the line.

We start by defining a utility to show the selection again:

def show_line_selection(uda, da, line_x=None, line_y=None):
    _, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 5))
    uda.ugrid.plot(ax=ax0, vmin=-20, vmax=90, cmap="terrain")
    da.plot(ax=ax1, x="mesh2d_s")
    if line_x is None:
        ax0.axhline(line_y, color="red")
    elif line_y is None:
        ax0.axvline(line_x, color="red")
    else:
        ax0.plot(line_x, line_y, color="red")
    ax0.set_aspect(1.0)

A single value for either x or y in .ugrid.sel will select values along a line:

da = uda.ugrid.sel(y=465_000.0)
show_line_selection(uda, da, line_y=465_000.0)
selection

Line segments that are not axis aligned can be selected with .ugrid.intersect_line:

da = uda.ugrid.intersect_line(start=(60_000.0, 400_000.0), end=(190_000.0, 475_000.0))
show_line_selection(uda, da, (60_000.0, 190_000.0), (400_000.0, 475_000.0))
selection

Linestrings can be selected with .ugrid.intersect_linestring:

linestring = shapely.geometry.LineString(
    [
        (60_000.0, 400_000.0),
        (190_000.0, 400_000.0),
        (120_000.0, 575_000.0),
        (250_000.0, 575_000.0),
    ]
)

da = uda.ugrid.intersect_linestring(linestring)
show_line_selection(uda, da, *shapely.get_coordinates(linestring).T)
selection

This will work for any type of shapely line:

ring = shapely.geometry.Point(155_000.0, 463_000).buffer(50_000.0).exterior
da = uda.ugrid.intersect_linestring(ring)
show_line_selection(uda, da, *shapely.get_coordinates(ring).T)
selection

Index selection#

We may also use ordinary index selection to create a subset. This does not require the .ugrid accessor. For example, to take only the first thousands faces:

subset = uda.isel(mesh2d_nFaces=np.arange(1000))
subset.ugrid.plot(vmin=-20, vmax=90, cmap="terrain", aspect=1, size=5)
selection
<matplotlib.collections.PolyCollection object at 0x7f6047a5bd40>

For a 2D topology, selecting faces by an index always results in a valid topology. However, selecting by node or edge does not give a guarantee that the result forms a valid 2D topology: e.g. if we only select two nodes, or only two edges from a face, the result cannot form a valid 2D face.

To avoid generating invalid topologies, xugrid always checks whether the result of a selection results in a valid 2D topology and raises an error if the result is invalid.

In general, index selection should only be performed on the “core” dimension of the UGRID topology. This is the edge dimension for 1D topologies, and the face dimension for 2D topologies.

Total running time of the script: (0 minutes 3.487 seconds)

Gallery generated by Sphinx-Gallery