Clip unstructured grid (ugrid)

Clip unstructured grid (ugrid)#

import CHM as pc
import uxarray as ux
# Large number of triangles (large meshes) needs out of core dask support

from dask_jobqueue import PBSCluster
from dask.distributed import Client

cluster = PBSCluster(
    cores=80,
    processes=40,
    memory='185GB',
    interface='ib0',
    local_directory='$TMPDIR',
    queue='development',
    walltime='06:00:00',
    job_script_prologue=["spack env activate analysis"]
)
print(cluster.job_script())
cluster.scale(jobs=25)
client = Client(cluster)
client
# ugrid expects mesh and vars to be seperate
ds = ux.open_mfdataset('output_2016.nc', 'output_2016.nc')
ds
<xarray.UxDataset> Size: 12TB
Dimensions:                (n_face: 28495325, n_max_face_nodes: 3,
                            n_node: 14638167, time: 8737)
Coordinates:
  * time                   (time) datetime64[ns] 70kB 2016-10-01 ... 2017-09-30
    Mesh2_face_x           (n_face) float64 228MB dask.array<chunksize=(28495325,), meta=np.ndarray>
    Mesh2_face_y           (n_face) float64 228MB dask.array<chunksize=(28495325,), meta=np.ndarray>
Dimensions without coordinates: n_face, n_max_face_nodes, n_node
Data variables: (12/51)
    Mesh2                  int32 4B ...
    Mesh2_face_nodes       (n_face, n_max_face_nodes) uint32 342MB dask.array<chunksize=(28495325, 3), meta=np.ndarray>
    Mesh2_node_x           (n_node) float64 117MB dask.array<chunksize=(14638167,), meta=np.ndarray>
    Mesh2_node_y           (n_node) float64 117MB dask.array<chunksize=(14638167,), meta=np.ndarray>
    Mesh2_node_z           (n_node) float64 117MB dask.array<chunksize=(14638167,), meta=np.ndarray>
    Mesh2_node_z_paraview  (n_node) float64 117MB dask.array<chunksize=(14638167,), meta=np.ndarray>
    ...                     ...
    param_Ninja7_V         (n_face) float64 228MB dask.array<chunksize=(2035381,), meta=np.ndarray>
    Elevation              (n_face) float64 228MB dask.array<chunksize=(2035381,), meta=np.ndarray>
    Slope                  (n_face) float64 228MB dask.array<chunksize=(2035381,), meta=np.ndarray>
    Aspect                 (n_face) float64 228MB dask.array<chunksize=(2035381,), meta=np.ndarray>
    Area                   (n_face) float64 228MB dask.array<chunksize=(2035381,), meta=np.ndarray>
    owner                  (n_face) int32 114MB dask.array<chunksize=(4070761,), meta=np.ndarray>
# save just mesh
# if this isn't done it can lead to file lock issues later
ds.chm.uxgrid_to_netcdf('mesh.nc')
del ds

# chunk depending on your data set
ds = ux.open_mfdataset('mesh.nc', 'output_2016.nc', chunks={'time':1})
# global id needed for output
# use the Slice to preserve the time dim
ds = ds[['t','global_id']].isel(time=slice(0,2)).chm.clip(shp_file_path='buffered_FABDEM-clip.shp')
ds.chm.uxgrid_to_netcdf('mesh_subset.nc')
ds.chm.vars_to_netcdf('vars_subset.nc')
(
    ds.t.isel(time=0).plot() + ds.t.isel(time=1).plot() 
).cols(1)