Regrid

The unstructured mesh is optimal for computations and storage because it reduces the number of computational elements. However, using structured netcdf, zarr, or tiff can be preferable for some analysis and plot making.

import CHM as pc
import uxarray as ux
import rioxarray as rio
import xarray as xr
import matplotlib.pyplot as plt

# If the mesh requires out of core compute, spin up a dask cluster

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

# Here, we are using a single timestep for a subset domain for air temp
# note that global_id is preserved as it must be in both the mesh and vars .nc in order to line up elements
ds = ux.open_mfdataset('mesh_subset.nc', 'vars_subset.nc')
ds

<xarray.UxDataset> Size: 5MB
Dimensions:       (n_face: 142546)
Coordinates:
    time          datetime64[ns] 8B ...
    Mesh2_face_x  (n_face) float64 1MB dask.array<chunksize=(142546,), meta=np.ndarray>
    Mesh2_face_y  (n_face) float64 1MB dask.array<chunksize=(142546,), meta=np.ndarray>
Dimensions without coordinates: n_face
Data variables:
    t             (n_face) float64 1MB dask.array<chunksize=(142546,), meta=np.ndarray>
    global_id     (n_face) uint64 1MB dask.array<chunksize=(142546,), meta=np.ndarray>

<xarray.UxDataset> Size: 5MB
Dimensions:       (n_face: 142546)
Coordinates:
    time          datetime64[ns] 8B ...
    Mesh2_face_x  (n_face) float64 1MB dask.array<chunksize=(142546,), meta=np.ndarray>
    Mesh2_face_y  (n_face) float64 1MB dask.array<chunksize=(142546,), meta=np.ndarray>
Dimensions without coordinates: n_face
Data variables:
    t             (n_face) float64 1MB dask.array<chunksize=(142546,), meta=np.ndarray>
    global_id     (n_face) uint64 1MB dask.array<chunksize=(142546,), meta=np.ndarray>

uxarray.UxDataset

• Dimensions:
  • n_face: 142546
