Debug and add to IntersectingPolygons #107
Conversation
| @@ -59,13 +59,78 @@ struct UnionIntersectingPolygons <: GeometryCorrection end | |||
| application_level(::UnionIntersectingPolygons) = GI.MultiPolygonTrait | |||
|
|
|||
| function (::UnionIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly) | |||
| union_multipoly = if GI.npolygon(multipoly) > 1 | |||
| union_multipoly = tuples(multipoly) | |||
There was a problem hiding this comment.
Here your could do:
| union_multipoly = tuples(multipoly) | |
| union_multipoly = tuples(multipoly; calc_extent=true) |
And use the extent in intersects.
There was a problem hiding this comment.
Amazing @rafaqz ! Yeah, I was wondering if it was that easy to get the extent, so that is great to know. However, what do you think about it being a user input? I think we might need to dispatch off of the Bools as Types stuff @asinghvi17 implemented. With calc_extent = true we have the following:
mp6 = GI.MultiPolygon([[[(-5.0, 10.0), (-5.0, 15.0), (15.0, 15.0), (15.0, 10.0), (-5.0, 10.0)]], [[(-5.0, -5.0), (-5.0, 0.0), (15.0, 0.0), (15.0, -5.0), (-5.0, -5.0)]], [[(10.0, -5.0), (10.0, 15.0), (15.0, 15.0), (15.0, -5.0), (10.0, -5.0)]], [[(-5.0, -5.0), (-5.0, 15.0), (0.0, 15.0), (0.0, -5.0), (-5.0, -5.0)]]])
typeof(GO.tuples(mp6; calc_extent = true))which returns
GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, Nothing}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, Nothing}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, Nothing}
while
typeof(GO.tuples(mp6; calc_extent = false))returns
GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}}, Nothing, Nothing}}, Nothing, Nothing}
Seems like that will cause type instability if we don't dispatch off of user input.
And if that is the case, I would rather have that be its own PR where I add extents to all of the polygon functions I. have written because I don't think I took advantage of it as much as I should have.
There was a problem hiding this comment.
This is a part of the problems listed in #97
There was a problem hiding this comment.
Yeah use BoolAsTypes! its just easier to use true for demonstration purposes
There was a problem hiding this comment.
Are you okay with this being in a different PR? If so, I might just merge and then work on that when I have a bit more time (and once I am done with #106 😆).
There was a problem hiding this comment.
For sure :)
Otherwise its all good so I jumped ahead to optimizing
There was a problem hiding this comment.
ok maybe I can PR with the extent stuff too if you have other priorities, seeing its mostly building on my code in other packages. Unfortunately I've been stuck with a lot of DD/Rasters work the last month and havent had time for looking at clipping performance here yet
There was a problem hiding this comment.
I am also happy to take a stab at it. It just isn't a priority like getting these bugs fixed or getting the points swapped over for me right now. But I am excited for one of us to work on it in the near future haha
| application_level(::DiffIntersectingPolygons) = GI.MultiPolygonTrait | ||
|
|
||
| function (::DiffIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly) | ||
| diff_multipoly = tuples(multipoly) |
There was a problem hiding this comment.
| diff_multipoly = tuples(multipoly) | |
| diff_multipoly = tuples(multipoly; calc_extent=true) |
There was a problem hiding this comment.
In here we should also add an efficiency step - if the extents are known already they should not be recalculated. This might be a bit difficult with how apply is currently implemented, though...
There was a problem hiding this comment.
That should work already? If not its a bug. Like Geointerface.extext will check the extent first and calculate if its not there. Maybe we dont copy it accross tho
There was a problem hiding this comment.
In general for apply I think you can't copy the extent, since there is some transformation going on there.
| poly_disjoint = true # assume current polygon is disjoint from others | ||
| for (next_idx, _) in Iterators.filter(last, Iterators.drop(Iterators.enumerate(keep_idx), curr_idx)) | ||
| next_poly = union_multipoly.geom[next_idx] | ||
| if intersects(curr_poly, next_poly) # if two polygons intersect |
There was a problem hiding this comment.
Then intersects needs an extent check (maybe it already has one?) to rule out non-intersecting extents
Co-authored-by: Anshul Singhvi <[email protected]>
There was a bug in the IntersectingPolygons fix. Unfortunately,
uniondepends on the mutipolygon being valid, so usingunionto make a multipolygon valid just doesn't work and it requires a lot more cross-comparison.I also added a
differencefix and found a bug in my polygon tracing.I think that this could be much improved with the edge sorting and other tricks as mentioned in #65 (even sorting the sub-polygons by extent would probably decrease computation for really large multipolygons) but at least it is a naive start.