Feature/cartesian spatial hash - Add support for computing barycentric coordinates from Cartesian coordinates

[fluidnumericsJoe]

Closes#1187

Overview

This PR provides an option to set the coordinate_system for the SpatialHash class. This is specifically used to control how the barycentric coordinates are calculated. When set to cartesian we use the cartesian coordinates to compute the barycentric coordinates.

This algorithm assumes faces are planar (linear), which is equivalent to using a local linearization of the spherical coordinates transformation to compute the barycentric coordinates. Determining if a point is in a face requires (1) valid barycentric coordinates and (2) the point is in the plane spanned by the planar face. Meeting these two criteria reduces to ensuring that the projection error of a query point onto the planar face is bounded by the local spherical coordinate linearization error. Because of this, the SpatialHash object now has a new attribute _ll_tolerance that is initialized in the private SpatialHash._local_linearization_tolerance method.

Additional private methods include

• uxarray.grid.neighbors._barycentric_coordinates_cartesian
• uxarray.grid.neighbors._triangle_area_cartesian
• uxarray.grid.neighbors._local_projection_error

For end users :

• SpatialHash.coordinate_system = 'cartesian' should be used when grids are global or cross the antimeridian.

DRAFT WORK IN PROGRESS

Expected Usage

PR Checklist

General

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

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

Examples

• Any new notebook examples added to docs/examples/ folder
• Clear the output of all cells before committing
• New notebook files added to docs/examples.rst toctree
• New notebook files added to new entry in docs/gallery.yml with appropriate thumbnail photo in docs/_static/thumbnails/

[fluidnumericsJoe]

