GridIntersect does not include cells with identical vertices, even with DISU

[mkennard-aquaveo]

If I am dealing with DISU and a vertex grid, I want GridIntersect to return all cells intersected, but currently it removes cells that share the same vertices as a previously intersected cell. This is fine for layered grids such as with DIS and DISV, but with DISU which has no layers, it is a problem.

GridIntersect ._intersect_point_shapely() includes these lines which cull the results and which I would like to bypass:

if verts in parsed_points:
    continue

GridIntersect._intersect_linestring_shapely() has a similar code block. I would like a way to make that optional so that I can get all intersected cells in the results. Perhaps a kwarg called something like "keep_all_layers" could be used to bypass that code block.

I could submit a pull request if desired.

[dbrakenhoff]

I can see the value of having that option (I even think I started out with a similar option at some point). There are some challenges in getting it implemented though and making it intuitive? Supporting DISU is a slightly odd fit since gridintersect is a 2D operation, which for a DISU grid I guess it means all model cells are projected onto one plane? I guess thats fine, it just requires the user to then parse the result for the desired cells afterwards.

Some questions

• Should initializing the gridintersect object with a DISU grid automatically turn on the option you're referring to?
• How you would differentiate between a point lying on a boundary between two cells (resulting in 2 intersection results for one point) and having multiple identical grid cells in the same plane because of DISU? That's what the current check is there for. The same goes for a linestring exactly on a boundary between two cells.
• Could the check be extended to check not just the vertices but also the grid cellid to address the above?
• If you submit a PR, I'd be happy to take a look! And while you're working on it anyway if you have the time, you could also add the GridIntersect._usg_grid_to_shape_generator?

Let me know what you think! Maybe some others more familiar with DISU have thoughts as well?

[jdhughes-dev]

@dbrakenhoff Could you take a look at this one? I see you have taken a look at it now. I also added @langevin-usgs to take a look.

[mkennard-aquaveo]

- Should initializing the gridintersect object with a DISU grid automatically turn on the option you're referring to?

I initialize the GridIntersect object the same for DISV and DISU, using a VertexGrid, so I think it's better to have the option be at the call to the intersect function.

- How you would differentiate between a point lying on a boundary between two cells (resulting in 2 intersection results for one point) and having multiple identical grid cells in the same plane because of DISU? That's what the current check is there for. The same goes for a linestring exactly on a boundary between two cells.

I think with the new option it should return all the cells, just as if the check wasn't there. If the user wants to cull that list down further by using top/bottom information and determining where the point or line is in 3D, that would be up to them. As you say, GridIntersect is a 2D operation and I'm not suggesting that should change.

- Could the check be extended to check not just the vertices but also the grid cellid to address the above?

We certainly wouldn't want a cellid to be repeated in the results, but I don't think that would happen anyway, would it?

- If you submit a PR, I'd be happy to take a look! And while you're working on it anyway if you have the time, you could also add the GridIntersect._usg_grid_to_shape_generator?

I'm not familiar with GridIntersect._usg_grid_to_shape_generator. Could you tell me more?

…

On Thu, Aug 19, 2021 at 9:52 AM Davíd Brakenhoff ***@***.***> wrote: I can see the value of having that option (I even think I started out with a similar option at some point). There are some challenges in getting it implemented though and making it intuitive? Supporting DISU is a slightly odd fit since gridintersect is a 2D operation, which for a DISU grid I guess it means all model cells are projected onto one plane? I guess thats fine, it just requires the user to then parse the result for the desired cells afterwards. Some questions - Should initializing the gridintersect object with a DISU grid automatically turn on the option you're referring to? - How you would differentiate between a point lying on a boundary between two cells (resulting in 2 intersection results for one point) and having multiple identical grid cells in the same plane because of DISU? That's what the current check is there for. The same goes for a linestring exactly on a boundary between two cells. - Could the check be extended to check not just the vertices but also the grid cellid to address the above? - If you submit a PR, I'd be happy to take a look! And while you're working on it anyway if you have the time, you could also add the GridIntersect._usg_grid_to_shape_generator? Let me know what you think! Maybe some others more familiar with DISU have thoughts as well? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#1198 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AKBA3XUDITXDTHFPZWFQ5XLT5USFNANCNFSM5CMNUFAQ> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email> .

-- Michael Kennard | Aquaveo (www.aquaveo.com)

[dbrakenhoff]

Ah, I thought that since you were working with DISU you would also be dealing with an UnstructuredGrid object (and not a VertexGrid). So in that case, ignore my suggestion about GridIntersect._usg_grid_to_shape_generator. I can pick that up later.

From your answers I gather what you're looking for is a lot simpler than what I was thinking. I think a simple "return all intersections" option is something we can include in the intersect method(s). Perhaps @langevin-usgs has some thoughts on this? But feel free to submit a PR with your proposed option!

[christianlangevin]

It makes sense to me to have an option to return all cells that touch the feature. Since you are working with a VertexGrid right now, that's a simple way to get the intersect capability to do what you want. Since you are using a VertexGrid to represent a 3D DISU model, I guess your cells are overlapping in the 2D plane. It's interesting, and I suppose not surprising, that the VertexGrid intersect works just fine in that case. My understanding @mkennard-aquaveo, is that you will submit a PR that allows an argument to return all cellids that intersect the line or point (I don't think this affects polygons)? I tried to think of a good argument name for that, but came up empty handed.

I guess the long term goal should be to come up with a better 3D intersect approach for an UnstructuredGrid object (though there is no reason we couldn't implement 3D intersect capabilities for structured and vertex grids). Perhaps it is implementing this proposed change if the feature is 2D, or perhaps providing a plane for the intersection? Or doing a full 3D intersection for a 3D feature? Would be good to follow terminology and approaches that is based on an existing library if there is one.

[mkennard-aquaveo]

OK, I submitted a pull request. I went with "keep_stacked_cells" as the name of the new kwarg. Happy to change it if you can think of a better name. I added a couple of tests that verify that it works.

[jdhughes-dev]

@mkennard-aquaveo should the issue be closed?

[mkennard-aquaveo]

I think so.

…

On Tue, Jul 19, 2022 at 12:27 PM jdhughes-usgs ***@***.***> wrote: @mkennard-aquaveo <https://github.com/mkennard-aquaveo> should the issue be closed? — Reply to this email directly, view it on GitHub <#1198 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AKBA3XV5BV3JD37Y7EIMXQDVU3XRTANCNFSM5CMNUFAQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- Michael Kennard | Aquaveo (www.aquaveo.com)