svector points

ext/GeometryOpsProjExt/segmentize.jl

@@ -1,14 +1,14 @@
 # This holds the `segmentize` geodesic functionality.
 
-import GeometryOps: GeodesicSegments, _fill_linear_kernel!
+import GeometryOps: GeodesicSegments, _fill_linear_kernel!, _sv_point
 import Proj
 
 function GeometryOps.GeodesicSegments(; max_distance, equatorial_radius::Real=6378137, flattening::Real=1/298.257223563, geodesic::Proj.geod_geodesic = Proj.geod_geodesic(equatorial_radius, flattening))
     return GeometryOps.GeodesicSegments{Proj.geod_geodesic}(geodesic, max_distance)
 end
 
 
-function GeometryOps._fill_linear_kernel!(method::GeodesicSegments{Proj.geod_geodesic}, new_coords::Vector, x1, y1, x2, y2)
+function GeometryOps._fill_linear_kernel!(::Type{T}, method::GeodesicSegments{Proj.geod_geodesic}, new_coords::Vector, x1, y1, x2, y2) where T
     geod_line = Proj.geod_inverseline(method.geodesic, y1, x1, y2, x2)
     # This is the distance in meters computed between the two points.
     # It's `s13` because `geod_inverseline` sets point 3 to the second input point.