@philipc2 - I recall you mentioning in an issue (that I can't seem to find) about problems along the antimeridian with spatial hashing. I think this would be the PR to get such a test in to verify that the Cartesian coordinates for spatial hashing resolves that problem.

[fluidnumericsJoe]

@philipc2 - In this latest version, I've added a test that uses a global grid and runs a few spatial hash queries on either side of the antimeridian. I can get this to pass just fine, after a few adjustments to the hash table setup and when using the barycentric coordinates calculated with Cartesian coordinates.

However, the changes to the hash table, which relate to the treatment face_bounds_lon and the resulting indexing through the hash table seem to be problematic. The simplest test case to consider is test_is_inside.

deftest_is_inside():
"""Verifies simple test for points inside and outside an element."""verts= [(0.0, 90.0), (-180, 0.0), (0.0, -90)]
uxgrid=ux.open_grid(verts, latlon=True)
# Verify that a point outside the element returns a face id of -1face_ids, bcoords=uxgrid.get_spatial_hash().query([90.0, 0.0])
assertface_ids[0] ==-1# Verify that a point inside the element returns a face id of 0face_ids, bcoords=uxgrid.get_spatial_hash().query([-90.0, 0.0])
assertface_ids[0] ==0assertnp.allclose(bcoords[0], [0.25, 0.5, 0.25], atol=1e-06)

In this test, there is a single face with vertices that do not cross the antimeridian; the face abuts the antimeridian. However, in this case I find that uxgrid.antimeridian_face_indices = [0] and the uxgrid.face_lon_bounds = [0, -180.0] . This forces the hash table indexing to index from x=0 eastward and "wrap around" from the west edge of the hash grid and continue eastward to the hash cell nearest x=-180.0

Previously, I had been sorting the face_lon_bounds which hid this issue. However, it causes problems with other grids in capturing faces that straddle the antimeridian.

There's another test that uses this same single face grid that is failing for the same reason at the moment. Last, the test involving the mpas dual grid is also failing for a similar reason (the face_lon_bounds are returning in the wrong order and cells are incorrectly marked as antimeridian).

I'm posting this in the hopes that I'm misunderstanding the output of face_lon_bounds and the antimeridian_face_indices but I've been stuck on this bit today. However, once this is resolved, I should be good to start modifying documentation and converting this draft to a PR.

Let me know if hopping on a screen share would be best to discuss

[fluidnumericsJoe]

Changing the vertex ordering so that vertices are listed counterclockwise in the test_is_inside test case now passes.

verts = [(-180, 0.0), (0.0, -90.0), (0.0, 90.0)]

[fluidnumericsJoe]

Now, only the mpas_dual test is failing. Is it possible that the dual grid creation is orienting vertices in the wrong direction (clockwise instead of counter-clockwise)??

[aaronzedwick]

Now, only the mpas_dual test is failing. Is it possible that the dual grid creation is orienting vertices in the wrong direction (clockwise instead of counter-clockwise)??

Yes, that is very possible. I think that is likely the issue.

[philipc2]

This is a good catch @fluidnumerics-joe

I've had some issues with those MPAS grids before, it does appear to be an issue with the Grid and not your functionality.

[philipc2]

I'm also curious if you've had a chance to try out our recently optimized Grid.get_faces_containing_point() functionality. I haven't been following to closely on the spatial hashing work, but I'm curious how the results and perfromance differs for your application.

We've recently updated it to allow for a batch of points to be submitted, as opposed to needing to iterate for each point.

[fluidnumericsJoe]

I'm also curious if you've had a chance to try out our recently optimized Grid.get_faces_containing_point() functionality. I haven't been following to closely on the spatial hashing work, but I'm curious how the results and perfromance differs for your application.

We've recently updated it to allow for a batch of points to be submitted, as opposed to needing to iterate for each point.

Oh man! I'm gonna have to try this out in Parcels. I'll put together a few simple tests once we're on the other side of this PR and get back to you with results :)

[fluidnumericsJoe]

Would you want me to attempt to fix the issue in this PR or should we resolve that separately ? I suppose we'd need a different test to check orientation of vertices to confirm if this is indeed the issue..

[aaronzedwick]

Would you want me to attempt to fix the issue in this PR or should we resolve that separately ? I suppose we'd need a different test to check orientation of vertices to confirm if this is indeed the issue..

I think a separate issue would be most appropriate, since it isn’t a problem with anything you are making here. I think raising an issue with an example and tackling it separately would work best.

[fluidnumericsJoe]

@aaronzedwick - Got it. I'll see if I can come up with something that doesn't involve the spatial hashing where this might also crop up.

[fluidnumericsJoe]

I did some testing and found that the failures in spatial hash query on the MPAS dual grid failures are not related to vertex ordering within elements.

I did the following to ensure that the vertices were ordered in the correct way

# Compute signed area for each face using shoelace formula
# This will help determine the orientation of the face (clockwise or counter-clockwise)
def signed_area(x, y):
return 0.5 * np.sum(x[:-1] * y[1:] - x[1:] * y[:-1])
n_faces = face_nodes.shape[0]
new_face_nodes = face_nodes.copy()
for i in range(n_faces):
node_ids = face_nodes[i, :]
node_ids = node_ids[~np.isnan(node_ids)].astype(int)
x = vertex_coords[0][node_ids]
y = vertex_coords[1][node_ids]
# Close the loop
x = np.append(x, x[0])
y = np.append(y, y[0])
area = signed_area(x, y)
# If area is negative, the face is CW -> reverse
if area < 0:
print(f"Face {i} is CW, reversing node order")
new_face_nodes[i, :len(node_ids)] = node_ids[::-1]
# Overwrite the mesh connectivity with corrected ordering
uxg["face_node_connectivity"].values[:] = new_face_nodes.astype(int)

While it did find that some elements indeed were clockwise and others were counter-clockwise, reversing the order of vertices did not resolve the issue. Digging into this a bit more

[philipc2]

Hi @fluidnumerics-joe

Wanted to follow up and see if there's anything you need from us.

If you think we should instead focus on the get_faces_containing_point() functionality (as you've explored in Parcels-code/Parcels#2103, feel free to close this PR