Dual Mesh Construction

[aaronzedwick]

Closes #825

Overview

Constructs the dual mesh of a grid using face_centers and node_face_connectivity

Expected Usage

import uxarray as ux

grid_path = "/path/to/grid.nc"

grid = ux.open_grid(grid_path)

dual = grid.compute_dual()

PR Checklist

General

• An issue is linked created and linked
• Add appropriate labels
• Filled out Overview and Expected Usage (if applicable) sections

Testing

• Adequate tests are created if there is new functionality
• Tests cover all possible logical paths in your function
• Tests are not too basic (such as simply calling a function and nothing else)

Documentation

• Docstrings have been added to all new functions
• Docstrings have updated with any function changes
• Internal functions have a preceding underscore (_) and have been added to docs/internal_api/index.rst
• User functions have been added to docs/user_api/index.rst

[aaronzedwick]

@philipc2 Thoughts on the API? Currently, I return a dual grid object as its own thing. This made sense to me, as changing the original grid would be difficult due to the different dimensions.

[philipc2]

@philipc2 Thoughts on the API? Currently, I return a dual grid object as its own thing. This made sense to me, as changing the original grid would be difficult due to the different dimensions.

This type of implementation is what I was envisioning, since we don't want anything to be replaced in-place when creating the dual grid.

Let's name it something like Grid.compute_dual(method=...), where we method could indicate what type of dual we'd like to compute.

[philipc2]

We should add functionality that can also convert a UxDataset or UxDataArray into it's "dual" form, which would involve making any face-centered data variable now node-centered, and vice-versa. This would however cause issues if the data variable is edge centered, since those edges would no longer exist.

[philipc2]

This was the block that we were able to get working in the Notebook.

import uxarray as ux
from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
import numpy as np

uxgrid = ux.open_grid("/Users/philipc/uxarray-gold/uxarray-gold/uxarray/test/meshfiles/geos-cs/c12/test-c12.native.nc4")


lonlat_t = [(lon, lat) for lon, lat in zip(uxgrid.node_lon.values, uxgrid.node_lat.values)]

# # Dictionary to track first occurrence and subsequent indices
occurrences = {}

# Iterate through the list and track occurrences
for index, tpl in enumerate(lonlat_t):
    if tpl in occurrences:
        occurrences[tpl].append((INT_DTYPE(index)))
    else:
        occurrences[tpl] = [INT_DTYPE(index)]
duplicate_dict = {}

for tpl, indices in occurrences.items():
    if len(indices) > 1:
        source_idx = indices[0]
        for duplicate_idx in indices[1:]:
            duplicate_dict[duplicate_idx] = source_idx
            
new_face_node_connectivity = uxgrid.face_node_connectivity.values.copy().ravel()

for idx, item in enumerate(new_face_node_connectivity):
    # O(1)
    if item in duplicate_dict:
        new_face_node_connectivity[idx] = duplicate_dict[item]


new_face_node_connectivity = new_face_node_connectivity.reshape((uxgrid.n_face, uxgrid.n_max_face_nodes))

node_face_conn = {node_i: [] for node_i in range(uxgrid.n_node)}
for face_i, face_nodes in enumerate(new_face_node_connectivity):
    for node_i in face_nodes:
        if node_i != ux.INT_FILL_VALUE:
            node_face_conn[node_i].append(face_i)

n_max_node_faces = -1   
for face_indicies in node_face_conn.values():
    if len(face_indicies) > n_max_node_faces:
        n_max_node_faces = len(face_indicies)
                     

node_face_connectivity = np.full((uxgrid.n_node, n_max_node_faces), INT_FILL_VALUE)

for node_idx, face_indices in enumerate(node_face_conn.values()):
    n_faces = len(face_indices)
    node_face_connectivity[node_idx, 0:n_faces] = face_indices
    
new_uxgrid = ux.Grid.from_topology(uxgrid.node_lon.values, 
                                   uxgrid.node_lat.values, 
                                   new_face_node_connectivity,
                                   node_face_connectivity=node_face_connectivity)

[aaronzedwick]

Thanks @philipc2! I will see if I can help figure out what's going on.

[aaronzedwick]

It is definitely something with how you are updating and returning the grid, I did a very basic implementation just as a function inside the grid class, called by compute dual at the start, and it works perfectly. Red is the original grid.

[rajeeja]

It is definitely something with how you are updating and returning the grid, I did a very basic implementation just as a function inside the grid class, called by compute dual at the start, and it works perfectly. Red is the original grid.

Right, it'd be better to compare the grid objects returned in each case, specifically the merged nodes part in both cases (notebook and installed versions)

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[philipc2]

It is definitely something with how you are updating and returning the grid, I did a very basic implementation just as a function inside the grid class, called by compute dual at the start, and it works perfectly. Red is the original grid.

Yeah, I'm almost wondering if there might be a bug somewhere in our Grid.from_topology(). Going to do some more digging.

[aaronzedwick]

It is definitely something with how you are updating and returning the grid, I did a very basic implementation just as a function inside the grid class, called by compute dual at the start, and it works perfectly. Red is the original grid.

Yeah, I'm almost wondering if there might be a bug somewhere in our Grid.from_topology(). Going to do some more digging.

It's working now though? I think it was just a problem with pip install . I have been having trouble getting it to work. That's what the notebook is for, just testing.

[philipc2]

Using UXarray we can construct the dual mesh using grid.compute_dual() , which returns a new grid object.

Reply via ReviewNB

[aaronzedwick]

You want me to replace everything with just this?

[philipc2]

Whoops, meant only the final portion.

[philipc2]

Can you use the .relabel() method with a legend for this instead?

Below is an example:

https://uxarray--859.org.readthedocs.build/en/859/user-guide/topological-aggregations.html#:~:text=%22edge%22)-,(%0A%20%20%20%20uxds.uxgrid.plot(line_color%3D%22black%22,500%2C%20title%3D%22Node%20to%20Edge%20Aggregation%20(Mean)%22),-Edge%20Aggregations

Reply via ReviewNB

[aaronzedwick]

.relabel() doesn't seem to work when combining plots of lines. It works in your notebook because everything you combined are points. Do you know how to make it work for actual grid plots? @philipc2

[rajeeja]

is this good to go?