• Coordinates: (3)
  • time
    ()
    datetime64[ns]
    ...

    standard_name :

    time

    long_name :

    Time

    [1 values with dtype=datetime64[ns]]

  • Mesh2_face_x
    (n_face)
    float64
    dask.array<chunksize=(142546,), meta=np.ndarray>

    standard_name :

    latitude

    long_name :

    Characteristics latitude of 2D mesh triangle (e.g. circumcenter coordinate).

    units :

    degrees_north

    Array Chunk Bytes 1.09 MiB 1.09 MiB Shape (142546,) (142546,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • Mesh2_face_y
    (n_face)
    float64
    dask.array<chunksize=(142546,), meta=np.ndarray>

    standard_name :

    latitude

    long_name :

    Characteristics latitude of 2D mesh triangle (e.g. circumcenter coordinate).

    units :

    degrees_north

    Array Chunk Bytes 1.09 MiB 1.09 MiB Shape (142546,) (142546,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

• Data variables: (2)
  • t
    (n_face)
    float64
    dask.array<chunksize=(142546,), meta=np.ndarray>

    _QuantizeBitGroomNumberOfSignificantDigits :

    4

    mesh :

    Mesh2

    location :

    face

    Array Chunk Bytes 1.09 MiB 1.09 MiB Shape (142546,) (142546,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • global_id
    (n_face)
    uint64
    dask.array<chunksize=(142546,), meta=np.ndarray>

    mesh :

    Mesh2

    location :

    face

    Array Chunk Bytes 1.09 MiB 1.09 MiB Shape (142546,) (142546,) Dask graph 1 chunks in 2 graph layers Data type uint64 numpy.ndarray

• Indexes: (0)
• Attributes: (0)

Show Grid Information

<uxarray.Grid>
Original Grid Type: UGRID
Grid Dimensions:
  * n_node: 72536
  * n_face: 142546
  * n_max_face_nodes: 3
  * n_max_node_faces: 11
Grid Coordinates (Spherical):
  * node_lon: (72536,)
  * node_lat: (72536,)
Grid Coordinates (Cartesian):
  * node_x: (72536,)
  * node_y: (72536,)
  * node_z: (72536,)
Grid Connectivity Variables:
  * node_face_connectivity: (72536, 11)
  * face_node_connectivity: (142546, 3)
Grid Descriptor Variables:
  * n_nodes_per_face: (142546,)

uxarray.UxDataset.uxgrid

• Dimensions:
  • n_face: 142546
  • n_max_face_nodes: 3
  • n_node: 72536
  • n_max_node_faces: 11
• Spherical Coordinates: (2)
  • node_lon
    (n_node)
    float64
    -115.0 -115.0 ... -115.4 -115.4

    standard_name :

    longitude

    long_name :

    Longitude of 2D mesh nodes.

    units :

    degrees_east

    array([-115.033806, -114.963399, -115.079893, ..., -115.419248, -115.417941,
           -115.419742], shape=(72536,))

  • node_lat
    (n_node)
    float64
    ...

    standard_name :

    latitude

    long_name :

    Latitude of 2D mesh nodes.

    units :

    degrees_north

    [72536 values with dtype=float64]

• Cartesian Coordinates: (3)
  • node_x
    (n_node)
    float64
    ...

    standard_name :

    x

    long name :

    Cartesian x location of the corner nodes of each face

    units :

    meters

    [72536 values with dtype=float64]

  • node_y
    (n_node)
    float64
    ...

    standard_name :

    y

    long name :

    Cartesian y location of the corner nodes of each face

    units :

    meters

    [72536 values with dtype=float64]

  • node_z
    (n_node)
    float64
    ...

    standard_name :

    z

    long name :

    Cartesian z location of the corner nodes of each face

    units :

    meters

    [72536 values with dtype=float64]

• Connectivity: (2)
  • node_face_connectivity
    (n_node, n_max_node_faces)
    int64
    343 1453 ... -9223372036854775808

    cf_role :

    node_face_connectivity

    long name :

    Faces that neighbor each node

    start_index :

    0

    _FillValue :

    -9223372036854775808

    array([[                 343,                 1453,                 1542,
            ..., -9223372036854775808, -9223372036854775808,
            -9223372036854775808],
           [                 827,                  831,                  898,
            ..., -9223372036854775808, -9223372036854775808,
            -9223372036854775808],
           [                1428,                 2622,                 2625,
            ..., -9223372036854775808, -9223372036854775808,
            -9223372036854775808],
           ...,
           [              141526,               142480,               142521,
            ..., -9223372036854775808, -9223372036854775808,
            -9223372036854775808],
           [              136501,               137622,               141534,
            ..., -9223372036854775808, -9223372036854775808,
            -9223372036854775808],
           [              142473,               142492,               142522,
            ..., -9223372036854775808, -9223372036854775808,
            -9223372036854775808]], shape=(72536, 11))

  • face_node_connectivity
    (n_face, n_max_face_nodes)
    int64
    765 1357 734 ... 72431 71088 72517

    cf_role :

    face_node_connectivity

    long_name :

    Maps every triangular face to its three corner nodes.

    _FillValue :

    -9223372036854775808

    array([[  765,  1357,   734],
           [  760,   822,   966],
           [ 1218,  1225,   155],
           ...,
           [72024, 72022, 72023],
           [72509, 72521, 72535],
           [72431, 71088, 72517]], shape=(142546, 3))

• Descriptors: (1)
  • n_nodes_per_face
    (n_face)
    int64
    ...

    cf_role :

    n_nodes_per_face

    long name :

    Number of nodes per face

    [142546 values with dtype=int64]

• Attributes: (1)

  source_grid_spec :

  UGRID

pyCHM is built on top of Uxarray, and thus can use the Uxarray regridding for small domains. Small is defined as anything that can be gridded on a single compute node as the underlying regridding is cython and not parallel across dask https://uxarray.readthedocs.io/en/latest/user-guide/remapping.html Which the uxarray regridder can be used directly, the .chm.regrid version handles the grid creation and other minutia.

# Regrid using the uxarray regridding
rg_01 = ds.chm.regrid(0.01) # in degrees
rg_001 = ds.chm.regrid(0.001)

fig, ax = plt.subplots(nrows=2,figsize=(8,10))
rg_01.t.plot(ax=ax[0], robust=True)
rg_001.t.plot(ax=ax[1],  robust=True)

<matplotlib.collections.QuadMesh at 0x1521d39f4920>

These can then me saved to a gridded format of choice that xarray can write

# GeoTiff
tiff = rg_01.rio.write_crs('+proj=longlat +datum=WGS84 +no_defs +type=crs')
tiff.t.rio.to_raster('test.tiff')

tiff = rio.open_rasterio('test.tiff')
tiff.plot()

<matplotlib.collections.QuadMesh at 0x1521ff86fb90>

However, with large meshes, using the chm.regrid will run out of memory. Then it’s time for the chm_ugrid2grid, which is an MPI backed regridder built on top of ESMF. It supports conservative and bilinear regridding. This tool doesn’t use Dask, so needs to be invoked with mpirun or similar. It writes to zarr or tiff.

%%bash
# the variable.nc MUST have global_id var and time dimension
# the mesh.nc MUST have global_id var
mpirun -n 5 chm_ugrid2grid vars_subset.nc --mesh mesh_subset.nc --dxdy 0.001 --method conservative --variables t --zarr-output subset.zarr --overwrite

[0] Initializing Zarr store at subset.zarr
[0] 20161001T000000 - t
[2] 20161001T000000 - t
[3] 20161001T000000 - t
[4] 20161001T000000 - t
[1] 20161001T000000 - t
[4] 20161001T010000 - t
[3] 20161001T010000 - t
[2] 20161001T010000 - t
[0] 20161001T010000 - t
[1] 20161001T010000 - t

%%bash
# If invariant with time variables like Elevation, Slope are used then it is important to set
# --timeoffset 0
# to obtain only 1 output of that variable
# for example"
mpirun -n 5 chm_ugrid2grid vars_subset.nc --mesh mesh_subset.nc --dxdy 0.001 --method conservative --variables Elevation Slope Aspect --zarr-output invariants.zarr --overwrite --timeoffset 0

ds = xr.open_zarr('subset.zarr', consolidated=False) # because of how the zarr is written in parallel it doesn't have consolidated metadata
ds

<xarray.Dataset> Size: 2MB
Dimensions:      (time: 2, latitude: 448, longitude: 613)
Coordinates:
  * time         (time) datetime64[ns] 16B 2016-10-01 2016-10-01T01:00:00
  * longitude    (longitude) float64 5kB -115.6 -115.6 -115.6 ... -115.0 -115.0
  * latitude     (latitude) float64 4kB 50.62 50.62 50.62 ... 51.06 51.06 51.06
Data variables:
    t            (time, latitude, longitude) float32 2MB dask.array<chunksize=(1, 448, 512), meta=np.ndarray>
    spatial_ref  int8 1B ...
Attributes:
    Conventions:     CF-1.8, GeoZarr-1.0
    geospatial_crs:  EPSG:4326

xarray.Dataset

• Dimensions:
  • time: 2
  • latitude: 448
  • longitude: 613
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2016-10-01 2016-10-01T01:00:00

    _ARRAY_DIMENSIONS :

    ['time']

    array(['2016-10-01T00:00:00.000000000', '2016-10-01T01:00:00.000000000'],
          dtype='datetime64[ns]')

  • longitude
    (longitude)
    float64
    -115.6 -115.6 ... -115.0 -115.0

    _ARRAY_DIMENSIONS :

    ['longitude']

    array([-115.566438, -115.565437, -115.564437, ..., -114.95598 , -114.954979,
           -114.953978], shape=(613,))

  • latitude
    (latitude)
    float64
    50.62 50.62 50.62 ... 51.06 51.06

    _ARRAY_DIMENSIONS :

    ['latitude']

    array([50.615141, 50.616143, 50.617144, ..., 51.060925, 51.061926, 51.062928],
          shape=(448,))

• Data variables: (2)
  • t
    (time, latitude, longitude)
    float32
    dask.array<chunksize=(1, 448, 512), meta=np.ndarray>

    _ARRAY_DIMENSIONS :

    ['time', 'latitude', 'longitude']

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 2.10 MiB 896.00 kiB Shape (2, 448, 613) (1, 448, 512) Dask graph 4 chunks in 2 graph layers Data type float32 numpy.ndarray

  • spatial_ref
    ()
    int8
    ...

    _ARRAY_DIMENSIONS :

    []

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    grid_mapping_name :

    latitude_longitude

    epsg_code :

    EPSG:4326

    semi_major_axis :

    6378137.0

    inverse_flattening :

    298.257223563

    longitude_of_prime_meridian :

    0.0

    latitude_of_prime_meridian :

    0.0

    [1 values with dtype=int8]

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2016-10-01 00:00:00', '2016-10-01 01:00:00'], dtype='datetime64[ns]', name='time', freq=None))

  • longitude
    PandasIndex

    PandasIndex(Index([-115.56643805770915, -115.56543730624433, -115.56443655477953,
           -115.56343580331472,  -115.5624350518499, -115.56143430038509,
           -115.56043354892029, -115.55943279745547, -115.55843204599066,
           -115.55743129452584,
           ...
           -114.96298492442777, -114.96198417296296, -114.96098342149814,
           -114.95998267003333, -114.95898191856853, -114.95798116710371,
            -114.9569804156389, -114.95597966417408, -114.95497891270928,
           -114.95397816124446],
          dtype='float64', name='longitude', length=613))

  • latitude
    PandasIndex

    PandasIndex(Index([ 50.61514095548525,  50.61614271665699,  50.61714447782873,
            50.61814623900047,   50.6191480001722, 50.620149761343946,
            50.62115152251568,  50.62215328368742,  50.62315504485915,
           50.624156806030896,
           ...
            51.05391234870652,  51.05491410987826,     51.05591587105,
            51.05691763222173,  51.05791939339347,  51.05892115456521,
            51.05992291573695, 51.060924676908684,  51.06192643808043,
            51.06292819925216],
          dtype='float64', name='latitude', length=448))

• Attributes: (2)

  Conventions :

  CF-1.8, GeoZarr-1.0

  geospatial_crs :

  EPSG:4326

fig, ax = plt.subplots(nrows=2,figsize=(8,10))
ds.t.isel(time=0).plot(ax=ax[0], robust=True)
ds.t.isel(time=1).plot(ax=ax[1],  robust=True)

<matplotlib.collections.QuadMesh at 0x1521fface3c0>