Point In Face

[aaronzedwick]

Closes #905

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

[aaronzedwick]

@ 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!

[philipc2]

The point_in_polygon functionality is going to be a great addition.

Can we implement a Grid.get_faces_containing_point() function that returns the indices of the faces that contain a given point?

uxgrid = ... 

# return the indicies of faces that contain the point (0, 90)
face_indices = uxgrid.get_faces_containing_point((0.0, 90.0))

I believe we have discussed this before, and a good approach would be to find the maximum radius of each face and use that to query our BallTree of the face centers. This should narrow down the possible candidates instead of needing to run point_in_polygon on the entire Grid each time.

[hongyuchen1030]

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!

Would it be okay to convert to spherical or cartesian inside the function based on what the user inputs, or should I require both?

I recommend minimizing coordinate conversions as much as possible, as they can introduce unnecessary floating-point errors. Since we already construct and store both Cartesian and spherical coordinates internally in the data array, it’s better to stick to one coordinate system. If the other system is needed, we should first retrieve it from the stored conversions rather than recomputing it.

[aaronzedwick]

The point_in_polygon functionality is going to be a great addition.

Can we implement a Grid.get_faces_containing_point() function that returns the indices of the faces that contain a given point?

uxgrid = ... 

# return the indicies of faces that contain the point (0, 90)
face_indices = uxgrid.get_faces_containing_point((0.0, 90.0))

I believe we have discussed this before, and a good approach would be to find the maximum radius of each face and use that to query our BallTree of the face centers. This should narrow down the possible candidates instead of needing to run point_in_polygon on the entire Grid each time.

This would return a single face always correct? Since faces don't overlap.

[philipc2]

The point_in_polygon functionality is going to be a great addition.
Can we implement a Grid.get_faces_containing_point() function that returns the indices of the faces that contain a given point?

uxgrid = ... 

# return the indicies of faces that contain the point (0, 90)
face_indices = uxgrid.get_faces_containing_point((0.0, 90.0))

I believe we have discussed this before, and a good approach would be to find the maximum radius of each face and use that to query our BallTree of the face centers. This should narrow down the possible candidates instead of needing to run point_in_polygon on the entire Grid each time.

This would return a single face always correct? Since faces don't overlap.

The majority of the time it would only return one, however for the case where the point is exactly on the nodes or edges (i.e. inclusive=True) there would be multiple indices returned, since those can be shared between multiple faces.

[aaronzedwick]

The point_in_polygon functionality is going to be a great addition.
Can we implement a Grid.get_faces_containing_point() function that returns the indices of the faces that contain a given point?

uxgrid = ... 

# return the indicies of faces that contain the point (0, 90)
face_indices = uxgrid.get_faces_containing_point((0.0, 90.0))

I believe we have discussed this before, and a good approach would be to find the maximum radius of each face and use that to query our BallTree of the face centers. This should narrow down the possible candidates instead of needing to run point_in_polygon on the entire Grid each time.

This would return a single face always correct? Since faces don't overlap.

The majority of the time it would only return one, however for the case where the point is exactly on the nodes or edges (i.e. inclusive=True) there would be multiple indices returned, since those can be shared between multiple faces.

Ah, that is true I hadn't thought of that. Thanks for the clarification, I will work on adding that to this PR

[aaronzedwick]

@hongyuchen1030 Could you give the point in the polygon function a review? The first test case test_point_inside is the main test I have been using, which is failing. Philip and I tried debugging it yesterday, but couldn't find any issue with my code. We wanted to see if you saw any obvious misuse of the intersection code.

[hongyuchen1030]

@hongyuchen1030 Could you give the point in the polygon function a review? The first test case test_point_inside is the main test I have been using, which is failing. Philip and I tried debugging it yesterday, but couldn't find any issue with my code. We wanted to see if you saw any obvious misuse of the intersection code.

Looking into it right now

My error messages are

        # Assert that the point is not in the polygon
        self.assertFalse(
>           point_in_polygon(polygon_xyz, polygon_lonlat, point_xyz, point_lonlat, ref_point_xyz, ref_point_lonlat))
E       TypeError: point_in_polygon() takes from 4 to 5 positional arguments but 6 were given

test_geometry.py:1623: TypeError

And only the first one, test_geometry.py::TestPointInPolygon::test_point_inside is passing, all others fail

[aaronzedwick]

@hongyuchen1030 Could you give the point in the polygon function a review? The first test case test_point_inside is the main test I have been using, which is failing. Philip and I tried debugging it yesterday, but couldn't find any issue with my code. We wanted to see if you saw any obvious misuse of the intersection code.

Looking into it right now

My error messages are

        # Assert that the point is not in the polygon
        self.assertFalse(
>           point_in_polygon(polygon_xyz, polygon_lonlat, point_xyz, point_lonlat, ref_point_xyz, ref_point_lonlat))
E       TypeError: point_in_polygon() takes from 4 to 5 positional arguments but 6 were given

test_geometry.py:1623: TypeError

And only the first one, test_geometry.py::TestPointInPolygon::test_point_inside is passing, all others fail

The other test cases are not updated to handle the new conditions of passing in both lon and lat. The passing test case is only passing because I commented out the assertion statement. Sorry for the confusion.

[hongyuchen1030]

@hongyuchen1030 Could you give the point in the polygon function a review? The first test case test_point_inside is the main test I have been using, which is failing. Philip and I tried debugging it yesterday, but couldn't find any issue with my code. We wanted to see if you saw any obvious misuse of the intersection code.

