UXARRAY / uxarray

Xarray extension for unstructured climate and global weather data analysis and visualization.
https://uxarray.readthedocs.io/
Apache License 2.0
156 stars 32 forks source link

Point in Polygon #905

Open aaronzedwick opened 2 months ago

aaronzedwick commented 2 months ago

Proposed new feature or change:

Add a function to check if a point is inside a polygon. Currently a pole point in polygon is implemented, possibly modify it and make a similar function? That function makes a reference edge using the point and the equator, which wouldn't work with all points.

philipc2 commented 2 months ago

We need to make sure that we implement the point-in-polygon algorithm for spherical geometry. There are dozens of implementations for standard planar geometry (shapely is the most obvious one), but these wouldn't be ideal.

aaronzedwick commented 2 months ago

We need to make sure that we implement the point-in-polygon algorithm for spherical geometry. There are dozens of implementations for standard planar geometry (shapely is the most obvious one), but these wouldn't be ideal.

That’s why I mentioned @hongyuchen1030 I believe it works on spherical geometry.

hongyuchen1030 commented 2 months ago

We need to make sure that we implement the point-in-polygon algorithm for spherical geometry. There are dozens of implementations for standard planar geometry (shapely is the most obvious one), but these wouldn't be ideal.

That’s why I mentioned @hongyuchen1030 I believe it works on spherical geometry.

Unfortunately, they don't. And it is a complicated problem to solve on the spherical geometry. I remember we have already had such a discussion a long time ago. For the recent one, please see my comments here https://github.com/UXARRAY/uxarray/issues/898#issuecomment-2307768927

How to detect if a point is inside a face on the spherical surface. (I know the immediate answer is ray-casting, but we are on the spherical surface, so the ray will go back!)

for now, I can only come up with the pole_in_polygon that checks pole point only. With some modification, it might be used for a general case. Again this problem is more complicate on the sphere than on the planar surface

aaronzedwick commented 2 months ago

We need to make sure that we implement the point-in-polygon algorithm for spherical geometry. There are dozens of implementations for standard planar geometry (shapely is the most obvious one), but these wouldn't be ideal.

That’s why I mentioned @hongyuchen1030 I believe it works on spherical geometry.

Unfortunately, they don't. And it is a complicated problem to solve on the spherical geometry. I remember we have already had such a discussion a long time ago. For the recent one, please see my comments here #898 (comment)

How to detect if a point is inside a face on the spherical surface. (I know the immediate answer is ray-casting, but we are on the spherical surface, so the ray will go back!)

for now, I can only come up with the pole_in_polygon that checks pole point only. With some modification, it might be used for a general case. Again this problem is more complicate on the sphere than on the planar surface

They don’t? I thought your pole inside polygon worked on a sphere. Would it not work in the general case? Any ideas on how to make it general?

hongyuchen1030 commented 2 months ago

We need to make sure that we implement the point-in-polygon algorithm for spherical geometry. There are dozens of implementations for standard planar geometry (shapely is the most obvious one), but these wouldn't be ideal.

That’s why I mentioned @hongyuchen1030 I believe it works on spherical geometry.

Unfortunately, they don't. And it is a complicated problem to solve on the spherical geometry. I remember we have already had such a discussion a long time ago. For the recent one, please see my comments here #898 (comment)

How to detect if a point is inside a face on the spherical surface. (I know the immediate answer is ray-casting, but we are on the spherical surface, so the ray will go back!)

for now, I can only come up with the pole_in_polygon that checks pole point only. With some modification, it might be used for a general case. Again this problem is more complicate on the sphere than on the planar surface

They don’t? I thought your pole inside polygon worked on a sphere. Would it not work in the general case? Any ideas on how to make it general?

Yes, the pole point inside worked on the sphere. Previously we only needed to query the pole point. I will create a new PR for the general case and put my algorithm there.

aaronzedwick commented 1 week ago

@hongyuchen1030 @philipc2 I think I will start to go ahead and implement something similar as to what is described in this paper here.

hongyuchen1030 commented 1 week ago

@hongyuchen1030 @philipc2 I think I will start to go ahead and implement something similar as to what is described in this paper here.

Me and Paul already talked about this paper and didn't think it's suitable for the Uxarray. Please keep in mind that it uses exact computation (CGAL). It will build a tree diagram as well Image

Before you go into any workaround method (like this one), I suggest learning from the computational geometry, especially the mesh arrangement.

aaronzedwick commented 1 week ago

@hongyuchen1030 @philipc2 I think I will start to go ahead and implement something similar as to what is described in this paper here.

Me and Paul already talked about this paper and didn't think it's suitable for the Uxarray. Please keep in mind that it uses exact computation (CGAL). It will build a tree diagram as well Image

Before you go into any workaround method (like this one), I suggest learning from the computational geometry, especially the mesh arrangement.

Can you elaborate further on what you mean by learning from the computation geometry? I’m happy to pursue any method you think is best. I just need to start working on this (or steps to lead towards completing this) for my remapping algorithm.

