Regrid to structured grid

[Chrismarsh]

I would like to regrid an unstructured mesh to a structured mesh for further processing to tiff, etc with xarray+rioxarray.

ds = ux.open_mfdataset("grid_kan.nc", 'run/output_2017_hrdps/output_2017_hrdps.nc')

xmin = ds.Mesh2_node_x.min().values
xmax = ds.Mesh2_node_x.max().values
ymin = ds.min().values
ymax = ds.Mesh2_node_y.max().values

x = np.abs(xmin-xmax)
y = np.abs(ymin-ymax)

dxdy = 0.01

numX, numY = int(x/dxdy), int(y/dxdy)

print(f'numX, numY = {numX}, {numY}')

print(f'xmin={xmin} xmax={xmax} ymin={ymin} ymax={ymax} ')

# cell centres
dxdy2 = dxdy/2.
x_center = np.linspace(start=xmin + dxdy2,
                        stop=xmax - dxdy2, num=numX)

y_center = np.linspace(start=ymin + dxdy2,
                        stop=ymax - dxdy2, num=numY)



target_grid = ux.Grid.from_structured(lat=y_center, lon=x_center)
target_grid.uxgrid

remap = ds.isel(time=0)["t"].remap.nearest_neighbor(
    destination_grid=target_grid.uxgrid, remap_to="faces"
)
remap

results in

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[13], [line 30](vscode-notebook-cell:?execution_count=13&line=30)
     24 y_center = np.linspace(start=ymin + dxdy2,
     25                         stop=ymax - dxdy2, num=numY)
     29 target_grid = ux.Grid.from_structured(lat=y_center, lon=x_center)
---> [30](vscode-notebook-cell:?execution_count=13&line=30) target_grid.uxgrid

AttributeError: 'Grid' object has no attribute 'uxgrid'

Is this not a supported use case?

[philipc2]

Since target_grid is already a Grid, you should do the following instead:

target_grid = ux.Grid.from_structured(lat=y_center, lon=x_center)

remap = ds.isel(time=0)["t"].remap.nearest_neighbor(
    destination_grid=target_grid, remap_to="faces"
)

[Chrismarsh]

That's what I had thought, but it doesn't work, so I was just following the docs

remap = ds.isel(time=0)["t"].remap.nearest_neighbor(
    destination_grid=target_grid, remap_to="faces"
)
remap

[...]

29 def _nearest_neighbor_query(
     30     source_grid: Grid,
     31     destination_grid: Grid,
   (...)
     35     return_distances: bool = False,
     36 ):
     37     """
     38     Query the nearest neighbors from a source grid for specified destination points.
     39 
   (...)
     63         Distances to the nearest source points, returned only if `return_distances` is True.
     64     """
