Effeciantly subset existing ugrid to a smaller ugrid given a bbox or polygon

[saeed-moghimi-noaa]

Hi There,

Is there a suggested tool to perform efficient subset of a large ugrid data set to a new ugrid (xarray) with the consistent connectivity given a bbox or polygon as an input using xugrid tool set? I would like to retain all the properties of the parent mesh such as node coordinates and connectivity and values.

something like:

new_ugrid = subset (large_input_ugrid, bbox)

Thanks.

[Huite]

Hi @saeed-moghimi-noaa

The user guide part on selection should give relevant information on this: https://deltares.github.io/xugrid/examples/selection.html

The .ugrid.sel method will return subsets with valid topologies.

Note that this currently selects by the face coordinates, which may not be what you desire. See issue #182 for a workaround (although to be clear: we should tackle it with a proper solution so it doesn't require a workaround).

Selecting by polygon can be done relatively easily, but also requires more on than step. I'll open an issue for the functionality.

For now, one could use this multi-step approach:

import xugrid as xu
import geopandas as gpd

# Read a polygon into a geodataframe
gdf = gdp.read_file("my-polygon.shp")

# Read the ugrid file containing a single variable
uda = xu.open_dataarray("my-ugrid.nc")

# Burn the vector geometry in to the unstructured topology
burned = xu.burn_vector_geometry(gdf, uda, column="id")  # let's assume there's an id colum with a value of 1

# Use the burned result to select only the values in the polygon, then drop the faces outside of the polygon
selection = uda.where(burned == 1).dropna(uda.ugrid.grid.face_dimension)

There's a bit more info on "burning" vector geometries in this part of the user guide: https://deltares.github.io/xugrid/examples/vector_conversion.html

In all cases, this results in a new, consistent, geometry. In case you want to preserve the original topology, you would use .where instead and "mask" all values outside of the selection by NaN.

[mgeraeds]

Hi Huite,

I have had the same question and applied this same concept as you state here, based on the earlier examples of burn_vector_geometry. I'd like to apply this principle to a complete dataset, and that doesn't seem to work (yet). I've done it exactly in the same way as mentioned above, the only difference is the last step - I do selection= uds.where(burned == 1, drop=True), since doing exactly what you describe above leads to out-of-memory errors very fast (6 GB dataset on a supercomputer with 336GB RAM). Seems like there's a memory leak or something, but have no clue why or where. I'm using a Dask SLURMCluster.

For variables on the faces, this works, which is great. However, it does not work for variables on the edges, it makes a new dimension mesh2d_nFaces. It would make sense to me that this function would be made to work for an entire dataset (including the variables that are not located on the faces) and I think it would be really useful, as it's difficult to unify the topology along all of the other non-face-based variables) otherwise. Or at least, as far as I know. Is there any way to do this as of now?

[Huite]

Hi @mgeraeds,

I've done it exactly in the same way as mentioned above, the only difference is the last step - I do selection= uds.where(burned == 1, drop=True)

Ah, good to know that that works better. I'm guessing xarray/dask is able to delay the operation if you call the drop immediately. You could also do .isel after extracting the face indices from the burned polygon data. That is guaranteed to be proper lazy evaluation (if your data is chunked / dask backed). I.e.:

selection = uda.isel({uda.ugrid.grid.face_dimension: burned.notnull()})

That might be the best, and I should maybe update my answer above.

Regarding edges: the edges aren't aware of the burned faces. So we'd have to associate the face data to the edge data somehow. Basically, you'd want #173.

Since that doesn't exist yet, you'll have to do it yourself. Here's a basic example with some of the sample data:

import xugrid as xu
import xarray as xr


provinces = xu.data.provinces_nl().to_crs(epsg=28992)
uda = xu.data.elevation_nl()

utrecht = provinces.loc[provinces["name"] == "Utrecht"]
burned = xu.burn_vector_geometry(utrecht, uda, column="level")

grid = uda.ugrid.grid
indexer = xr.DataArray(grid.edge_face_connectivity, dims=(grid.edge_dimension, "two"))
burned_edges = burned.isel({grid.face_dimension: indexer})
condition = burned_edges.notnull().any("two")

selection = uda.where(condition, drop=True)

Dependent on which edges you're interested, maybe you want condition = burned_edges.notnull().all("two") instead.

To get a feel for it, this is condition.ugrid.plot(cmap="turbo"):

Ideally in the future you'd be able to do:

condition = burned.ugrid.to_edge().notnull().any()

But that will require some thinking.

[JoerivanEngelen]

@saeed-moghimi-noaa
FYI: I converted this to a Q&A, as your question will also be interesting to other users.
Could you confirm whether your question is answered sufficiently by @Huite's reply?