Looking into it right now
My error messages are

        # Assert that the point is not in the polygon
        self.assertFalse(
>           point_in_polygon(polygon_xyz, polygon_lonlat, point_xyz, point_lonlat, ref_point_xyz, ref_point_lonlat))
E       TypeError: point_in_polygon() takes from 4 to 5 positional arguments but 6 were given

test_geometry.py:1623: TypeError

And only the first one, test_geometry.py::TestPointInPolygon::test_point_inside is passing, all others fail

The other test cases are not updated to handle the new conditions of passing in both lon and lat. The passing test case is only passing because I commented out the assertion statement. Sorry for the confusion.

By any chance you know which specific face is not passing?

[aaronzedwick]

@hongyuchen1030 Could you give the point in the polygon function a review? The first test case test_point_inside is the main test I have been using, which is failing. Philip and I tried debugging it yesterday, but couldn't find any issue with my code. We wanted to see if you saw any obvious misuse of the intersection code.

Looking into it right now
My error messages are

        # Assert that the point is not in the polygon
        self.assertFalse(
>           point_in_polygon(polygon_xyz, polygon_lonlat, point_xyz, point_lonlat, ref_point_xyz, ref_point_lonlat))
E       TypeError: point_in_polygon() takes from 4 to 5 positional arguments but 6 were given

test_geometry.py:1623: TypeError

And only the first one, test_geometry.py::TestPointInPolygon::test_point_inside is passing, all others fail

The other test cases are not updated to handle the new conditions of passing in both lon and lat. The passing test case is only passing because I commented out the assertion statement. Sorry for the confusion.

By any chance you know which specific face is not passing?

Half the faces don't pass, and the other half do. The first 3 at least fail, if I remember correctly. Interestingly, with your changes on a gca_gca intersection to make it cartesian only, it works better if the reference point is changed to be next to the pole but not the pole itself. However, there hasn't been any improvement when doing the same thing with the current implementation of the intersection code. So something seems to have changed from the changes you are making in that PR.

[philipc2]

@hongyuchen1030

uxarray/test/test_intersections.py

Lines 95 to 114 in f783f61

	def test_GCA_GCA_pole(self):
	face_lonlat = np.deg2rad(np.array([-175, 26.5]))
	# this fails when the pole is set to exactly -90.0
	ref_point_lonlat = np.deg2rad(np.array([0.0, -89.9]))
	face_xyz = np.array(_lonlat_rad_to_xyz(*face_lonlat))
	ref_point_xyz = np.array(_lonlat_rad_to_xyz(*ref_point_lonlat))
	edge_a_lonlat = np.deg2rad(np.array((-175, -24.5)))
	edge_b_lonlat = np.deg2rad(np.array((-173, 25.7)))
	edge_a_xyz = np.array(_lonlat_rad_to_xyz(*edge_a_lonlat))
	edge_b_xyz = np.array(_lonlat_rad_to_xyz(*edge_b_lonlat))
	gca_a_xyz = np.array([face_xyz, ref_point_xyz])
	gca_b_xyz = np.array([edge_a_xyz, edge_b_xyz])
	# The edge should intersect
	self.assertTrue(len(gca_gca_intersection(gca_a_xyz, gca_b_xyz)))

I've added the test case here in #1112

The correct intersection is computed up to -89.9999999999999, but after that it returns no intersection, such as when the point is set to exactly -90.0

[aaronzedwick]

@hongyuchen1030

uxarray/test/test_intersections.py

Lines 95 to 114 in f783f61

	def test_GCA_GCA_pole(self):
	face_lonlat = np.deg2rad(np.array([-175, 26.5]))
	# this fails when the pole is set to exactly -90.0
	ref_point_lonlat = np.deg2rad(np.array([0.0, -89.9]))
	face_xyz = np.array(_lonlat_rad_to_xyz(*face_lonlat))
	ref_point_xyz = np.array(_lonlat_rad_to_xyz(*ref_point_lonlat))
	edge_a_lonlat = np.deg2rad(np.array((-175, -24.5)))
	edge_b_lonlat = np.deg2rad(np.array((-173, 25.7)))
	edge_a_xyz = np.array(_lonlat_rad_to_xyz(*edge_a_lonlat))
	edge_b_xyz = np.array(_lonlat_rad_to_xyz(*edge_b_lonlat))
	gca_a_xyz = np.array([face_xyz, ref_point_xyz])
	gca_b_xyz = np.array([edge_a_xyz, edge_b_xyz])
	# The edge should intersect
	self.assertTrue(len(gca_gca_intersection(gca_a_xyz, gca_b_xyz)))

I've added the test case here in #1112

The correct intersection is computed up to -89.9999999999999, but after that it returns no intersection, such as when the point is set to exactly -90.0

Thanks for this Philip.

[aaronzedwick]

@hongyuchen1030 using your branch and the code from this branch, it works properly when the reference point is not on the pole, so whatever the issue seems to be fixed in that PR. Thanks for taking the time to look at this for me though!

[hongyuchen1030]

@hongyuchen1030 using your branch and the code from this branch, it works properly when the reference point is not on the pole, so whatever the issue seems to be fixed in that PR. Thanks for taking the time to look at this for me though!

I think me and @philipc2 had resolved this through the recently merged branch #1112 , Thanks again for @philipc2's help. With the updated algorithm, it should work for the the pole point as well now.