hongyuchen1030 commented 1 week ago

Can you elaborate further on what you mean by learning from the computation geometry? I’m happy to pursue any method you think is best. I just need to start working on this (or steps to lead towards completing this) for my remapping algorithm.

Have you tried to look into the mesh arrangement materials that been posted on https://github.com/UXARRAY/uxarray/issues/983?

I just need to start working on this (or steps to lead towards completing this) for my remapping algorithm.

Have you looked into the source codes and the detailed implementation of such algorithms in this paper? There are lots of state-of-art algorithms. However, not every one of them fits the requirement of Uxarray. We are looking for algorithms that works on the spherical surface, accurate enough but without loosing much performance. Again, this packages utilize CGAL exact computation so I don't think it will work for us. If you are interested, you can also learn about the CGAL library.

aaronzedwick commented 1 week ago

@hongyuchen1030 I have, however it seems like a large scale project to tackle, that would stall finishing bilinear remapping for awhile. If you think it is worth it however, I would be happy to start working on that.

Yes I have been looking at the source code for the paper, although I have only had a day to go through it so I don’t have a complete and thorough understanding of how it works yet. However in the paper they seemed to suggest their solution was robust and efficient speed wise.

Ok! Sounds good, I will look into CGAL. Just for clarity, is there a specific reason you think CGAL is incompatible with our project?

hongyuchen1030 commented 1 week ago

@hongyuchen1030 I have, however it seems like a large scale project to tackle, that would stall finishing bilinear remapping for awhile. If you think it is worth it however, I would be happy to start working on that.

I understand that it's a large-scale project, and tackling it now could delay the completion of bilinear remapping. However, I encourage you to read the MeshArrangement paper — you'll realize that we're essentially reinventing the wheel right now but in an unsystematic way. Point-in-polygon on a sphere is a complex problem, and that paper calls it "robust" because of CGAL's exact computation, and "efficient" due to using a tree structure to reduce query complexity to logN. But these claims come from a geoscience perspective, not from computational geometry, so they might not be as solid as they seem.

That said, do you know how long their algorithm would take to run?

In the meantime, if you need a quicker solution for your bilinear mapping, here are two options we’ve discussed before:

  1. Write a helper function that converts each point to the "pole point" and adjusts other node coordinates. Then, call the pole_point_in_polygon function, which is robust for pole points.

  2. Use what Paul has in TempestRemap, but keep in mind that it's not robust. You can test different cases and patch as needed.

Ultimately, I expect we'll have to sort everything out with MeshArrangement. Uxarray is already redoing much of this work. It's a big task, and you could wait until we decide to use it systematically in the future.

Ok! Sounds good, I will look into CGAL. Just for clarity, is there a specific reason you think CGAL is incompatible with our project?

As we've highlighted many times in our meetings, Uxarray does not favor using C++ or any compiled code, and performance is absolutely critical for this project. While CGAL offers robustness, its slow performance doesn’t align with the speed requirements we need, especially for handling large-scale datasets efficiently.

aaronzedwick commented 1 week ago

@hongyuchen1030 I have, however it seems like a large scale project to tackle, that would stall finishing bilinear remapping for awhile. If you think it is worth it however, I would be happy to start working on that.

I understand that it's a large-scale project, and tackling it now could delay the completion of bilinear remapping. However, I encourage you to read the MeshArrangement paper — you'll realize that we're essentially reinventing the wheel right now but in an unsystematic way. Point-in-polygon on a sphere is a complex problem, and that paper calls it "robust" because of CGAL's exact computation, and "efficient" due to using a tree structure to reduce query complexity to logN. But these claims come from a geoscience perspective, not from computational geometry, so they might not be as solid as they seem.

That said, do you know how long their algorithm would take to run?

In the meantime, if you need a quicker solution for your bilinear mapping, here are two options we’ve discussed before:

  1. Write a helper function that converts each point to the "pole point" and adjusts other node coordinates. Then, call the pole_point_in_polygon function, which is robust for pole points.
  2. Use what Paul has in TempestRemap, but keep in mind that it's not robust. You can test different cases and patch as needed.

Ultimately, I expect we'll have to sort everything out with MeshArrangement. Uxarray is already redoing much of this work. It's a big task, and you could wait until we decide to use it systematically in the future.

Ok! Sounds good, I will look into CGAL. Just for clarity, is there a specific reason you think CGAL is incompatible with our project?

As we've highlighted many times in our meetings, Uxarray does not favor using C++ or any compiled code, and performance is absolutely critical for this project. While CGAL offers robustness, its slow performance doesn’t align with the speed requirements we need, especially for handling large-scale datasets efficiently.

Okay! That makes sense! Thanks for the clarification! I was looking into doing something similar to Paul's, but it didn't seem robust like you said, especially as I started trying different implementations. That's why I went looking for something else. Using pole_in_point_polygon sounds like it could be a good solution for now. I will look into that, along with the mesh arrangement. Thanks for the discussion and time!