@@ -17,11 +17,11 @@ function GeometryOps._fill_linear_kernel!(method::GeodesicSegments{Proj.geod_geo
         n_segments = ceil(Int, distance / method.max_distance)
         for i in 1:(n_segments - 1)
             y, x, _ = Proj.geod_position(geod_line, i / n_segments * distance)
-            push!(new_coords, (x, y))
+            push!(new_coords, GeometryOps._sv_point((x, y), T))
         end
     end
     # End the line with the original coordinate,
     # to avoid any multiplication errors.
-    push!(new_coords, (x2, y2))
+    push!(new_coords, GeometryOps._sv_point((x2, y2), T))
     return nothing
 end

src/GeometryOps.jl

@@ -6,16 +6,21 @@ using GeoInterface
 using GeometryBasics
 import Tables
 using LinearAlgebra, Statistics
-import GeometryBasics.StaticArrays
 import Base.@kwdef
 
 using GeoInterface.Extents: Extents
 
 const GI = GeoInterface
 const GB = GeometryBasics
+const SA = GeometryBasics.StaticArrays
 
 const TuplePoint{T} = Tuple{T, T} where T <: AbstractFloat
-const Edge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T
+const TupleEdge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T
+
+const SVPoint{T} = GI.Point{false, false, SA.SVector{2, T}, Nothing} where T <: AbstractFloat
+const SVEdge{T} = Tuple{SVPoint{T}, SVPoint{T}} where T
+
+const Edge{T} = Union{TupleEdge{T}, SVEdge{T}} where T
 
 include("primitives.jl")
 include("utils.jl")
@@ -50,6 +55,7 @@ include("transformations/flip.jl")
 include("transformations/reproject.jl")
 include("transformations/segmentize.jl")
 include("transformations/simplify.jl")
+include("transformations/svpoints.jl")
 include("transformations/tuples.jl")
 include("transformations/transform.jl")
 include("transformations/correction/geometry_correction.jl")

src/methods/centroid.jl

@@ -79,26 +79,23 @@ function centroid_and_length(
     xcentroid = T(0)
     ycentroid = T(0)
     length = T(0)
-    point₁ = GI.getpoint(geom, 1)
+    x1, y1 = _tuple_point(GI.getpoint(geom, 1), T)
     # Loop over line segments of line string
     for point₂ in GI.getpoint(geom)
+        x2, y2 = _tuple_point(point₂, T)
         # Calculate length of line segment
-        length_component = sqrt(
-            (GI.x(point₂) - GI.x(point₁))^2 +
-            (GI.y(point₂) - GI.y(point₁))^2
-        )
+        length_component = distance((x1, y1), (x2, y2), T)
         # Accumulate the line segment length into `length`
         length += length_component
         # Weighted average of line segment centroids
-        xcentroid += (GI.x(point₁) + GI.x(point₂)) * (length_component / 2)
-        ycentroid += (GI.y(point₁) + GI.y(point₂)) * (length_component / 2)
-        #centroid = centroid .+ ((point₁ .+ point₂) .* (length_component / 2))
+        xcentroid += (x1 + x2) * (length_component / 2)
+        ycentroid += (y1 + y2) * (length_component / 2)
         # Advance the point buffer by 1 point to move to next line segment
-        point₁ = point₂
+        x1, y1 = x2, y2
     end
     xcentroid /= length
     ycentroid /= length
-    return (xcentroid, ycentroid), length
+    return _sv_point((xcentroid, ycentroid), T), length
 end
 
 """
@@ -109,9 +106,10 @@ Returns the centroid and area of a given geometry.
 function centroid_and_area(geom, ::Type{T}=Float64; threaded=false) where T
     target = TraitTarget{Union{GI.PolygonTrait,GI.LineStringTrait,GI.LinearRingTrait}}()
     init = (zero(T), zero(T)), zero(T)
-    applyreduce(_combine_centroid_and_area, target, geom; threaded, init) do g
+    centroid, area = applyreduce(_combine_centroid_and_area, target, geom; threaded, init) do g
         _centroid_and_area(GI.trait(g), g, T)
     end
+    return _sv_point(centroid, T), area
 end
 
 function _centroid_and_area(
@@ -146,14 +144,14 @@ function _centroid_and_area(
 end
 function _centroid_and_area(::GI.PolygonTrait, geom, ::Type{T}) where T
     # Exterior ring's centroid and area
-    (xcentroid, ycentroid), area = centroid_and_area(GI.getexterior(geom), T)
+    (xcentroid, ycentroid), area = _centroid_and_area(GI.LinearRingTrait(), GI.getexterior(geom), T)
     # Weight exterior centroid by area
     xcentroid *= area
     ycentroid *= area
     # Loop over any holes within the polygon
     for hole in GI.gethole(geom)
         # Hole polygon's centroid and area
-        (xinterior, yinterior), interior_area = centroid_and_area(hole, T)
+        (xinterior, yinterior), interior_area = _centroid_and_area(GI.LinearRingTrait(), hole, T)
         # Accumulate the area component into `area`
         area -= interior_area
         # Weighted average of centroid components

src/methods/clipping/clipping_processor.jl

@@ -15,7 +15,7 @@ polygons, or not an endpoint of a chain. =#
 #= This is the struct that makes up a_list and b_list. Many values are only used if point is
 an intersection point (ipt). =#
 @kwdef struct PolyNode{T <: AbstractFloat}
-    point::Tuple{T,T}          # (x, y) values of given point
+    point::GI.Point{false, false, SA.SVector{2, T}, Nothing} # GI.Point(SV[x y]) point values
     inter::Bool = false        # If ipt, true, else 0
     neighbor::Int = 0          # If ipt, index of equivalent point in a_list or b_list, else 0
     idx::Int = 0               # If crossing point, index within sorted a_idx_list
@@ -86,7 +86,7 @@ function _build_a_list(::Type{T}, poly_a, poly_b) where T
     # Loop through points of poly_a
     local a_pt1
     for (i, a_p2) in enumerate(GI.getpoint(poly_a))
-        a_pt2 = (T(GI.x(a_p2)), T(GI.y(a_p2)))
+        a_pt2 = _sv_point(a_p2, T)
         if i <= 1 || (a_pt1 == a_pt2)  # don't repeat points
             a_pt1 = a_pt2
             continue
@@ -99,7 +99,7 @@ function _build_a_list(::Type{T}, poly_a, poly_b) where T
         local b_pt1
         prev_counter = a_count
         for (j, b_p2) in enumerate(GI.getpoint(poly_b))
-            b_pt2 = _tuple_point(b_p2, T)
+            b_pt2 = _sv_point(b_p2, T)
             if j <= 1 || (b_pt1 == b_pt2)  # don't repeat points
                 b_pt1 = b_pt2
                 continue
@@ -190,7 +190,7 @@ function _build_b_list(::Type{T}, a_idx_list, a_list, n_b_intrs, poly_b) where T
     # Loop over points in poly_b and add each point and intersection point
     local b_pt1
     for (i, b_p2) in enumerate(GI.getpoint(poly_b))
-        b_pt2 = _tuple_point(b_p2, T)
+        b_pt2 = _sv_point(b_p2, T)
         if i ≤ 1 || (b_pt1 == b_pt2)  # don't repeat points
             b_pt1 = b_pt2
             continue
@@ -381,7 +381,9 @@ function _flag_ent_exit!(::GI.LinearRingTrait, poly, pt_list, delay_cross_f, del
     # Determine if non-overlapping line midpoint is inside or outside of polygon
     npts = length(pt_list)
     next_idx = start_idx < npts ? (start_idx + 1) : 1
-    start_val = (pt_list[start_idx].point .+ pt_list[next_idx].point) ./ 2
+    start_x = (GI.x(pt_list[start_idx].point) + GI.x(pt_list[next_idx].point)) / 2
+    start_y = (GI.y(pt_list[start_idx].point) + GI.y(pt_list[next_idx].point)) / 2
+    start_val = _sv_point((start_x, start_y))
     start_idx = next_idx - 1  # reset for iterating below
     status = !_point_filled_curve_orientation(start_val, poly; in = true, on = false, out = false)
     # Loop over points and mark entry and exit status
@@ -400,7 +402,9 @@ function _flag_ent_exit!(::GI.LinearRingTrait, poly, pt_list, delay_cross_f, del
                     start_crossing, end_crossing = delay_cross_f(status)
                 else  # delayed bouncing
                     next_idx = ii < npts ? (ii + 1) : 1
-                    next_val = (curr_pt.point .+ pt_list[next_idx].point) ./ 2
+                    next_x = (GI.x(curr_pt.point) + GI.x(pt_list[next_idx].point)) / 2
+                    next_y = (GI.y(curr_pt.point) + GI.y(pt_list[next_idx].point)) / 2
+                    next_val = _sv_point((next_x, next_y))
                     pt_in_poly = _point_filled_curve_orientation(next_val, poly; in = true, on = false, out = false)
                     #= start and end crossing status are the same and depend on if adjacent
                     edges of pt_list are within poly =#
@@ -526,11 +530,6 @@ function _trace_polynodes(::Type{T}, a_list, b_list, a_idx_list, f_step) where T
     return return_polys
 end
 
-# Get type of polygons that will be made
-# TODO: Increase type options
-_get_poly_type(::Type{T}) where T =
-    GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{Tuple{T, T}}, Nothing, Nothing}}, Nothing, Nothing}
-
 #=
     _find_non_cross_orientation(a_list, b_list, a_poly, b_poly)
 
@@ -567,12 +566,12 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol
     # Remove set of holes from all polygons
     for i in 1:n_polys
         n_new_per_poly = 0
-        for curr_hole in Iterators.map(tuples, hole_iterator) # loop through all holes
+        for curr_hole in Iterators.map(Base.Fix2(svpoints, T), hole_iterator) # loop through all holes
             # loop through all pieces of original polygon (new pieces added to end of list)
             for j in Iterators.flatten((i:i, (n_polys + 1):(n_polys + n_new_per_poly)))
                 curr_poly = return_polys[j]
                 remove_poly_idx[j] && continue
-                curr_poly_ext = GI.nhole(curr_poly) > 0 ? GI.Polygon(StaticArrays.SVector(GI.getexterior(curr_poly))) : curr_poly
+                curr_poly_ext = GI.nhole(curr_poly) > 0 ? GI.Polygon(SA.SVector(GI.getexterior(curr_poly))) : curr_poly
                 in_ext, on_ext, out_ext = _line_polygon_interactions(curr_hole, curr_poly_ext; closed_line = true)
                 if in_ext  # hole is at least partially within the polygon's exterior
                     new_hole, new_hole_poly, n_new_pieces = _combine_holes!(T, curr_hole, curr_poly, return_polys, remove_hole_idx)
@@ -595,7 +594,7 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol
                         end
                     end
                 # polygon is completly within hole
-                elseif coveredby(curr_poly_ext, GI.Polygon(StaticArrays.SVector(curr_hole)))
+                elseif coveredby(curr_poly_ext, GI.Polygon(SA.SVector(curr_hole)))
                     remove_poly_idx[j] = true
                 end
             end
@@ -623,16 +622,16 @@ changes.
 function _combine_holes!(::Type{T}, new_hole, curr_poly, return_polys, remove_hole_idx) where T
     n_new_polys = 0
     empty!(remove_hole_idx)
-    new_hole_poly = GI.Polygon(StaticArrays.SVector(new_hole))
+    new_hole_poly = GI.Polygon(SA.SVector(new_hole))
     # Combine any existing holes in curr_poly with new hole
     for (k, old_hole) in enumerate(GI.gethole(curr_poly))
-        old_hole_poly = GI.Polygon(StaticArrays.SVector(old_hole))
+        old_hole_poly = GI.Polygon(SA.SVector(old_hole))
         if intersects(new_hole_poly, old_hole_poly)
             # If the holes intersect, combine them into a bigger hole
             hole_union = union(new_hole_poly, old_hole_poly, T; target = GI.PolygonTrait())[1]
             push!(remove_hole_idx, k + 1)
             new_hole = GI.getexterior(hole_union)
-            new_hole_poly = GI.Polygon(StaticArrays.SVector(new_hole))
+            new_hole_poly = GI.Polygon(SA.SVector(new_hole))
             n_pieces = GI.nhole(hole_union)
             if n_pieces > 0  # if the hole has a hole, then this is a new polygon piece! 
                 append!(return_polys, [GI.Polygon([h]) for h in GI.gethole(hole_union)])

src/methods/clipping/cut.jl

@@ -69,7 +69,7 @@ function _cut(::Type{T}, ::GI.PolygonTrait, poly, ::GI.LineTrait, line) where T
     n_intr_pts = length(intr_list)
     # If an impossible number of intersection points, return original polygon
     if n_intr_pts < 2 || isodd(n_intr_pts)
-        return [tuples(poly)]
+        return [svpoints(poly, T)]
     end
     # Cut polygon by line
     cut_coords = _cut(T, ext_poly, line, poly_list, intr_list, n_intr_pts)
@@ -103,7 +103,7 @@ function _cut(::Type{T}, geom, line, geom_list, intr_list, n_intr_pts) where T
     _flag_ent_exit!(GI.LineTrait(), line, geom_list)
     # Add first point to output list
     return_coords = [[geom_list[1].point]]
-    cross_backs = [(T(Inf),T(Inf))]
+    cross_backs = [_sv_point((Inf, Inf), T)]
     poly_idx = 1
     n_polys = 1
     # Walk around original polygon to find split polygons

src/methods/clipping/difference.jl

@@ -59,10 +59,10 @@ function _difference(
         # add case for if they polygons are the same (all intersection points!)
         # add a find_first check to find first non-inter poly!
         if b_in_a && !a_in_b  # b in a and can't be the same polygon
-            poly_a_b_hole = GI.Polygon([tuples(ext_a), tuples(ext_b)])
+            poly_a_b_hole = GI.Polygon([svpoints(ext_a, T), svpoints(ext_b, T)])
             push!(polys, poly_a_b_hole)
         elseif !b_in_a && !a_in_b # polygons don't intersect
-            push!(polys, tuples(poly_a))
+            push!(polys, svpoints(poly_a, T))
             return polys
         end
     end
@@ -73,7 +73,7 @@ function _difference(
     end
     if GI.nhole(poly_b) != 0
         for hole in GI.gethole(poly_b)
-            hole_poly = GI.Polygon(StaticArrays.SVector(hole))
+            hole_poly = GI.Polygon(SA.SVector(hole))
             new_polys = intersection(hole_poly, poly_a, T; target = GI.PolygonTrait)
             if length(new_polys) > 0
                 append!(polys, new_polys)
@@ -114,10 +114,10 @@ function _difference(
     ::GI.MultiPolygonTrait, multipoly_b;
     kwargs...,
 ) where T
-    polys = [tuples(poly_a, T)]
+    polys = [svpoints(poly_a, T)]
     for poly_b in GI.getpolygon(multipoly_b)
         isempty(polys) && break
-        polys = mapreduce(p -> difference(p, poly_b; target), append!, polys)
+        polys = mapreduce(p -> difference(p, poly_b, T; target), append!, polys)
     end
     return polys
 end
@@ -133,12 +133,12 @@ function _difference(
     fix_multipoly = UnionIntersectingPolygons(), kwargs...,
 ) where T
     if !isnothing(fix_multipoly) # Fix multipoly_a to prevent returning an invalid multipolygon
-        multipoly_a = fix_multipoly(multipoly_a)
+        multipoly_a = fix_multipoly(multipoly_a, T)
     end
     polys = Vector{_get_poly_type(T)}()
     sizehint!(polys, GI.npolygon(multipoly_a))
     for poly_a in GI.getpolygon(multipoly_a)
-        append!(polys, difference(poly_a, poly_b; target))
+        append!(polys, difference(poly_a, poly_b, T; target))
     end
     return polys
 end
@@ -155,12 +155,12 @@ function _difference(
     fix_multipoly = UnionIntersectingPolygons(), kwargs...,
 ) where T
     if !isnothing(fix_multipoly) # Fix multipoly_a to prevent returning an invalid multipolygon
-        multipoly_a = fix_multipoly(multipoly_a)
+        multipoly_a = fix_multipoly(multipoly_a, T)
         fix_multipoly = nothing
     end
     local polys
     for (i, poly_b) in enumerate(GI.getpolygon(multipoly_b))
-        polys = difference(i == 1 ? multipoly_a : GI.MultiPolygon(polys), poly_b; target, fix_multipoly)
+        polys = difference(i == 1 ? multipoly_a : GI.MultiPolygon(polys), poly_b, T; target, fix_multipoly)
     end
     return polys
 end