Point In Face & Get Faces Containing Point #1056

aaronzedwick · 2024-11-05T20:59:44Z

Overview

Expected Usage

from uxarray.grid.geometry import point_in_polygon

# Defined polygon
polygon = [ [-10,  10, -10, 10], [10, 10, -10, -10]]

# Point to check
point = [10, 10]

point_in_polygon(polygon, point, inclusive=True)

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

…/uxarray into zedwick/point-in-polygon

aaronzedwick · 2024-11-27T15:07:32Z

@ philipc2 do you think this should be an internal function or exposed to the user?

philipc2

Please include an ASV benchmark. I'd suggest doing a parameterized benchmark for the 120 and 480 km MPAS grids.

philipc2

Please use the optimized functions from #1072 and try to write the function entirely in Numba. This may require us to pass in both the cartesian and spherical versions of point & polygon. Let me know if you have any questions!

uxarray/grid/geometry.py

uxarray/grid/grid.py

rajeeja

point_in_face(faces_edge_cartesian[0], point, inclusive=True)

this always requires calculation of edges, can we have another version of this -

point_in_face(face_id, point or
point_in_face(face_conn, point)

For the notebook example, we have 4 faces with id's say 0 to 4 and connectivity as shown below:


xarray.DataArray'face_node_connectivity'n_face: 4n_max_face_nodes: 6

array([[ 0,  1,  2,  3,  4,  5],
       [15,  7,  6, 12,  0,  5],
       [ 0, 12,  8,  9, 13,  1],
       [ 4, 14, 11, 10, 15,  5]])

getting edges adds another step that might not always be necessary and can be done internally for these other signatures, the two functions above would be more user friendly IMO.

aaronzedwick · 2025-02-04T00:07:23Z

point_in_face(faces_edge_cartesian[0], point, inclusive=True)

this always requires calculation of edges, can we have another version of this -

point_in_face(face_id, point or

point_in_face(face_conn, point)

For the notebook example, we have 4 faces with id's say 0 to 4 and connectivity as shown below:
xarray.DataArray'face_node_connectivity'n_face: 4n_max_face_nodes: 6

array([[ 0,  1,  2,  3,  4,  5],
       [15,  7,  6, 12,  0,  5],
       [ 0, 12,  8,  9, 13,  1],
       [ 4, 14, 11, 10, 15,  5]])
getting edges adds another step that might not always be necessary and can be done internally for these other signatures, the two functions above would be more user friendly IMO.

This function is not a grid function, so passing just the face connectivity isn’t enough. You have to pass the node coordinates of the face. Either way, you have to parse the coordinates of the face you want to use.

We can pass the node coordinates as nodes into it if you like, but I originally went with this implementation for two reasons. First, the conversion will in theory impact performance. And second, if the user has to construct the node coordinates and send them in anyway, it would require less effort to use the edge coordinates Cartesian function to create the array for the user. It is easier and ensures the correct array shape is being passed to the function. Also, that is the shape of array the ‘faces_containing_point’ use, as well as all of Hongyu’s code.

I definitely understand wanting to make them more user friendly, but I see ‘get_faces_containing_point’ as the more common user function. And as the point in face function is meant to be called without being attached to the grid, you will always have to construct some kind of coordinate information and pass it in, whether it is the nodes that make up faces, or the nodes that make up edges. Otherwise the function won’t have access to that information.

rajeeja

bug fix for conversion _lonlat_rad_to_xyz in coordinates.py, max_face_radius and get_faces_containing_points are critical changes. good work

philipc2

Please fix the ASV benchmark, it is still failing on our CI.

benchmarks/mpas_ocean.py

uxarray/grid/geometry.py

philipc2 · 2025-02-04T20:06:25Z

Some perfromance issues, but they appear to be due to our Grid.isel(), which I am already working on optimizations for.

This is 10,000 queries on a 30km MPAS grid

rajeeja · 2025-02-04T20:08:18Z

Some perfromance issues, but they appear to be due to our Grid.isel(), which I am already working on optimizations for.

This is 10,000 queries on a 30km MPAS grid

do you generate this locally or asv generates it in some automated fashion?

philipc2 · 2025-02-04T20:11:33Z

Some perfromance issues, but they appear to be due to our Grid.isel(), which I am already working on optimizations for.
This is 10,000 queries on a 30km MPAS grid

do you generate this locally or asv generates it in some automated fashion?

import uxarray as ux
import numpy as np
import numba as nb

print(nb.__version__)
import cProfile

profiler = cProfile.Profile()


grid_path = "/Users/philipc/PycharmProjects/uxarray/unstructured-grid-viz-cookbook/meshfiles/x1.655362.grid.nc"

uxgrid=ux.open_grid(grid_path)

face_lon = uxgrid.face_lon.values.astype(np.float64)
face_lat = uxgrid.face_lat.values.astype(np.float64)



profiler.enable()
# Profile something here
for i in range(10000):
    point = np.array([face_lon[i], face_lat[i]])
    res = uxgrid.get_faces_containing_point(point)
profiler.disable()

# Save stats
profiler.dump_stats('point_in_polygon_10000.prof')

Then I ran snakeviz point_in_polygon_10000.prof

philipc2 · 2025-02-04T20:16:41Z

uxarray/grid/grid.py

+        if len(point) == 2:
+            point_xyz = np.array(_lonlat_rad_to_xyz(*np.deg2rad(point)))
+            point_lonlat = point
+        elif len(point) == 3:
+            point_xyz = point
+            point_lonlat = np.array(_xyz_to_lonlat_deg(*point_xyz), dtype=np.float64)
+        else:
+            raise ValueError(
+                "Point must either be in spherical or cartesian coordinates."
+            )


I ran into an issue where if my points were float32 the Numba functionality would fail, since it would be trying to call Numba functions between float32 and float64 variables, which isn't possible.

Maybe doing a np.asaray(point, dtype=np.float64) could fix it. May you add a test for this to replicate the functionality?

In some of my tests I have found that using lonlat with a dtype=np.float64 can cause issues on the nodes, it will not find all the faces nearby. I think this is possible due to rounding precission issues maybe?

Casting to a higher precision shouldn't usually lead to issues.

It does, doing dtype=np.float32 ensures that Hongyus algorithm recognizes that it is on the great circle arc. But when it is float64 it doesn't seem to realize it is on the arc.

Do you have a minimal example? That's interesting.

When I query with Cartesian coordinates in f64, I get the correct indices.

point = np.array([grid.node_x[0].values, grid.node_y[0].values, grid.node_z[0].values], dtype=np.float64)

Though, with f32 I get a Numba error.

point = np.array([grid.node_x[0].values, grid.node_y[0].values, grid.node_z[0].values], dtype=np.float32)

No implementation of function Function(<built-in function dot>) found for signature: >>> dot(array(float64, 1d, C), array(float32, 1d, C))

This does appear to be a precision issue when using spherical coordinates, since we internally convert them to cartesian. We should try to avoid as many of these conversions as possible. @hongyuchen1030 and I have encountered this a lot in our previous implementations.

The error you are getting happens within Hongyus code. I have had it error out sometimes when using cartesian coordinates that aren't in float64

This could be fixed by ensuring that the only xyz coordinates are f64. But that still won't fix the precision issue of the spherical coordinates when they are changed to f64.

I wonder if adding a small eps to increase the radius of each bounding circle might help?

EDIT: I see you did that below.

philipc2

Can we include a test for points along a line of constant latitude?

    def test_point_along_arc():
        node_lon = np.array([-40, -40, 40, 40])
        node_lat = np.array([-20, 20, 20, -20])
        face_node_connectivity = np.array([[0, 1, 2, 3]], dtype=np.int64)

        uxgrid = ux.Grid.from_topology(node_lon, node_lat, face_node_connectivity)

        # point at exactly 20 degrees latitude
        out1 = uxgrid.get_faces_containing_point((0, 20))

        # point at 25.41 degrees latitude (max along the great circle arc)
        out2 = uxgrid.get_faces_containing_point((0, 25.41))

        nt.assert_array_equal(out1, out2)

…/uxarray into zedwick/point-in-polygon

uxarray/grid/grid.py

uxarray/grid/geometry.py

philipc2 · 2025-02-06T20:22:26Z

uxarray/grid/grid.py

+        if len(point) == 2:
+            point_xyz = np.array(_lonlat_rad_to_xyz(*np.deg2rad(point)))
+            point_lonlat = point
+        elif len(point) == 3:
+            point_xyz = point
+            point_lonlat = np.array(_xyz_to_lonlat_deg(*point_xyz), dtype=np.float64)
+        else:
+            raise ValueError(
+                "Point must either be in spherical or cartesian coordinates."
+            )


When I query with Cartesian coordinates in f64, I get the correct indices.

point = np.array([grid.node_x[0].values, grid.node_y[0].values, grid.node_z[0].values], dtype=np.float64)

Though, with f32 I get a Numba error.

point = np.array([grid.node_x[0].values, grid.node_y[0].values, grid.node_z[0].values], dtype=np.float32)

No implementation of function Function(<built-in function dot>) found for signature: >>> dot(array(float64, 1d, C), array(float32, 1d, C))

philipc2 · 2025-02-06T20:24:35Z

uxarray/grid/grid.py

+        if len(point) == 2:
+            point_xyz = np.array(_lonlat_rad_to_xyz(*np.deg2rad(point)))
+            point_lonlat = point
+        elif len(point) == 3:
+            point_xyz = point
+            point_lonlat = np.array(_xyz_to_lonlat_deg(*point_xyz), dtype=np.float64)
+        else:
+            raise ValueError(
+                "Point must either be in spherical or cartesian coordinates."
+            )


This does appear to be a precision issue when using spherical coordinates, since we internally convert them to cartesian. We should try to avoid as many of these conversions as possible. @hongyuchen1030 and I have encountered this a lot in our previous implementations.

aaronzedwick and others added 20 commits October 24, 2024 13:30

Initial work

75b88e8

Switched to projection method

998aa99

Merge branch 'main' into zedwick/point-in-polygon

7d7be67

Updated test cases, fixed spherical bug

e38ca24

Added test cases

92bab4a

Increased performance

354f821

Update test_geometry.py

5fc59aa

Merge branch 'main' into zedwick/point-in-polygon

c4a4979

updated api

3b09973

Merge branch 'zedwick/point-in-polygon' of https://github.com/UXARRAY…

f7cb51c

…/uxarray into zedwick/point-in-polygon

Merge branch 'main' into zedwick/point-in-polygon

22e84d5

Updated detection method

4043c8b

Merge branch 'main' into zedwick/point-in-polygon

95a0dea

Merge branch 'zedwick/point-in-polygon' of https://github.com/UXARRAY…

5dbed72

…/uxarray into zedwick/point-in-polygon

Implemented new ray casting method

c1eae89

Merge branch 'main' into zedwick/point-in-polygon

586c8ad

Added node crossing and point on edge checks

b33326d

Merge branch 'main' into zedwick/point-in-polygon

0acaa2b

Added docstrings, test case

807ce07

Merge branch 'zedwick/point-in-polygon' of https://github.com/UXARRAY…

b4eba94

…/uxarray into zedwick/point-in-polygon

aaronzedwick changed the title ~~DRAFT: Point In Polygon~~ Point In Polygon Nov 27, 2024

aaronzedwick marked this pull request as ready for review November 27, 2024 15:13

Update test_geometry.py

fa6fd8c

aaronzedwick requested review from philipc2, hongyuchen1030 and rajeeja November 27, 2024 15:15

Fixed point conversion bug

0092258

philipc2 requested changes Nov 27, 2024

View reviewed changes

uxarray/grid/geometry.py Outdated Show resolved Hide resolved

aaronzedwick linked an issue Jan 29, 2025 that may be closed by this pull request

Implement Grid.get_faces_containing_point() #1141

Open

philipc2 reviewed Jan 29, 2025

View reviewed changes

uxarray/grid/grid.py Outdated Show resolved Hide resolved

philipc2 reviewed Jan 29, 2025

View reviewed changes

uxarray/grid/grid.py Outdated Show resolved Hide resolved

Merge branch 'main' into zedwick/point-in-polygon

ebfb0dd

philipc2 removed the run-benchmark Run ASV benchmark workflow label Jan 31, 2025

rajeeja requested changes Feb 3, 2025

View reviewed changes

Merge branch 'main' into zedwick/point-in-polygon

0638cfb

philipc2 and others added 2 commits February 3, 2025 19:30

Merge branch 'main' into zedwick/point-in-polygon

bcb9f17

Merge branch 'main' into zedwick/point-in-polygon

37feaee

rajeeja requested review from rajeeja and philipc2 February 4, 2025 19:16

rajeeja approved these changes Feb 4, 2025

View reviewed changes

philipc2 requested changes Feb 4, 2025

View reviewed changes

philipc2 reviewed Feb 4, 2025

View reviewed changes

benchmarks/mpas_ocean.py Outdated Show resolved Hide resolved

philipc2 reviewed Feb 4, 2025

View reviewed changes

uxarray/grid/geometry.py Outdated Show resolved Hide resolved

philipc2 reviewed Feb 4, 2025

View reviewed changes

philipc2 and others added 6 commits February 4, 2025 14:18

Merge branch 'main' into zedwick/point-in-polygon

fdb28e4

Merge branch 'main' into zedwick/point-in-polygon

568689c

Merge branch 'main' into zedwick/point-in-polygon

44406e2

Added tests, updated numba, fixed bug

3175c53

Merge branch 'zedwick/point-in-polygon' of https://github.com/UXARRAY…

683e2cc

…/uxarray into zedwick/point-in-polygon

Merge branch 'main' into zedwick/point-in-polygon

ef22d37

philipc2 requested changes Feb 6, 2025

View reviewed changes

aaronzedwick and others added 2 commits February 7, 2025 08:23

Merge branch 'main' into zedwick/point-in-polygon

6557c93

Addressed review comments

a93b79a

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Point In Face & Get Faces Containing Point #1056

Point In Face & Get Faces Containing Point #1056

aaronzedwick commented Nov 5, 2024 •

edited

Loading

aaronzedwick commented Nov 27, 2024

philipc2 left a comment

philipc2 left a comment

rajeeja left a comment

aaronzedwick commented Feb 4, 2025 •

edited

Loading

rajeeja left a comment

philipc2 left a comment

philipc2 commented Feb 4, 2025

rajeeja commented Feb 4, 2025

philipc2 commented Feb 4, 2025

philipc2 Feb 4, 2025

aaronzedwick Feb 4, 2025

philipc2 Feb 4, 2025

aaronzedwick Feb 5, 2025

philipc2 Feb 5, 2025

philipc2 Feb 6, 2025

philipc2 Feb 6, 2025

aaronzedwick Feb 6, 2025

aaronzedwick Feb 6, 2025 •

edited

Loading

philipc2 Feb 7, 2025 •

edited

Loading

philipc2 left a comment

philipc2 Feb 6, 2025

philipc2 Feb 6, 2025

Point In Face & Get Faces Containing Point #1056

Are you sure you want to change the base?

Point In Face & Get Faces Containing Point #1056

Conversation

aaronzedwick commented Nov 5, 2024 • edited Loading

Overview

Expected Usage

PR Checklist

aaronzedwick commented Nov 27, 2024

philipc2 left a comment

Choose a reason for hiding this comment

philipc2 left a comment

Choose a reason for hiding this comment

rajeeja left a comment

Choose a reason for hiding this comment

aaronzedwick commented Feb 4, 2025 • edited Loading

rajeeja left a comment

Choose a reason for hiding this comment

philipc2 left a comment

Choose a reason for hiding this comment

philipc2 commented Feb 4, 2025

rajeeja commented Feb 4, 2025

philipc2 commented Feb 4, 2025

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

aaronzedwick Feb 6, 2025 • edited Loading

Choose a reason for hiding this comment

philipc2 Feb 7, 2025 • edited Loading

Choose a reason for hiding this comment

philipc2 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

aaronzedwick commented Nov 5, 2024 •

edited

Loading

aaronzedwick commented Feb 4, 2025 •

edited

Loading

aaronzedwick Feb 6, 2025 •

edited

Loading

philipc2 Feb 7, 2025 •

edited

Loading