---> [65](https://vscode-remote+ssh-002dremote-002bppp5-002dvscode.vscode-resource.vscode-cdn.net/fs/homeu2/eccc/crd/ords/ccrp/chm003/spack/var/spack/environments/analysis/.spack-env/view/lib/python3.12/site-packages/uxarray/remap/nearest_neighbor.py:65)     source_tree = source_grid._get_scipy_kd_tree(coordinates=KDTREE_DIM_MAP[source_dim])
     66     destination_points = _prepare_points(destination_grid, destination_dim)
     67     distances, nearest_indices = source_tree.query(destination_points, k=k, workers=-1)

AttributeError: 'NoneType' object has no attribute '_get_scipy_kd_tree'

[Chrismarsh]

Oh, it's because something is weird with the isel. ds is correctly defined

xarray.UxDataset> Size: 2TB
Dimensions:                             (n_face: 365193, Three: 3,
                                         n_node: 207826, time: 8737)
Coordinates:
  * time                                (time) datetime64[ns] 70kB 2017-10-01...
    Mesh2_face_x                        (n_face) float64 3MB dask.array<chunksize=(365193,), meta=np.ndarray>
    Mesh2_face_y                        (n_face) float64 3MB dask.array<chunksize=(365193,), meta=np.ndarray>
Dimensions without coordinates: n_face, Three, n_node
Data variables: (12/141)
    Mesh2                               int32 4B ...
    Mesh2_face_nodes                    (n_face, Three) uint32 4MB dask.array<chunksize=(365193, 3), meta=np.ndarray>
    Mesh2_node_x                        (n_node) float64 2MB dask.array<chunksize=(207826,), meta=np.ndarray>
    Mesh2_node_y                        (n_node) float64 2MB dask.array<chunksize=(207826,), meta=np.ndarray>
    Mesh2_node_z                        (n_node) float64 2MB dask.array<chunksize=(207826,), meta=np.ndarray>
    Mesh2_node_z_paraview               (n_node) float64 2MB dask.array<chunksize=(207826,), meta=np.ndarray>
    ...                                  ...
    param_Ninja7_V                      (n_face) float64 3MB dask.array<chunksize=(365193,), meta=np.ndarray>
    Elevation                           (n_face) float64 3MB dask.array<chunksize=(365193,), meta=np.ndarray>
    Slope                               (n_face) float64 3MB dask.array<chunksize=(365193,), meta=np.ndarray>
    Aspect                              (n_face) float64 3MB dask.array<chunksize=(365193,), meta=np.ndarray>
    Area                                (n_face) float64 3MB dask.array<chunksize=(365193,), meta=np.ndarray>

but once I do .isel(time=0) soemthing goes wrong. It can no longer be html printed in jupyter

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File ~/spack/var/spack/environments/analysis/.spack-env/view/lib/python3.12/site-packages/IPython/core/formatters.py:347, in BaseFormatter.__call__(self, obj)
    345     method = get_real_method(obj, self.print_method)
    346     if method is not None:
--> 347         return method()
    348     return None
    349 else:

File ~/spack/var/spack/environments/analysis/.spack-env/view/lib/python3.12/site-packages/uxarray/core/dataarray.py:89, in UxDataArray._repr_html_(self)
     87 if OPTIONS["display_style"] == "text":
     88     return f"<pre>{escape(repr(self))}</pre>"
---> 89 return array_repr(self)

File ~/spack/var/spack/environments/analysis/.spack-env/view/lib/python3.12/site-packages/uxarray/formatting_html.py:204, in array_repr(arr)
    200     sections.append(xrfm.index_section(indexes))
    202 sections.append(xrfm.attr_section(arr.attrs))
--> 204 return _obj_repr_with_grid(arr, header_components, sections)

File ~/spack/var/spack/environments/analysis/.spack-env/view/lib/python3.12/site-packages/uxarray/formatting_html.py:125, in _obj_repr_with_grid(obj, header_components, sections)
    122 header = f"<div class='xr-header'>{''.join(h for h in header_components)}</div>"
    123 sections = "".join(f"<li class='xr-section-item'>{s}</li>" for s in sections)
--> 125 grid_html_repr = grid_repr(
    126     obj.uxgrid,
    127     max_items_collapse=0,
    128     header_name=f"uxarray.{type(obj).__name__}.uxgrid",
    129 )
    131 icons_svg, css_style = xrfm._load_static_files()
    132 obj_repr_html = (
    133     "<div>"
    134     f"{icons_svg}<style>{css_style}</style>"
   (...)
    140     "</div>"
    141 )

File ~/spack/var/spack/environments/analysis/.spack-env/view/lib/python3.12/site-packages/uxarray/formatting_html.py:79, in grid_repr(grid, max_items_collapse, header_name)
     76 """HTML repr for ``Grid`` class."""
     77 header_components = _grid_header(grid, header_name)
---> 79 sections = _grid_sections(grid, max_items_collapse)
     81 return xrfm._obj_repr(grid, header_components, sections)

File ~/spack/var/spack/environments/analysis/.spack-env/view/lib/python3.12/site-packages/uxarray/formatting_html.py:23, in _grid_sections(grid, max_items_collapse)
     21 def _grid_sections(grid, max_items_collapse=15):
     22     cartesian_coordinates = list(
---> 23         [coord for coord in ugrid.CARTESIAN_COORDS if coord in grid._ds]
     24     )
     25     spherical_coordinates = list(
     26         [coord for coord in ugrid.SPHERICAL_COORDS if coord in grid._ds]
     27     )
     28     descriptor = list(
     29         [desc for desc in descriptors.DESCRIPTOR_NAMES if desc in grid._ds]
     30     )

AttributeError: 'NoneType' object has no attribute '_ds'

[Chrismarsh]

Ahhh ok @philipc2, one problem is that this

ds['swe'].isel(time=0)

works, but

ds.isel(time=0)['swe']

does not. But I would expect that to work based on my xarray usage. Subseting a dataset in time and then by variable or by variable and then in time, should be equivalent reductions.

Having fixed that I now get

ValueError: conflicting sizes for dimension 'n_face': length 2684 on 'swe' and length 365193 on {'n_face': 'Mesh2_face_x'}

when I call the regrid:

remap = ds['swe'].isel(time=0).remap.nearest_neighbor(
    destination_grid=target_grid, remap_to="faces")
remap

[Chrismarsh]

Further to this, remap_to="nodes" works

[Chrismarsh]

if I use .topological_mean(destination="face") to go from the remapped nodes to faces, it then "works". Feels a bit hacky. I remain unsure what I'm doing wrong.