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!, SVPoint_2D
 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.SVPoint_2D((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.SVPoint_2D((x2, y2), T))
     return nothing
 end

src/GeometryOps.jl

@@ -14,9 +14,7 @@ using GeoInterface.Extents: Extents
 
 const GI = GeoInterface
 const GB = GeometryBasics
-
-const TuplePoint{T} = Tuple{T, T} where T <: AbstractFloat
-const Edge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T
+const SA = GeometryBasics.StaticArrays
 
 include("types.jl")
 include("primitives.jl")
@@ -55,6 +53,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 = TuplePoint_2D(GI.getpoint(geom, 1), T)
     # Loop over line segments of line string
     for point₂ in GI.getpoint(geom)
+        x2, y2 = TuplePoint_2D(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 SVPoint_2D((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 SVPoint_2D(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::SVPoint{2, T, false, false}
     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
@@ -90,7 +90,7 @@ function _build_a_list(::Type{T}, poly_a, poly_b; exact) 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 = SVPoint_2D(a_p2, T)
         if i <= 1 || (a_pt1 == a_pt2)  # don't repeat points
             a_pt1 = a_pt2
             continue
@@ -103,7 +103,7 @@ function _build_a_list(::Type{T}, poly_a, poly_b; exact) 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 = SVPoint_2D(b_p2, T)
             if j <= 1 || (b_pt1 == b_pt2)  # don't repeat points
                 b_pt1 = b_pt2
                 continue
@@ -194,7 +194,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 = SVPoint_2D(b_p2, T)
         if i ≤ 1 || (b_pt1 == b_pt2)  # don't repeat points
             b_pt1 = b_pt2
             continue
@@ -543,7 +543,7 @@ function _trace_polynodes(::Type{T}, a_list, b_list, a_idx_list, f_step, poly_a,
     n_a_pts, n_b_pts = length(a_list), length(b_list)
     total_pts = n_a_pts + n_b_pts
     n_cross_pts = length(a_idx_list)
-    return_polys = Vector{_get_poly_type(T)}(undef, 0)
+    return_polys = Vector{PolyType2D{T}}(undef, 0)
     # Keep track of number of processed intersection points
     visited_pts = 0
     processed_pts = 0
@@ -601,11 +601,6 @@ function _trace_polynodes(::Type{T}, a_list, b_list, a_idx_list, f_step, poly_a,
     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; exact)
 
@@ -642,12 +637,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; exact, 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)
@@ -670,7 +665,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
@@ -698,16 +693,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/coverage.jl

@@ -99,12 +99,12 @@ function _coverage(::Type{T}, ring, xmin, xmax, ymin, ymax; exact) where T
         end
     end
     ring_cw = isclockwise(ring)
-    p1 = _tuple_point(GI.getpoint(ring, start_idx), T)
+    p1 = TuplePoint_2D(GI.getpoint(ring, start_idx), T)
     # Must rotate clockwise for the algorithm to work
     point_idx = ring_cw ? Iterators.flatten((start_idx + 1:GI.npoint(ring), 1:start_idx)) :
         Iterators.flatten((start_idx - 1:-1:1, GI.npoint(ring):-1:start_idx))
     for i in point_idx
-        p2 = _tuple_point(GI.getpoint(ring, i), T)
+        p2 = TuplePoint_2D(GI.getpoint(ring, i), T)
         # Determine if edge points are within the cell
         p1_in_cell = _point_in_cell(p1, xmin, xmax, ymin, ymax)
         p2_in_cell = _point_in_cell(p2, xmin, xmax, ymin, ymax)

src/methods/clipping/cut.jl

@@ -69,7 +69,7 @@ function _cut(::Type{T}, ::GI.PolygonTrait, poly, ::GI.LineTrait, line; exact) w
     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; exact)
@@ -103,7 +103,7 @@ function _cut(::Type{T}, geom, line, geom_list, intr_list, n_intr_pts; exact) wh
     _flag_ent_exit!(GI.LineTrait(), line, geom_list; exact)
     # Add first point to output list
     return_coords = [[geom_list[1].point]]
-    cross_backs = [(T(Inf),T(Inf))]
+    cross_backs = [SVPoint_2D((Inf, Inf), T)]
     poly_idx = 1
     n_polys = 1
     # Walk around original polygon to find split polygons

src/methods/clipping/difference.jl

@@ -62,10 +62,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
@@ -76,7 +76,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)
@@ -115,10 +115,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
@@ -134,12 +134,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)}()
+    polys = Vector{PolyType2D{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
@@ -156,7 +156,7 @@ 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