From 0b62b4969e50f8b6794ec39c1cc1b172bbd81661 Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Sun, 14 Apr 2024 13:00:54 -0700
Subject: [PATCH 01/16] Start switching to StaticArray points

---
 src/GeometryOps.jl                         |  4 ++-
 src/methods/clipping/clipping_processor.jl | 10 ++++----
 src/methods/clipping/difference.jl         |  2 +-
 src/methods/clipping/union.jl              | 12 ++++-----
 src/transformations/svpoints.jl            | 29 ++++++++++++++++++++++
 src/transformations/transform.jl           | 12 ++++++---
 src/transformations/tuples.jl              |  6 ++++-
 src/utils.jl                               |  5 ++++
 8 files changed, 62 insertions(+), 18 deletions(-)
 create mode 100644 src/transformations/svpoints.jl

diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index dc52c7d17..9ace6c03a 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -6,13 +6,14 @@ using GeoInterface
 using GeometryBasics
 import Tables
 using LinearAlgebra, Statistics
-import GeometryBasics.StaticArrays
+# 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
@@ -50,6 +51,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")
diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index ac3d3e06d..e85d15d01 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -571,7 +571,7 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol
             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)
@@ -594,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
@@ -622,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)])
diff --git a/src/methods/clipping/difference.jl b/src/methods/clipping/difference.jl
index f7477e0d2..dcc2bb24f 100644
--- a/src/methods/clipping/difference.jl
+++ b/src/methods/clipping/difference.jl
@@ -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)
diff --git a/src/methods/clipping/union.jl b/src/methods/clipping/union.jl
index e390ef2f1..62be99b6d 100644
--- a/src/methods/clipping/union.jl
+++ b/src/methods/clipping/union.jl
@@ -114,8 +114,8 @@ function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b)
         _add_union_holes_contained_polys!(polys, poly_b, poly_a)
     else  # Polygons intersect, but neither is contained in the other
         n_a_holes = GI.nhole(poly_a)
-        ext_poly_a = GI.Polygon(StaticArrays.SVector(GI.getexterior(poly_a)))
-        ext_poly_b = GI.Polygon(StaticArrays.SVector(GI.getexterior(poly_b)))
+        ext_poly_a = GI.Polygon(SA.SVector(GI.getexterior(poly_a)))
+        ext_poly_b = GI.Polygon(SA.SVector(GI.getexterior(poly_b)))
         #= Start with poly_b when comparing with holes from poly_a and then switch to poly_a
         to compare with holes from poly_b. For current_poly, use ext_poly_b to avoid
         repeating overlapping holes in poly_a and poly_b =#
@@ -135,7 +135,7 @@ function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b)
                 when current_poly is poly_a this includes poly_a holes so overlapping holes
                 between poly_a and poly_b within the overlap are added, in addition to all
                 holes in non-overlapping regions =#
-                h_poly = GI.Polygon(StaticArrays.SVector(ih))
+                h_poly = GI.Polygon(SA.SVector(ih))
                 new_holes = difference(h_poly, current_poly; target = GI.PolygonTrait())
                 append!(polys[1].geom, (GI.getexterior(new_h) for new_h in new_holes))
             end
@@ -155,7 +155,7 @@ function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly)
     union_poly = polys[1]
     ext_int_ring = GI.getexterior(interior_poly)
     for (i, ih) in enumerate(GI.gethole(exterior_poly))
-        poly_ih = GI.Polygon(StaticArrays.SVector(ih))
+        poly_ih = GI.Polygon(SA.SVector(ih))
         in_ih, on_ih, out_ih = _line_polygon_interactions(ext_int_ring, poly_ih; closed_line = true)
         if in_ih  # at least part of interior polygon exterior is within the ith hole
             if !on_ih && !out_ih
@@ -180,13 +180,13 @@ function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly)
             else
                 #= interior polygon's exterior is outside of the ith hole - the interior
                 polygon could either be disjoint from the hole, or contain the hole =#
-                ext_int_poly = GI.Polygon(StaticArrays.SVector(ext_int_ring))
+                ext_int_poly = GI.Polygon(SA.SVector(ext_int_ring))
                 in_int, _, _ = _line_polygon_interactions(ih, ext_int_poly; closed_line = true)
                 if in_int
                     #= interior polygon contains the hole - overlapping holes between the
                     interior and exterior polygons will be added =#
                     for jh in GI.gethole(interior_poly)
-                        poly_jh = GI.Polygon(StaticArrays.SVector(jh))
+                        poly_jh = GI.Polygon(SA.SVector(jh))
                         if intersects(poly_ih, poly_jh)
                             new_holes = intersection(poly_ih, poly_jh; target = GI.PolygonTrait())
                             append!(union_poly.geom, (GI.getexterior(new_h) for new_h in new_holes))
diff --git a/src/transformations/svpoints.jl b/src/transformations/svpoints.jl
new file mode 100644
index 000000000..8fdf7c04d
--- /dev/null
+++ b/src/transformations/svpoints.jl
@@ -0,0 +1,29 @@
+# # SVector conversion
+
+"""
+    sv_points(obj)
+
+Convert all points in `obj` to GI.Points wrapping StaticVectors, wherever the are nested.
+
+Returns a similar object or collection of objects using GeoInterface.jl geometries wrapping
+GI.Points containing a StaticVector.
+
+# Keywords
+
+$APPLY_KEYWORDS
+"""
+function sv_points(geom, ::Type{T} = Float64; kw...) where T
+    if _ismeasured(geom)
+        return apply(PointTrait(), geom; kw...) do p
+            GI.Point(SA.SVector{4}(T(GI.x(p)), T(GI.y(p)), T(GI.z(p)), T(GI.m(p))))
+        end
+    elseif _is3d(geom)
+        return apply(PointTrait(), geom; kw...) do p
+            GI.Point(SA.SVector{3}(T(GI.x(p)), T(GI.y(p)), T(GI.z(p))))
+        end
+    else
+        return apply(PointTrait(), geom; kw...) do p
+            GI.Point(SA.SVector{2}(T(GI.x(p)), T(GI.y(p))))
+        end
+    end
+end
diff --git a/src/transformations/transform.jl b/src/transformations/transform.jl
index b09e85482..17c26772b 100644
--- a/src/transformations/transform.jl
+++ b/src/transformations/transform.jl
@@ -44,14 +44,18 @@ ctor{2, Int64}[[2, 1], [4, 3], [6, 5], [2, 1]], nothing, nothing), GeoInterface.
 }[[4, 3], [6, 5], [7, 6], [4, 3]], nothing, nothing)], nothing, nothing)
 ```
 """
-function transform(f, geom; kw...) 
-    if _is3d(geom)
+function transform(f, geom, ::Type{T} = Float64; kw...) where T
+    if _ismeasured(geom)
         return apply(PointTrait(), geom; kw...) do p
-            f(StaticArrays.SVector{3}((GI.x(p), GI.y(p), GI.z(p))))
+            GI.Point(T.(f(SA.SVector{4}(GI.x(p), GI.y(p), GI.z(p), GI.m(p)))))
+        end
+    elseif _is3d(geom)
+        return apply(PointTrait(), geom; kw...) do p
+            GI.Point(T.(f(SA.SVector{3}((GI.x(p), GI.y(p), GI.z(p))))))
         end
     else
         return apply(PointTrait(), geom; kw...) do p
-            f(StaticArrays.SVector{2}((GI.x(p), GI.y(p))))
+            GI.Point(T.(f(SA.SVector{2}((GI.x(p), GI.y(p))))))
         end
     end
 end
diff --git a/src/transformations/tuples.jl b/src/transformations/tuples.jl
index 4bef4f91d..7ac511857 100644
--- a/src/transformations/tuples.jl
+++ b/src/transformations/tuples.jl
@@ -13,7 +13,11 @@ geometries wrapping `Tuple` points.
 $APPLY_KEYWORDS
 """
 function tuples(geom, ::Type{T} = Float64; kw...) where T
-    if _is3d(geom)
+    if _ismeasured(geom)
+        return apply(PointTrait(), geom; kw...) do p
+            (T(GI.x(p)), T(GI.y(p)), T(GI.z(p)), T(GI.m(p)))
+        end
+    elseif _is3d(geom)
         return apply(PointTrait(), geom; kw...) do p
             (T(GI.x(p)), T(GI.y(p)), T(GI.z(p)))
         end
diff --git a/src/utils.jl b/src/utils.jl
index f59d313b7..74339060e 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,4 +1,9 @@
 # # Utility functions
+_ismeasured(geom)::Bool = _ismeasured(GI.trait(geom), geom)
+_ismeasured(::GI.AbstractGeometryTrait, geom)::Bool = GI.ismeasured(geom)
+_ismeasured(::GI.FeatureTrait, feature)::Bool = _ismeasured(GI.geometry(feature))
+_ismeasured(::GI.FeatureCollectionTrait, fc)::Bool = _ismeasured(GI.getfeature(fc, 1))
+_ismeasured(::Nothing, geom)::Bool = _ismeasured(first(geom)) # Otherwise step into an itererable
 
 _is3d(geom)::Bool = _is3d(GI.trait(geom), geom)
 _is3d(::GI.AbstractGeometryTrait, geom)::Bool = GI.is3d(geom)

From 06a4a8e23cac676a2d7774cfcbd67bb99e2cfdf7 Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Mon, 15 Apr 2024 17:02:43 -0700
Subject: [PATCH 02/16] Update basic transforms to use svector points

---
 Project.toml                               |  1 +
 ext/GeometryOpsProjExt/segmentize.jl       |  8 ++---
 src/methods/clipping/clipping_processor.jl |  5 ----
 src/transformations/flip.jl                | 12 +++++---
 src/transformations/segmentize.jl          | 24 +++++++--------
 src/transformations/simplify.jl            | 34 +++++++++++-----------
 src/transformations/svpoints.jl            |  6 ++--
 src/transformations/transform.jl           | 10 ++-----
 src/utils.jl                               | 11 +++++++
 test/transformations/flip.jl               |  4 +--
 test/transformations/transform.jl          |  2 +-
 11 files changed, 61 insertions(+), 56 deletions(-)

diff --git a/Project.toml b/Project.toml
index bf1184a99..7d0734baa 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
diff --git a/ext/GeometryOpsProjExt/segmentize.jl b/ext/GeometryOpsProjExt/segmentize.jl
index 308c97188..d9d053207 100644
--- a/ext/GeometryOpsProjExt/segmentize.jl
+++ b/ext/GeometryOpsProjExt/segmentize.jl
@@ -1,6 +1,6 @@
 # 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))
@@ -8,7 +8,7 @@ function GeometryOps.GeodesicSegments(; max_distance, equatorial_radius::Real=63
 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
\ No newline at end of file
diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 386782c79..c7c6e226f 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -526,11 +526,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)
 
diff --git a/src/transformations/flip.jl b/src/transformations/flip.jl
index a8b73766f..bc9a9d72a 100644
--- a/src/transformations/flip.jl
+++ b/src/transformations/flip.jl
@@ -14,14 +14,18 @@ original type).
 
 $APPLY_KEYWORDS
 """
-function flip(geom; kw...) 
-    if _is3d(geom)
+function flip(geom; kw...)
+    if _ismeasured(geom)
         return apply(PointTrait(), geom; kw...) do p
-            (GI.y(p), GI.x(p), GI.z(p))
+            GI.Point(SA.SVector{4}(GI.y(p), GI.x(p), GI.z(p), GI.m(p)))
+        end
+    elseif _is3d(geom)
+        return apply(PointTrait(), geom; kw...) do p
+            GI.Point(SA.SVector{3}(GI.y(p), GI.x(p), GI.z(p)))
         end
     else
         return apply(PointTrait(), geom; kw...) do p
-            (GI.y(p), GI.x(p))
+            GI.Point(SA.SVector{2}(GI.y(p), GI.x(p)))
         end
     end
 end
diff --git a/src/transformations/segmentize.jl b/src/transformations/segmentize.jl
index 1e9a47946..09d523fe1 100644
--- a/src/transformations/segmentize.jl
+++ b/src/transformations/segmentize.jl
@@ -168,47 +168,47 @@ This is useful for plotting geometries with a limited number of vertices, or for
 
 Returns a geometry of similar type to the input geometry, but resampled.
 """
-function segmentize(geom; max_distance, threaded::Union{Bool, BoolsAsTypes} = _False())
-    return segmentize(LinearSegments(; max_distance), geom; threaded = _booltype(threaded))
+function segmentize(geom, ::Type{T} = Float64; max_distance, threaded::Union{Bool, BoolsAsTypes} = _False()) where T <: AbstractFloat
+    return segmentize(LinearSegments(; max_distance), geom, T; threaded = _booltype(threaded))
 end
-function segmentize(method::SegmentizeMethod, geom; threaded::Union{Bool, BoolsAsTypes} = _False())
+function segmentize(method::SegmentizeMethod, geom, ::Type{T} = Float64; threaded::Union{Bool, BoolsAsTypes} = _False()) where T <: AbstractFloat
     @assert method.max_distance > 0 "`max_distance` should be positive and nonzero!  Found $(method.max_distance)."
-    segmentize_function = Base.Fix1(_segmentize, method)
+    segmentize_function(g) = _segmentize(T, method, g)
     return apply(segmentize_function, Union{GI.LinearRingTrait, GI.LineStringTrait}, geom; threaded)
 end
 
-_segmentize(method, geom) = _segmentize(method, geom, GI.trait(geom))
+_segmentize(::Type{T}, method, geom) where T = _segmentize(T, method, geom, GI.trait(geom))
 #= 
 This is a method which performs the common functionality for both linear and geodesic algorithms, 
 and calls out to the "kernel" function which we've defined per linesegment.
 =#
-function _segmentize(method::Union{LinearSegments, GeodesicSegments}, geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait})
+function _segmentize(::Type{T}, method::Union{LinearSegments, GeodesicSegments}, geom, ::Union{GI.LineStringTrait, GI.LinearRingTrait}) where T
     first_coord = GI.getpoint(geom, 1)
     x1, y1 = GI.x(first_coord), GI.y(first_coord)
-    new_coords = NTuple{2, Float64}[]
+    new_coords = _get_point_type(T)[] #NTuple{2, Float64}[]
     sizehint!(new_coords, GI.npoint(geom))
-    push!(new_coords, (x1, y1))
+    push!(new_coords, _sv_point((x1, y1), T))
     for coord in Iterators.drop(GI.getpoint(geom), 1)
         x2, y2 = GI.x(coord), GI.y(coord)
-        _fill_linear_kernel!(method, new_coords, x1, y1, x2, y2)
+        _fill_linear_kernel!(T, method, new_coords, x1, y1, x2, y2)
         x1, y1 = x2, y2
     end 
     return rebuild(geom, new_coords)
 end
 
-function _fill_linear_kernel!(method::LinearSegments, new_coords::Vector, x1, y1, x2, y2)
+function _fill_linear_kernel!(::Type{T}, method::LinearSegments, new_coords::Vector, x1, y1, x2, y2) where T
     dx, dy = x2 - x1, y2 - y1
     distance = hypot(dx, dy) # this is a more stable way to compute the Euclidean distance
     if distance > method.max_distance
         n_segments = ceil(Int, distance / method.max_distance)
         for i in 1:(n_segments - 1)
             t = i / n_segments
-            push!(new_coords, (x1 + t * dx, y1 + t * dy))
+            push!(new_coords, _sv_point((x1 + t * dx, y1 + t * dy), T))
         end
     end
     # End the line with the original coordinate,
     # to avoid any multiplication errors.
-    push!(new_coords, (x2, y2))
+    push!(new_coords, _sv_point((x2, y2), T))
     return nothing
 end
 #=
diff --git a/src/transformations/simplify.jl b/src/transformations/simplify.jl
index 01c54e3da..73f163fd6 100644
--- a/src/transformations/simplify.jl
+++ b/src/transformations/simplify.jl
@@ -184,43 +184,43 @@ GI.npoint(simple)
 6
 ```
 """
-simplify(alg::SimplifyAlg, data; kw...) = _simplify(alg, data; kw...)
+simplify(alg::SimplifyAlg, data, ::Type{T} = Float64; kw...) where T = _simplify(T, alg, data; kw...)
 # Default algorithm is DouglasPeucker
 simplify(
-    data; prefilter_alg = nothing,
+    data, ::Type{T} = Float64; prefilter_alg = nothing,
     calc_extent=false, threaded=false, crs=nothing, kw...,
- ) = _simplify(DouglasPeucker(; kw...), data; prefilter_alg, calc_extent, threaded, crs)
+ ) where T = _simplify(T, DouglasPeucker(; kw...), data; prefilter_alg, calc_extent, threaded, crs)
 
 
 #= For each algorithm, apply simplication to all curves, multipoints, and
 points, reconstructing everything else around them. =#
-function _simplify(alg::SimplifyAlg, data; prefilter_alg=nothing, kw...)
-    simplifier(geom) = _simplify(GI.trait(geom), alg, geom; prefilter_alg)
+function _simplify(::Type{T}, alg::SimplifyAlg, data; prefilter_alg=nothing, kw...) where T
+    simplifier(geom) = _simplify(T, GI.trait(geom), alg, geom; prefilter_alg)
     return apply(simplifier, _SIMPLIFY_TARGET, data; kw...)
 end
 
 
 ## For Point and MultiPoint traits we do nothing
-_simplify(::GI.PointTrait, alg, geom; kw...) = geom
-_simplify(::GI.MultiPointTrait, alg, geom; kw...) = geom
+_simplify(::Type{T}, ::GI.PointTrait, alg, geom; kw...) where T = svpoints(geom, T)
+_simplify(::Type{T}, ::GI.MultiPointTrait, alg, geom; kw...) where T = svpoints(geom, T)
 
 ## For curves, rings, and polygon we simplify
 function _simplify(
-    ::GI.AbstractCurveTrait, alg, geom;
+    ::Type{T}, ::GI.AbstractCurveTrait, alg, geom;
     prefilter_alg, preserve_endpoint = true,
-)
+) where T
     points = if isnothing(prefilter_alg)
-        tuple_points(geom)
+        _sv_points(T, geom)
     else
-        _simplify(prefilter_alg, tuple_points(geom), preserve_endpoint)
+        _simplify(prefilter_alg, _sv_points(T, geom), preserve_endpoint)
     end
     return rebuild(geom, _simplify(alg, points, preserve_endpoint))
 end
 
-function _simplify(::GI.PolygonTrait, alg, geom;  kw...)
+function _simplify(::Type{T}, ::GI.PolygonTrait, alg, geom;  kw...) where T
     ## Force treating children as LinearRing
     simplifier(g) = _simplify(
-        GI.LinearRingTrait(), alg, g;
+        T, GI.LinearRingTrait(), alg, g;
         kw..., preserve_endpoint = false,
     )
     rebuilder(g) = rebuild(g, simplifier(g))
@@ -435,7 +435,7 @@ end
 
 # Calculates double the area of a triangle given its vertices
 _triangle_double_area(p1, p2, p3) =
-    abs(p1[1] * (p2[2] - p3[2]) + p2[1] * (p3[2] - p1[2]) + p3[1] * (p1[2] - p2[2]))
+    abs(GI.x(p1) * (GI.y(p2) - GI.y(p3)) + GI.x(p2) * (GI.y(p3) - GI.y(p1)) + GI.x(p3) * (GI.y(p1) - GI.y(p2)))
 
 
 # # Shared utils
@@ -494,10 +494,10 @@ function _build_tolerances(f, points)
     return real_tolerances
 end
 
-function tuple_points(geom)
-    points = Array{Tuple{Float64,Float64}}(undef, GI.npoint(geom))
+function _sv_points(::Type{T}, geom) where T
+    points = Array{_get_point_type(T)}(undef, GI.npoint(geom))
     for (i, p) in enumerate(GI.getpoint(geom))
-        points[i] = (GI.x(p), GI.y(p))
+        points[i] = _sv_point(p)
     end
     return points
 end
diff --git a/src/transformations/svpoints.jl b/src/transformations/svpoints.jl
index 8fdf7c04d..d323c4791 100644
--- a/src/transformations/svpoints.jl
+++ b/src/transformations/svpoints.jl
@@ -1,7 +1,7 @@
-# # SVector conversion
+# # SVector Points conversion
 
 """
-    sv_points(obj)
+    svpoints(obj)
 
 Convert all points in `obj` to GI.Points wrapping StaticVectors, wherever the are nested.
 
@@ -12,7 +12,7 @@ GI.Points containing a StaticVector.
 
 $APPLY_KEYWORDS
 """
-function sv_points(geom, ::Type{T} = Float64; kw...) where T
+function svpoints(geom, ::Type{T} = Float64; kw...) where T
     if _ismeasured(geom)
         return apply(PointTrait(), geom; kw...) do p
             GI.Point(SA.SVector{4}(T(GI.x(p)), T(GI.y(p)), T(GI.z(p)), T(GI.m(p))))
diff --git a/src/transformations/transform.jl b/src/transformations/transform.jl
index 17c26772b..086061d50 100644
--- a/src/transformations/transform.jl
+++ b/src/transformations/transform.jl
@@ -25,10 +25,7 @@ julia> f = CoordinateTransformations.Translation(3.5, 1.5)
 Translation(3.5, 1.5)
 
 julia> GO.transform(f, geom)
-GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Float64}}, Nothing, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.Linea
-rRing{false, false, Vector{StaticArraysCore.SVector{2, Float64}}, Nothing, Nothing}[GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Float64}}, Nothing, Nothing}(StaticArraysCo
-re.SVector{2, Float64}[[4.5, 3.5], [6.5, 5.5], [8.5, 7.5], [4.5, 3.5]], nothing, nothing), GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Float64}}, Nothing, Nothing}(StaticA
-rraysCore.SVector{2, Float64}[[6.5, 5.5], [8.5, 7.5], [9.5, 8.5], [6.5, 5.5]], nothing, nothing)], nothing, nothing)
+GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}}, Nothing, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.LinearRing{false, false, Vector{GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}}, Nothing, Nothing}[GeoInterface.Wrappers.LinearRing{false, false, Vector{GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}[GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([4.5, 3.5], nothing), GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([6.5, 5.5], nothing), GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([8.5, 7.5], nothing), GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([4.5, 3.5], nothing)], nothing, nothing), GeoInterface.Wrappers.LinearRing{false, false, Vector{GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}[GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([6.5, 5.5], nothing), GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([8.5, 7.5], nothing), GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([9.5, 8.5], nothing), GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([6.5, 5.5], nothing)], nothing, nothing)], nothing, nothing)
 ```
 
 With Rotations.jl you need to actuall multiply the Rotation
@@ -38,10 +35,7 @@ by the `SVector` point, which is easy using an anonymous function.
 julia> using Rotations
 
 julia> GO.transform(p -> one(RotMatrix{2}) * p, geom)
-GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Int64}}, Nothing, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.LinearR
-ing{false, false, Vector{StaticArraysCore.SVector{2, Int64}}, Nothing, Nothing}[GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Int64}}, Nothing, Nothing}(StaticArraysCore.SVe
-ctor{2, Int64}[[2, 1], [4, 3], [6, 5], [2, 1]], nothing, nothing), GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Int64}}, Nothing, Nothing}(StaticArraysCore.SVector{2, Int64
-}[[4, 3], [6, 5], [7, 6], [4, 3]], nothing, nothing)], nothing, nothing)
+GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}}, Nothing, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.LinearRing{false, false, Vector{GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}}, Nothing, Nothing}[GeoInterface.Wrappers.LinearRing{false, false, Vector{GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}[GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([1.0, 2.0], nothing), GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([3.0, 4.0], nothing), GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([5.0, 6.0], nothing), GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([1.0, 2.0], nothing)], nothing, nothing), GeoInterface.Wrappers.LinearRing{false, false, Vector{GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}[GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([3.0, 4.0], nothing), GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([5.0, 6.0], nothing), GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([6.0, 7.0], nothing), GeoInterface.Wrappers.Point{false, false, StaticArraysCore.SVector{2, Float64}, Nothing}([3.0, 4.0], nothing)], nothing, nothing)], nothing, nothing)
 ```
 """
 function transform(f, geom, ::Type{T} = Float64; kw...) where T
diff --git a/src/utils.jl b/src/utils.jl
index 74339060e..c4f2b15bf 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -131,3 +131,14 @@ function point_in_extent(p, extent::Extents.Extent)
     (x1, x2), (y1, y1) = extent.X, extent.Y
     return x1 <= GI.x(p) && y1 <= GI.y(p) && x2 >= GI.x(p) && y2 >= GI.y(p)
 end
+
+_get_point_type(::Type{T}) where T =
+    GI.Point{false, false, SA.SVector{2, T}, Nothing}
+
+_sv_point(p) = GI.Point(SA.SVector{2, Float64}(GI.x(p), GI.y(p)))
+_sv_point(p, ::Type{T}) where T = GI.Point(SA.SVector{2, T}(GI.x(p), GI.y(p)))
+# 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}
+
diff --git a/test/transformations/flip.jl b/test/transformations/flip.jl
index eeeb2c848..32c289fdc 100644
--- a/test/transformations/flip.jl
+++ b/test/transformations/flip.jl
@@ -8,6 +8,6 @@ import GeometryOps as GO
                        GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])])
 
 
-    @test GO.flip(geom) == GI.Polygon([GI.LinearRing([(2, 1), (4, 3), (6, 5), (2, 1)]), 
-                                       GI.LinearRing([(4, 3), (6, 5), (7, 6), (4, 3)])])
+    @test GO.equals(GO.flip(geom), GI.Polygon([GI.LinearRing([(2, 1), (4, 3), (6, 5), (2, 1)]), 
+                                       GI.LinearRing([(4, 3), (6, 5), (7, 6), (4, 3)])]))
 end
diff --git a/test/transformations/transform.jl b/test/transformations/transform.jl
index cf4704cc0..4cdc5c5cd 100644
--- a/test/transformations/transform.jl
+++ b/test/transformations/transform.jl
@@ -10,5 +10,5 @@ using CoordinateTransformations
     translated = GI.Polygon([GI.LinearRing([[4.5, 3.5], [6.5, 5.5], [8.5, 7.5], [4.5, 3.5]]), 
                              GI.LinearRing([[6.5, 5.5], [8.5, 7.5], [9.5, 8.5], [6.5, 5.5]])])
     f = CoordinateTransformations.Translation(3.5, 1.5)
-    @test GO.transform(f, geom) == translated
+    @test GO.equals(GO.transform(f, geom), translated)
 end

From afab2c8d029b205b31b4dd60cdb3eaf6a08c39c7 Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Tue, 16 Apr 2024 10:00:06 -0700
Subject: [PATCH 03/16] Update non poly fixes to use svector points

---
 src/transformations/correction/closed_ring.jl | 23 ++++++++-----------
 .../correction/geometry_correction.jl         |  9 ++++----
 .../correction/intersecting_polygons.jl       |  8 +++----
 3 files changed, 18 insertions(+), 22 deletions(-)

diff --git a/src/transformations/correction/closed_ring.jl b/src/transformations/correction/closed_ring.jl
index b8076b821..119586b76 100644
--- a/src/transformations/correction/closed_ring.jl
+++ b/src/transformations/correction/closed_ring.jl
@@ -48,26 +48,21 @@ struct ClosedRing <: GeometryCorrection end
 
 application_level(::ClosedRing) = GI.PolygonTrait
 
-function (::ClosedRing)(::GI.PolygonTrait, polygon)
-    exterior = _close_linear_ring(GI.getexterior(polygon))
-
+function (::ClosedRing)(::Type{T}, ::GI.PolygonTrait, polygon) where T
+    exterior = _close_linear_ring(T, GI.getexterior(polygon))
+    
     holes = map(GI.gethole(polygon)) do hole
-        _close_linear_ring(hole) # TODO: make this more efficient, or use tuples!
+        _close_linear_ring(T, hole) # TODO: make this more efficient, or use tuples!
     end
 
     return GI.Wrappers.Polygon([exterior, holes...])
 end
 
-function _close_linear_ring(ring)
-    if GI.getpoint(ring, 1) == GI.getpoint(ring, GI.npoint(ring))
-        # the ring is closed, all hail the ring
-        return ring
-    else
-        # Assemble the ring as a vector
-        tups = tuples.(GI.getpoint(ring))
+function _close_linear_ring(::Type{T}, ring) where T
+    ring = svpoints(ring, T)
+    if !equals(GI.getpoint(ring, 1), GI.getpoint(ring, GI.npoint(ring)))
         # Close the ring
-        push!(tups, tups[1])
-        # Return an actual ring
-        return GI.LinearRing(tups)
+        push!(ring.geom, ring.geom[1])
     end
+    return ring
 end
\ No newline at end of file
diff --git a/src/transformations/correction/geometry_correction.jl b/src/transformations/correction/geometry_correction.jl
index 5cccb5ea2..8bab902da 100644
--- a/src/transformations/correction/geometry_correction.jl
+++ b/src/transformations/correction/geometry_correction.jl
@@ -40,18 +40,19 @@ abstract type GeometryCorrection end
 
 application_level(gc::GeometryCorrection) = error("Not implemented yet for $(gc)")
 
-(gc::GeometryCorrection)(geometry) = gc(GI.trait(geometry), geometry)
+(gc::GeometryCorrection)(geometry, ::Type{T} = Float64) where T = gc(T, GI.trait(geometry), geometry)
 
-(gc::GeometryCorrection)(trait::GI.AbstractGeometryTrait, geometry) = error("Not implemented yet for $(gc) and $(trait).")
+(gc::GeometryCorrection)(_, trait::GI.AbstractGeometryTrait, geometry) = error("Not implemented yet for $(gc) and $(trait).")
 
-function fix(geometry; corrections = GeometryCorrection[ClosedRing(),], kwargs...)
+function fix(geometry, ::Type{T} = Float64; corrections = GeometryCorrection[ClosedRing(),], kwargs...) where T
     traits = application_level.(corrections)
     final_geometry = geometry
     for Trait in (GI.PointTrait, GI.MultiPointTrait, GI.LineStringTrait, GI.LinearRingTrait, GI.MultiLineStringTrait, GI.PolygonTrait, GI.MultiPolygonTrait)
         available_corrections = findall(x -> x == Trait, traits)
         isempty(available_corrections) && continue
         @debug "Correcting for $(Trait)"
-        net_function = reduce(∘, corrections[available_corrections])
+        @show T
+        net_function = reduce(∘, Base.Fix2.(corrections[available_corrections], T))
         final_geometry = apply(net_function, Trait, final_geometry; kwargs...)
     end
     return final_geometry
diff --git a/src/transformations/correction/intersecting_polygons.jl b/src/transformations/correction/intersecting_polygons.jl
index e8858d45f..17403cc49 100644
--- a/src/transformations/correction/intersecting_polygons.jl
+++ b/src/transformations/correction/intersecting_polygons.jl
@@ -58,8 +58,8 @@ struct UnionIntersectingPolygons <: GeometryCorrection end
 
 application_level(::UnionIntersectingPolygons) = GI.MultiPolygonTrait
 
-function (::UnionIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly)
-    union_multipoly = tuples(multipoly)
+function (::UnionIntersectingPolygons)(::Type{T}, ::GI.MultiPolygonTrait, multipoly) where T
+    union_multipoly = tuples(multipoly, T)
     n_polys = GI.npolygon(multipoly)
     if n_polys > 1
         keep_idx = trues(n_polys)  # keep track of sub-polygons to remove
@@ -101,8 +101,8 @@ struct DiffIntersectingPolygons <: GeometryCorrection end
 
 application_level(::DiffIntersectingPolygons) = GI.MultiPolygonTrait
 
-function (::DiffIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly)
-    diff_multipoly = tuples(multipoly)
+function (::DiffIntersectingPolygons)(::Type{T}, ::GI.MultiPolygonTrait, multipoly) where T
+    diff_multipoly = tuples(multipoly, T)
     n_starting_polys = GI.npolygon(multipoly)
     n_polys = n_starting_polys
     if n_polys > 1

From 29cd3466dc4f3613ea017281380ef843cd0938dc Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Tue, 16 Apr 2024 11:45:23 -0700
Subject: [PATCH 04/16] Update polygon clipping to use svector points

---
 src/GeometryOps.jl                            |  7 ++-
 src/methods/clipping/clipping_processor.jl    | 18 ++++---
 src/methods/clipping/cut.jl                   |  4 +-
 src/methods/clipping/difference.jl            | 16 +++---
 src/methods/clipping/intersection.jl          | 40 +++++++-------
 src/methods/clipping/union.jl                 | 52 +++++++++----------
 src/methods/geom_relations/overlaps.jl        |  2 +-
 .../correction/geometry_correction.jl         |  1 -
 .../correction/intersecting_polygons.jl       |  8 +--
 src/utils.jl                                  |  6 +--
 10 files changed, 81 insertions(+), 73 deletions(-)

diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index 9ace6c03a..43f11a392 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -16,7 +16,12 @@ 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")
diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index c7c6e226f..6930bf4ca 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/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 =#
@@ -562,7 +566,7 @@ 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]
diff --git a/src/methods/clipping/cut.jl b/src/methods/clipping/cut.jl
index 2b9a76498..1339dea68 100644
--- a/src/methods/clipping/cut.jl
+++ b/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
diff --git a/src/methods/clipping/difference.jl b/src/methods/clipping/difference.jl
index dcc2bb24f..9f02baa55 100644
--- a/src/methods/clipping/difference.jl
+++ b/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
@@ -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
diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl
index 832e3f787..c64b2e2a4 100644
--- a/src/methods/clipping/intersection.jl
+++ b/src/methods/clipping/intersection.jl
@@ -72,9 +72,9 @@ function _intersection(
     if isempty(polys) # no crossing points, determine if either poly is inside the other
         a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b)
         if a_in_b
-            push!(polys, GI.Polygon([tuples(ext_a)]))
+            push!(polys, GI.Polygon([svpoints(ext_a, T)]))
         elseif b_in_a
-            push!(polys, GI.Polygon([tuples(ext_b)]))
+            push!(polys, GI.Polygon([svpoints(ext_b, T)]))
         end
     end
     remove_idx = falses(length(polys))
@@ -117,11 +117,11 @@ function _intersection(
     fix_multipoly = UnionIntersectingPolygons(), kwargs...,
 ) where T
     if !isnothing(fix_multipoly) # Fix multipoly_b to prevent duplicated intersection regions
-        multipoly_b = fix_multipoly(multipoly_b)
+        multipoly_b = fix_multipoly(multipoly_b, T)
     end
     polys = Vector{_get_poly_type(T)}()
     for poly_b in GI.getpolygon(multipoly_b)
-        append!(polys, intersection(poly_a, poly_b; target))
+        append!(polys, intersection(poly_a, poly_b, T; target))
     end
     return polys
 end
@@ -134,7 +134,7 @@ _intersection(
     ::GI.MultiPolygonTrait, multipoly_a,
     ::GI.PolygonTrait, poly_b;
     kwargs...,
-) where T = intersection(poly_b, multipoly_a; target , kwargs...)
+) where T = intersection(poly_b, multipoly_a, T; target , kwargs...)
 
 #= Multipolygon with multipolygon intersection - note that all intersection regions between
 any sub-polygons of `multipoly_a` and any of the sub-polygons of `multipoly_b` are counted
@@ -148,13 +148,13 @@ function _intersection(
     fix_multipoly = UnionIntersectingPolygons(), kwargs...,
 ) where T
     if !isnothing(fix_multipoly) # Fix both multipolygons to prevent duplicated regions
-        multipoly_a = fix_multipoly(multipoly_a)
-        multipoly_b = fix_multipoly(multipoly_b)
+        multipoly_a = fix_multipoly(multipoly_a, T)
+        multipoly_b = fix_multipoly(multipoly_b, T)
         fix_multipoly = nothing
     end
     polys = Vector{_get_poly_type(T)}()
     for poly_a in GI.getpolygon(multipoly_a)
-        append!(polys, intersection(poly_a, multipoly_b; target, fix_multipoly))
+        append!(polys, intersection(poly_a, multipoly_b, T; target, fix_multipoly))
     end
     return polys
 end
@@ -196,7 +196,7 @@ were possible given geometry extents or if none are found, return an empty list
 GI.Points. =#
 function _intersection_points(::Type{T}, ::GI.AbstractTrait, a, ::GI.AbstractTrait, b) where T
     # Initialize an empty list of points
-    result = GI.Point[]
+    result = Vector{_get_point_type(T)}()
     # Check if the geometries extents even overlap
     Extents.intersects(GI.extent(a), GI.extent(b)) || return result
     # Create a list of edges from the two input geometries
@@ -220,7 +220,7 @@ function _intersection_points(::Type{T}, ::GI.AbstractTrait, a, ::GI.AbstractTra
                 on_b_edge = (!b_closed && j == npoints_b && 0 <= β <= 1) ||
                     (0 <= β < 1)
                 if on_a_edge && on_b_edge
-                    push!(result, GI.Point(point))
+                    push!(result, point)
                 end
             end
         end
@@ -241,9 +241,9 @@ second isn't. Finally, if the intersection is line_over, then both points are va
 are the two points that define the endpoints of the overlapping region between the two
 lines.
 
-Also note again that each intersection is a tuple of two tuples. The first is the
-intersection point (x,y) while the second is the ratio along the initial lines (α, β) for
-that point. 
+Also note again that each intersection is a tuple of two elements. The first is the
+intersection point GI.Point(SV[x,y]) while the second is the ratio along the initial lines
+(α, β) for that point. 
 
 Calculation derivation can be found here: https://stackoverflow.com/questions/563198/ =#
 function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge) where T
@@ -285,10 +285,10 @@ function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge) where T
             T(β_num / r_cross_s)
         end
         # Calculate intersection point using α and β
-        x = T(px + α * rx)
-        y = T(py + α * ry)
+        x = px + α * rx
+        y = py + α * ry
         if 0 ≤ α ≤ 1 && 0 ≤ β ≤ 1  # only a valid intersection is 0 ≤ α, β ≤ 1
-            intr1 = (x, y), (α, β)
+            intr1 = _sv_point((x, y), T), (α, β)
             line_orient = (α == 0 || α == 1 || β == 0 || β == 1) ? line_hinge : line_cross
         end
     elseif sx * Δqp_y == sy * Δqp_x  # if parallel lines are collinear
@@ -304,23 +304,23 @@ function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge) where T
         n_intrs = 0
         if 0 ≤ a1_β ≤ 1
             n_intrs += 1
-            intr1 = (T.(a1), (zero(T), T(a1_β)))
+            intr1 = (_sv_point(a1, T), (zero(T), T(a1_β)))
         end
         if 0 ≤ a2_β ≤ 1
             n_intrs += 1
-            new_intr = (T.(a2), (one(T), T(a2_β)))
+            new_intr = (_sv_point(a2, T), (one(T), T(a2_β)))
             n_intrs == 1 && (intr1 = new_intr)
             n_intrs == 2 && (intr2 = new_intr)
         end
         if 0 < b1_α < 1 
             n_intrs += 1
-            new_intr = (T.(b1), (T(b1_α), zero(T)))
+            new_intr = (_sv_point(b1, T), (T(b1_α), zero(T)))
             n_intrs == 1 && (intr1 = new_intr)
             n_intrs == 2 && (intr2 = new_intr)
         end
         if 0 < b2_α < 1
             n_intrs += 1
-            new_intr = (T.(b2), (T(b2_α), one(T)))
+            new_intr = (_sv_point(b2, T), (T(b2_α), one(T)))
             n_intrs == 1 && (intr1 = new_intr)
             n_intrs == 2 && (intr2 = new_intr)
         end
diff --git a/src/methods/clipping/union.jl b/src/methods/clipping/union.jl
index 62be99b6d..6420c01f4 100644
--- a/src/methods/clipping/union.jl
+++ b/src/methods/clipping/union.jl
@@ -59,12 +59,12 @@ function _union(
     if n_pieces == 0 # no crossing points, determine if either poly is inside the other
         a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b)
         if a_in_b
-            push!(polys, GI.Polygon([tuples(ext_b)]))
+            push!(polys, GI.Polygon([svpoints(ext_b, T)]))
         elseif b_in_a
-            push!(polys,  GI.Polygon([tuples(ext_a)]))
+            push!(polys,  GI.Polygon([svpoints(ext_a, T)]))
         else
-            push!(polys, tuples(poly_a))
-            push!(polys, tuples(poly_b))
+            push!(polys, svpoints(poly_a, T))
+            push!(polys, svpoints(poly_b, T))
             return polys
         end
     elseif n_pieces > 1
@@ -78,7 +78,7 @@ function _union(
     end
     # Add in holes
     if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0
-        _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b)
+        _add_union_holes!(T, polys, a_in_b, b_in_a, poly_a, poly_b)
     end
     # Remove uneeded collinear points on same edge
     for p in polys
@@ -107,11 +107,11 @@ _union_step(x, _) = x ? (-1) : 1
 #= Add holes from two polygons to the exterior polygon formed by their union. If adding the
 the holes reveals that the polygons aren't actually intersecting, return the original
 polygons. =#
-function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b)
+function _add_union_holes!(::Type{T}, polys, a_in_b, b_in_a, poly_a, poly_b) where T
     if a_in_b
-        _add_union_holes_contained_polys!(polys, poly_a, poly_b)
+        _add_union_holes_contained_polys!(T, polys, poly_a, poly_b)
     elseif b_in_a
-        _add_union_holes_contained_polys!(polys, poly_b, poly_a)
+        _add_union_holes_contained_polys!(T, polys, poly_b, poly_a)
     else  # Polygons intersect, but neither is contained in the other
         n_a_holes = GI.nhole(poly_a)
         ext_poly_a = GI.Polygon(SA.SVector(GI.getexterior(poly_a)))
@@ -128,7 +128,7 @@ function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b)
                 #= if the hole isn't in the overlapping region between the two polygons, add
                 the hole to the resulting polygon as we know it can't interact with any
                 other holes =#
-                push!(polys[1].geom, ih)
+                push!(polys[1].geom, svpoints(ih, T))
             else
                 #= if the hole is at least partially in the overlapping region, take the
                 difference of the hole from the polygon it didn't originate from - note that
@@ -136,7 +136,7 @@ function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b)
                 between poly_a and poly_b within the overlap are added, in addition to all
                 holes in non-overlapping regions =#
                 h_poly = GI.Polygon(SA.SVector(ih))
-                new_holes = difference(h_poly, current_poly; target = GI.PolygonTrait())
+                new_holes = difference(h_poly, current_poly, T; target = GI.PolygonTrait())
                 append!(polys[1].geom, (GI.getexterior(new_h) for new_h in new_holes))
             end
             if i == n_a_holes
@@ -151,7 +151,7 @@ end
 #= Add holes holes to the union of two polygons where one of the original polygons was
 inside of the other. If adding the the holes reveal that the polygons aren't actually
 intersecting, return the original polygons.=#
-function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly)
+function _add_union_holes_contained_polys!(::Type{T}, polys, interior_poly, exterior_poly) where T
     union_poly = polys[1]
     ext_int_ring = GI.getexterior(interior_poly)
     for (i, ih) in enumerate(GI.gethole(exterior_poly))
@@ -161,21 +161,21 @@ function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly)
             if !on_ih && !out_ih
                 #= interior polygon is completly within the ith hole - polygons aren't
                 touching and do not actually form a union =#
-                polys[1] = tuples(interior_poly)
-                push!(polys, tuples(exterior_poly))
+                polys[1] = svpoints(interior_poly, T)
+                push!(polys, svpoints(exterior_poly, T))
                 return polys
             else
                 #= interior polygon is partially within the ith hole - area of interior
                 polygon reduces the size of the hole =#
-                new_holes = difference(poly_ih, interior_poly; target = GI.PolygonTrait())
+                new_holes = difference(poly_ih, interior_poly, T; target = GI.PolygonTrait())
                 append!(union_poly.geom, (GI.getexterior(new_h) for new_h in new_holes))
             end
         else  # none of interior polygon exterior is within the ith hole
             if !out_ih
                 #= interior polygon's exterior is the same as the ith hole - polygons do
                 form a union, but do not overlap so all holes stay in final polygon =#
-                append!(union_poly.geom, Iterators.drop(GI.gethole(exterior_poly), i))
-                append!(union_poly.geom, GI.gethole(interior_poly))
+                append!(union_poly.geom, Iterators.map(Base.Fix2(svpoints, T), Iterators.drop(GI.gethole(exterior_poly), i)))
+                append!(union_poly.geom, Iterators.map(Base.Fix2(svpoints, T), GI.gethole(interior_poly)))
                 return polys
             else
                 #= interior polygon's exterior is outside of the ith hole - the interior
@@ -188,14 +188,14 @@ function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly)
                     for jh in GI.gethole(interior_poly)
                         poly_jh = GI.Polygon(SA.SVector(jh))
                         if intersects(poly_ih, poly_jh)
-                            new_holes = intersection(poly_ih, poly_jh; target = GI.PolygonTrait())
+                            new_holes = intersection(poly_ih, poly_jh, T; target = GI.PolygonTrait())
                             append!(union_poly.geom, (GI.getexterior(new_h) for new_h in new_holes))
                         end
                     end
                 else
                     #= interior polygon and the exterior polygon are disjoint - add the ith
                     hole as it is not covered by the interior polygon =#
-                    push!(union_poly.geom, ih)
+                    push!(union_poly.geom, svpoints(ih, T))
                 end
             end
         end
@@ -214,21 +214,21 @@ function _union(
     fix_multipoly = UnionIntersectingPolygons(), kwargs...,
 ) where T
     if !isnothing(fix_multipoly) # Fix multipoly_b to prevent repeated regions in the output
-        multipoly_b = fix_multipoly(multipoly_b)
+        multipoly_b = fix_multipoly(multipoly_b, T)
     end
-    polys = [tuples(poly_a, T)]
+    polys = [svpoints(poly_a, T)]
     for poly_b in GI.getpolygon(multipoly_b)
         if intersects(polys[1], poly_b)
             # If polygons intersect and form a new polygon, swap out polygon
-            new_polys = union(polys[1], poly_b; target)
+            new_polys = union(polys[1], poly_b, T; target)
             if length(new_polys) > 1 # case where they intersect by just one point
-                push!(polys, tuples(poly_b, T))  # add poly_b to list
+                push!(polys, svpoints(poly_b, T))  # add poly_b to list
             else
                 polys[1] = new_polys[1]
             end
         else
             # If they don't intersect, poly_b is now a part of the union as its own polygon
-            push!(polys, tuples(poly_b, T))
+            push!(polys, svpoints(poly_b, T))
         end
     end
     return polys
@@ -241,7 +241,7 @@ _union(
     ::GI.MultiPolygonTrait, multipoly_a,
     ::GI.PolygonTrait, poly_b;
     kwargs...,
-) where T = union(poly_b, multipoly_a; target, kwargs...)
+) where T = union(poly_b, multipoly_a, T; target, kwargs...)
 
 #= Multipolygon with multipolygon union - note that all of the sub-polygons of `multipoly_a`
 and the sub-polygons of `multipoly_b` are included and combined together where there are
@@ -254,13 +254,13 @@ function _union(
     fix_multipoly = UnionIntersectingPolygons(), kwargs...,
 ) where T
     if !isnothing(fix_multipoly) # Fix multipoly_b to prevent repeated regions in the output
-        multipoly_b = fix_multipoly(multipoly_b)
+        multipoly_b = fix_multipoly(multipoly_b, T)
         fix_multipoly = nothing
     end
     multipolys = multipoly_b
     local polys
     for poly_a in GI.getpolygon(multipoly_a)
-        polys = union(poly_a, multipolys; target, fix_multipoly)
+        polys = union(poly_a, multipolys, T; target, fix_multipoly)
         multipolys = GI.MultiPolygon(polys)
     end
     return polys
diff --git a/src/methods/geom_relations/overlaps.jl b/src/methods/geom_relations/overlaps.jl
index 945fe1bb3..5b3329e47 100644
--- a/src/methods/geom_relations/overlaps.jl
+++ b/src/methods/geom_relations/overlaps.jl
@@ -216,7 +216,7 @@ end
 outside of the other edge, return true. Else false. =#
 function _overlaps(
     (a1, a2)::Edge,
-    (b1, b2)::Edge
+    (b1, b2)::Edge,
 )
     # meets in more than one point
     seg_val, _, _ = _intersection_point(Float64, (a1, a2), (b1, b2))
diff --git a/src/transformations/correction/geometry_correction.jl b/src/transformations/correction/geometry_correction.jl
index 8bab902da..14b01e63f 100644
--- a/src/transformations/correction/geometry_correction.jl
+++ b/src/transformations/correction/geometry_correction.jl
@@ -51,7 +51,6 @@ function fix(geometry, ::Type{T} = Float64; corrections = GeometryCorrection[Clo
         available_corrections = findall(x -> x == Trait, traits)
         isempty(available_corrections) && continue
         @debug "Correcting for $(Trait)"
-        @show T
         net_function = reduce(∘, Base.Fix2.(corrections[available_corrections], T))
         final_geometry = apply(net_function, Trait, final_geometry; kwargs...)
     end
diff --git a/src/transformations/correction/intersecting_polygons.jl b/src/transformations/correction/intersecting_polygons.jl
index 17403cc49..2668b88c1 100644
--- a/src/transformations/correction/intersecting_polygons.jl
+++ b/src/transformations/correction/intersecting_polygons.jl
@@ -59,7 +59,7 @@ struct UnionIntersectingPolygons <: GeometryCorrection end
 application_level(::UnionIntersectingPolygons) = GI.MultiPolygonTrait
 
 function (::UnionIntersectingPolygons)(::Type{T}, ::GI.MultiPolygonTrait, multipoly) where T
-    union_multipoly = tuples(multipoly, T)
+    union_multipoly = svpoints(multipoly, T)
     n_polys = GI.npolygon(multipoly)
     if n_polys > 1
         keep_idx = trues(n_polys)  # keep track of sub-polygons to remove
@@ -72,7 +72,7 @@ function (::UnionIntersectingPolygons)(::Type{T}, ::GI.MultiPolygonTrait, multip
                 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
-                        new_polys = union(curr_poly, next_poly; target = GI.PolygonTrait())
+                        new_polys = union(curr_poly, next_poly, T; target = GI.PolygonTrait())
                         n_new_polys = length(new_polys)
                         if n_new_polys == 1  # if polygons combined
                             poly_disjoint = false
@@ -102,7 +102,7 @@ struct DiffIntersectingPolygons <: GeometryCorrection end
 application_level(::DiffIntersectingPolygons) = GI.MultiPolygonTrait
 
 function (::DiffIntersectingPolygons)(::Type{T}, ::GI.MultiPolygonTrait, multipoly) where T
-    diff_multipoly = tuples(multipoly, T)
+    diff_multipoly = svpoints(multipoly, T)
     n_starting_polys = GI.npolygon(multipoly)
     n_polys = n_starting_polys
     if n_polys > 1
@@ -119,7 +119,7 @@ function (::DiffIntersectingPolygons)(::Type{T}, ::GI.MultiPolygonTrait, multipo
                     !keep_idx[curr_piece_idx] && continue
                     curr_poly = diff_multipoly.geom[curr_piece_idx]
                     if intersects(curr_poly, next_poly)  # if two polygons intersect
-                        new_polys = difference(curr_poly, next_poly; target = GI.PolygonTrait())
+                        new_polys = difference(curr_poly, next_poly, T; target = GI.PolygonTrait())
                         n_new_pieces = length(new_polys) - 1
                         if n_new_pieces < 0  # current polygon is covered by next_polygon
                             keep_idx[curr_piece_idx] = false
diff --git a/src/utils.jl b/src/utils.jl
index c4f2b15bf..7e328a14f 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -62,7 +62,7 @@ Convert any geometry or collection of geometries into a flat
 vector of `Tuple{Tuple{Float64,Float64},Tuple{Float64,Float64}}` edges.
 """
 function to_edges(x, ::Type{T} = Float64) where T
-    edges = Vector{Edge{T}}(undef, _nedge(x))
+    edges = Vector{TupleEdge{T}}(undef, _nedge(x))
     _to_edges!(edges, x, 1)
     return edges
 end
@@ -95,7 +95,7 @@ end
 _tuple_point(p) = GI.x(p), GI.y(p)
 _tuple_point(p, ::Type{T}) where T = T(GI.x(p)), T(GI.y(p))
 
-function to_extent(edges::Vector{Edge})
+function to_extent(edges::Vector{<:Edge})
     x, y = extrema(first, edges)
     Extents.Extent(X=x, Y=y)
 end
@@ -140,5 +140,5 @@ _sv_point(p, ::Type{T}) where T = GI.Point(SA.SVector{2, T}(GI.x(p), GI.y(p)))
 # 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}
+    GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{_get_point_type(T)}, Nothing, Nothing}}, Nothing, Nothing}
 

From 9d5c2f868bb25cbc21e05e6423d2edd85a3a2310 Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Tue, 16 Apr 2024 12:43:59 -0700
Subject: [PATCH 05/16] Update centroid to use svpoints

---
 src/methods/centroid.jl            | 24 +++++++++-----------
 test/methods/centroid.jl           | 36 +++++++++++++++---------------
 test/transformations/segmentize.jl |  4 +++-
 3 files changed, 32 insertions(+), 32 deletions(-)

diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl
index 9dd0e4167..8cb01146b 100644
--- a/src/methods/centroid.jl
+++ b/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
diff --git a/test/methods/centroid.jl b/test/methods/centroid.jl
index 5b807eb1d..123501207 100644
--- a/test/methods/centroid.jl
+++ b/test/methods/centroid.jl
@@ -3,16 +3,16 @@
     l1 = LG.LineString([[0.0, 0.0], [10.0, 0.0], [10.0, 10.0]])
     c1, len1 = GO.centroid_and_length(l1)
     c1_from_LG = LG.centroid(l1)
-    @test c1[1] ≈ GI.x(c1_from_LG)
-    @test c1[2] ≈ GI.y(c1_from_LG)
+    @test GI.x(c1) ≈ GI.x(c1_from_LG)
+    @test GI.y(c1) ≈ GI.y(c1_from_LG)
     @test len1 ≈ 20.0
 
     # Spiral line string
     l2 = LG.LineString([[0.0, 0.0], [2.5, -2.5], [-5.0, -3.0], [-4.0, 6.0], [10.0, 10.0], [12.0, -14.56]])
     c2, len2 = GO.centroid_and_length(l2)
     c2_from_LG = LG.centroid(l2)
-    @test c2[1] ≈ GI.x(c2_from_LG) 
-    @test c2[2] ≈ GI.y(c2_from_LG)
+    @test GI.x(c2) ≈ GI.x(c2_from_LG) 
+    @test GI.y(c2) ≈ GI.y(c2_from_LG)
     @test len2 ≈ 59.3090856788928
 
     # Test that non-closed line strings throw an error for centroid_and_area
@@ -22,15 +22,14 @@
     r1 = LG.LinearRing([[0.0, 0.0], [3456.0, 7894.0], [6291.0, 1954.0], [0.0, 0.0]])
     c3 = GO.centroid(r1)
     c3_from_LG = LG.centroid(r1)
-    @test c3[1] ≈ GI.x(c3_from_LG) 
-    @test c3[2] ≈ GI.y(c3_from_LG)
+    @test GI.x(c3) ≈ GI.x(c3_from_LG) 
+    @test GI.y(c3) ≈ GI.y(c3_from_LG)
 end
 
 @testset "Polygons" begin
     # Basic rectangle
     p1 = AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))")
     c1 = GO.centroid(p1)
-    c1 .≈ (5, 5)
     @test GI.x(c1) ≈ 5
     @test GI.y(c1) ≈ 5
 
@@ -41,8 +40,8 @@ end
     ]])
     c2, area2 = GO.centroid_and_area(p2)
     c2_from_LG = LG.centroid(p2)
-    @test c2[1] ≈ GI.x(c2_from_LG)
-    @test c2[2] ≈ GI.y(c2_from_LG)
+    @test GI.x(c2) ≈ GI.x(c2_from_LG)
+    @test GI.y(c2) ≈ GI.y(c2_from_LG)
     @test area2 ≈ LG.area(p2)
 
     # Randomly generated polygon with lots of sides
@@ -55,8 +54,8 @@ end
     ]])
     c3, area3 = GO.centroid_and_area(p3)
     c3_from_LG = LG.centroid(p3)
-    @test c3[1] ≈ GI.x(c3_from_LG)
-    @test c3[2] ≈ GI.y(c3_from_LG)
+    @test GI.x(c3) ≈ GI.x(c3_from_LG)
+    @test GI.y(c3) ≈ GI.y(c3_from_LG)
     @test area3 ≈ LG.area(p3)
 
     # Polygon with one hole
@@ -66,8 +65,8 @@ end
     ])
     c4, area4 = GO.centroid_and_area(p4)
     c4_from_LG = LG.centroid(p4)
-    @test c4[1] ≈ GI.x(c4_from_LG)
-    @test c4[2] ≈ GI.y(c4_from_LG)
+    @test GI.x(c4) ≈ GI.x(c4_from_LG)
+    @test GI.y(c4) ≈ GI.y(c4_from_LG)
     @test area4 ≈ LG.area(p4)
 
     # Polygon with two holes
@@ -78,8 +77,8 @@ end
     ])
     c5 = GO.centroid(p5)
     c5_from_LG = LG.centroid(p5)
-    @test c5[1] ≈ GI.x(c5_from_LG)
-    @test c5[2] ≈ GI.y(c5_from_LG)
+    @test GI.x(c5) ≈ GI.x(c5_from_LG)
+    @test GI.y(c5) ≈ GI.y(c5_from_LG)
 
     # Same polygon as P5 but using a GeoInterface polygon
     p6 = GI.Polygon([
@@ -88,7 +87,8 @@ end
         [(-3.0, -9.0), (3.0, -9.0), (3.0, -8.5), (-3.0, -8.5), (-3.0, -9.0)],
     ])
     c6 = GO.centroid(p6)
-    @test all(c5 .≈ c6)
+    @test GI.x(c5) ≈ GI.x(c6)
+    @test GI.y(c5) ≈ GI.y(c6)
 end
 @testset "MultiPolygons" begin
     # Combine poylgons made above
@@ -104,7 +104,7 @@ end
     ])
     c1, area1 = GO.centroid_and_area(m1)
     c1_from_LG = LG.centroid(m1)
-    @test c1[1] ≈ GI.x(c1_from_LG)
-    @test c1[2] ≈ GI.y(c1_from_LG)
+    @test GI.x(c1) ≈ GI.x(c1_from_LG)
+    @test GI.y(c1) ≈ GI.y(c1_from_LG)
     @test area1 ≈ LG.area(m1)
 end
diff --git a/test/transformations/segmentize.jl b/test/transformations/segmentize.jl
index f34c8b554..b82f27ec3 100644
--- a/test/transformations/segmentize.jl
+++ b/test/transformations/segmentize.jl
@@ -55,7 +55,9 @@ end
     ar = GO.area(lr)
     for max_distance in exp10.(LinRange(log10(0.01), log10(1), 10))
         segmentized = GO.segmentize(GO.LinearSegments(; max_distance), lr)
-        @test all(GO.centroid(segmentized) .≈ ct)
+        seg_ct = GO.centroid(segmentized)
+        @test GI.x(seg_ct) ≈ GI.x(seg_ct)
+        @test GI.y(seg_ct) ≈ GI.y(seg_ct)
         @test GO.area(segmentized) ≈ ar
     end
 end

From 8884057fd98a6a05f4209bfb1f0c0df76eaff2e8 Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Tue, 16 Apr 2024 13:10:31 -0700
Subject: [PATCH 06/16] Update polygonize to use svector point

---
 src/methods/polygonize.jl | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/methods/polygonize.jl b/src/methods/polygonize.jl
index 78a7ac244..6b3477cb8 100644
--- a/src/methods/polygonize.jl
+++ b/src/methods/polygonize.jl
@@ -67,26 +67,26 @@ If `xs` and `ys` are passed in they are used as the pixel center points.
 # Keywords
 - `minpoints`: ignore polygons with less than `minpoints` points. 
 """
-polygonize(A::AbstractMatrix; kw...) = polygonize(axes(A)..., A; kw...) 
+polygonize(A::AbstractMatrix, ::Type{T} = Float64; kw...) where T = polygonize(axes(A)..., A, T; kw...) 
 
-function polygonize(xs, ys, A::AbstractMatrix; minpoints=10)
+function polygonize(xs, ys, A::AbstractMatrix, ::Type{T} = Float64; minpoints=10) where T
     ## This function uses a lazy map to get contours.  
     contours = Iterators.map(get_contours(A)) do contour
-        poly = map(contour) do xy
+        map(contour) do xy
             x, y = Tuple(xy)
-            Point2f(xs[x], ys[y])
+            _sv_point((xs[x], ys[y]), T)
         end
     end
+    @show typeof(contours)
     ## If we filter off the minimum points, then it's a hair more efficient
     ## not to convert contours with length < missingpoints to polygons.
     if minpoints > 1
         contours = Iterators.filter(contours) do contour
             length(contour) > minpoints
         end
-       return map(Polygon, contours)
-    else
-        return map(Polygon, contours)
     end
+    return [GI.Polygon([c]) for c in contours]
+    # map(GI.Polygon, contours)
 end
 
 ## rotate direction clockwise

From 72d2796977f7f26283901e142e478dcabb0d2e16 Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Tue, 16 Apr 2024 13:22:15 -0700
Subject: [PATCH 07/16] Cleanup changes

---
 Project.toml                           | 1 -
 src/GeometryOps.jl                     | 1 -
 src/methods/geom_relations/overlaps.jl | 2 +-
 src/methods/polygonize.jl              | 1 -
 4 files changed, 1 insertion(+), 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index 7d0734baa..bf1184a99 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,7 +8,6 @@ CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index 43f11a392..0a2d70058 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -6,7 +6,6 @@ using GeoInterface
 using GeometryBasics
 import Tables
 using LinearAlgebra, Statistics
-# import GeometryBasics.StaticArrays
 import Base.@kwdef
 
 using GeoInterface.Extents: Extents
diff --git a/src/methods/geom_relations/overlaps.jl b/src/methods/geom_relations/overlaps.jl
index 5b3329e47..945fe1bb3 100644
--- a/src/methods/geom_relations/overlaps.jl
+++ b/src/methods/geom_relations/overlaps.jl
@@ -216,7 +216,7 @@ end
 outside of the other edge, return true. Else false. =#
 function _overlaps(
     (a1, a2)::Edge,
-    (b1, b2)::Edge,
+    (b1, b2)::Edge
 )
     # meets in more than one point
     seg_val, _, _ = _intersection_point(Float64, (a1, a2), (b1, b2))
diff --git a/src/methods/polygonize.jl b/src/methods/polygonize.jl
index 6b3477cb8..6215bbdbe 100644
--- a/src/methods/polygonize.jl
+++ b/src/methods/polygonize.jl
@@ -86,7 +86,6 @@ function polygonize(xs, ys, A::AbstractMatrix, ::Type{T} = Float64; minpoints=10
         end
     end
     return [GI.Polygon([c]) for c in contours]
-    # map(GI.Polygon, contours)
 end
 
 ## rotate direction clockwise

From b679282cc55f4a6acdba0c3515a369c1dab91afc Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Tue, 16 Apr 2024 14:30:26 -0700
Subject: [PATCH 08/16] Remove show call

---
 src/methods/polygonize.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/methods/polygonize.jl b/src/methods/polygonize.jl
index 6215bbdbe..ed0cd5bb9 100644
--- a/src/methods/polygonize.jl
+++ b/src/methods/polygonize.jl
@@ -77,7 +77,6 @@ function polygonize(xs, ys, A::AbstractMatrix, ::Type{T} = Float64; minpoints=10
             _sv_point((xs[x], ys[y]), T)
         end
     end
-    @show typeof(contours)
     ## If we filter off the minimum points, then it's a hair more efficient
     ## not to convert contours with length < missingpoints to polygons.
     if minpoints > 1

From e9847217ebdd587063acd41380e35c10431babab Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Mon, 10 Jun 2024 14:05:21 -0400
Subject: [PATCH 09/16] Switch SVPoints to subtype of static array

---
 src/GeometryOps.jl                         |  8 --------
 src/methods/clipping/clipping_processor.jl |  2 +-
 src/methods/clipping/intersection.jl       |  8 ++++----
 src/primitives.jl                          | 18 ++++++++++++++++++
 src/transformations/flip.jl                |  6 +++---
 src/transformations/simplify.jl            |  2 +-
 src/transformations/svpoints.jl            |  6 +++---
 src/transformations/transform.jl           |  6 +++---
 src/utils.jl                               |  6 ++----
 9 files changed, 35 insertions(+), 27 deletions(-)

diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index d9ee8ca23..6ac68be67 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -16,14 +16,6 @@ const GI = GeoInterface
 const GB = GeometryBasics
 const SA = GeometryBasics.StaticArrays
 
-const TuplePoint{T} = Tuple{T, T} where T <: AbstractFloat
-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")
 
diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 225a8f565..a2f067eb9 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/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::GI.Point{false, false, SA.SVector{2, T}, Nothing} # GI.Point(SV[x y]) point values
+    point::SVPoint_2D{T}
     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
diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl
index c4cde14b3..211c34f1a 100644
--- a/src/methods/clipping/intersection.jl
+++ b/src/methods/clipping/intersection.jl
@@ -183,7 +183,7 @@ end
         ::Nothing,
     }
 
-Return a list of intersection points between two geometries of type GI.Point.
+Return a list of intersection points between two geometries of type SVPoint.
 If no intersection point was possible given geometry extents, returns an empty
 list.
 """
@@ -243,14 +243,14 @@ are the two points that define the endpoints of the overlapping region between t
 lines.
 
 Also note again that each intersection is a tuple of two elements. The first is the
-intersection point GI.Point(SV[x,y]) while the second is the ratio along the initial lines
+intersection point SVPoint((x, y)) while the second is the ratio along the initial lines
 (α, β) for that point. 
 
 Calculation derivation can be found here: https://stackoverflow.com/questions/563198/ =#
 function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge; exact) where T
     # Default answer for no intersection
     line_orient = line_out
-    intr1 = (_sv_point((zero(T), zero(T))), (zero(T), zero(T)))
+    intr1 = (_sv_point((zero(T), zero(T)), T), (zero(T), zero(T)))
     intr2 = intr1
     no_intr_result = (line_orient, intr1, intr2)
     # Seperate out line segment points
@@ -430,7 +430,7 @@ function _find_cross_intersection(::Type{T}, a1, a2, b1, b2, a_ext, b_ext) where
     else
         (a1y + α * Δay + b1y + β * Δby) / 2
     end
-    pt = _sv_point((x, y))
+    pt = _sv_point((x, y), T)
     # Check if point is within segment envelopes and adjust to endpoint if not
     if !_point_in_extent(pt, a_ext) || !_point_in_extent(pt, b_ext)
         pt, α, β = _nearest_endpoint(T, a1, a2, b1, b2)
diff --git a/src/primitives.jl b/src/primitives.jl
index 78b582811..49e832461 100644
--- a/src/primitives.jl
+++ b/src/primitives.jl
@@ -72,6 +72,24 @@ struct _False <: BoolsAsTypes end
 @inline _booltype(x::Bool)::BoolsAsTypes = x ? _True() : _False()
 @inline _booltype(x::BoolsAsTypes)::BoolsAsTypes = x
 
+const TuplePoint{T} = Tuple{T, T} where T <: AbstractFloat
+const TupleEdge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T
+
+struct SVPoint{N,T,Z,M} <: GeometryBasics.StaticArraysCore.StaticVector{N,T}
+    vals::NTuple{N,T}  # TODO: Should Z and M be booleans or BoolsAsTypes?
+end
+const SVPoint_2D{T} = SVPoint{2, T, _False, _False} where T <: AbstractFloat
+const SVPoint_3D{T} = SVPoint{3, T, _True, _False} where T <: AbstractFloat
+const SVPoint_4D{T} = SVPoint{4, T, _True, _True} where T <: AbstractFloat
+
+Base.getindex(p::SVPoint, i::Int64) = p.vals[i]
+SVPoint(vals::NTuple{2, T}) where T = SVPoint_2D{T}(vals)
+SVPoint(vals::NTuple{3, T}) where T = SVPoint_3D{T}(vals)
+SVPoint(vals::NTuple{4, T}) where T = SVPoint_4D{T}(vals)
+
+const SVEdge{T} = Tuple{SVPoint{N,T,Z,M}, SVPoint{N,T,Z,M}} where {N,T,Z,M}
+const Edge{T} = Union{TupleEdge{T}, SVEdge{T}} where T
+
 """
     TraitTarget{T}
 
diff --git a/src/transformations/flip.jl b/src/transformations/flip.jl
index bc9a9d72a..98a54f8f0 100644
--- a/src/transformations/flip.jl
+++ b/src/transformations/flip.jl
@@ -17,15 +17,15 @@ $APPLY_KEYWORDS
 function flip(geom; kw...)
     if _ismeasured(geom)
         return apply(PointTrait(), geom; kw...) do p
-            GI.Point(SA.SVector{4}(GI.y(p), GI.x(p), GI.z(p), GI.m(p)))
+            svpoints((GI.y(p), GI.x(p), GI.z(p), GI.m(p)))
         end
     elseif _is3d(geom)
         return apply(PointTrait(), geom; kw...) do p
-            GI.Point(SA.SVector{3}(GI.y(p), GI.x(p), GI.z(p)))
+            svpoints((GI.y(p), GI.x(p), GI.z(p)))
         end
     else
         return apply(PointTrait(), geom; kw...) do p
-            GI.Point(SA.SVector{2}(GI.y(p), GI.x(p)))
+            svpoints((GI.y(p), GI.x(p)))
         end
     end
 end
diff --git a/src/transformations/simplify.jl b/src/transformations/simplify.jl
index aaaf3b289..cd550ab3b 100644
--- a/src/transformations/simplify.jl
+++ b/src/transformations/simplify.jl
@@ -500,7 +500,7 @@ end
 function _sv_points(::Type{T}, geom) where T
     points = Array{_get_point_type(T)}(undef, GI.npoint(geom))
     for (i, p) in enumerate(GI.getpoint(geom))
-        points[i] = _sv_point(p)
+        points[i] = _sv_point(p, T)
     end
     return points
 end
diff --git a/src/transformations/svpoints.jl b/src/transformations/svpoints.jl
index d323c4791..4150683c5 100644
--- a/src/transformations/svpoints.jl
+++ b/src/transformations/svpoints.jl
@@ -15,15 +15,15 @@ $APPLY_KEYWORDS
 function svpoints(geom, ::Type{T} = Float64; kw...) where T
     if _ismeasured(geom)
         return apply(PointTrait(), geom; kw...) do p
-            GI.Point(SA.SVector{4}(T(GI.x(p)), T(GI.y(p)), T(GI.z(p)), T(GI.m(p))))
+            SVPoint{4, T, _True, _True}((GI.x(p), GI.y(p), GI.z(p), GI.m(p)))
         end
     elseif _is3d(geom)
         return apply(PointTrait(), geom; kw...) do p
-            GI.Point(SA.SVector{3}(T(GI.x(p)), T(GI.y(p)), T(GI.z(p))))
+            SVPoint{3, T, _True, _False}((GI.x(p), GI.y(p), GI.z(p)))
         end
     else
         return apply(PointTrait(), geom; kw...) do p
-            GI.Point(SA.SVector{2}(T(GI.x(p)), T(GI.y(p))))
+            SVPoint{2, T, _False, _False}((GI.x(p), GI.y(p)))
         end
     end
 end
diff --git a/src/transformations/transform.jl b/src/transformations/transform.jl
index 086061d50..93c9c870a 100644
--- a/src/transformations/transform.jl
+++ b/src/transformations/transform.jl
@@ -41,15 +41,15 @@ GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearR
 function transform(f, geom, ::Type{T} = Float64; kw...) where T
     if _ismeasured(geom)
         return apply(PointTrait(), geom; kw...) do p
-            GI.Point(T.(f(SA.SVector{4}(GI.x(p), GI.y(p), GI.z(p), GI.m(p)))))
+            f(SVPoint_4D{T}(p))
         end
     elseif _is3d(geom)
         return apply(PointTrait(), geom; kw...) do p
-            GI.Point(T.(f(SA.SVector{3}((GI.x(p), GI.y(p), GI.z(p))))))
+            f(SVPoint_3D{T}(p))
         end
     else
         return apply(PointTrait(), geom; kw...) do p
-            GI.Point(T.(f(SA.SVector{2}((GI.x(p), GI.y(p))))))
+            f(SVPoint_2D{T}(p))
         end
     end
 end
diff --git a/src/utils.jl b/src/utils.jl
index e04c39ef0..7b8f35824 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -132,11 +132,9 @@ function _point_in_extent(p, extent::Extents.Extent)
     return x1 ≤ GI.x(p) ≤ x2 && y1 ≤ GI.y(p) ≤ y2
 end
 
-_get_point_type(::Type{T}) where T =
-    GI.Point{false, false, SA.SVector{2, T}, Nothing}
+_get_point_type(::Type{T}) where T = SVPoint_2D{T}
 
-_sv_point(p) = GI.Point(SA.SVector{2, Float64}(GI.x(p), GI.y(p)))
-_sv_point(p, ::Type{T}) where T = GI.Point(SA.SVector{2, T}(GI.x(p), GI.y(p)))
+_sv_point(p, ::Type{T}) where T = SVPoint_2D{T}(_tuple_point(p))
 # Get type of polygons that will be made
 # TODO: Increase type options
 _get_poly_type(::Type{T}) where T =

From 8e545081e7ca53b09be3f773d0cfefd919a70da1 Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Tue, 11 Jun 2024 10:22:46 -0400
Subject: [PATCH 10/16] Experiment with SVPoint constructors

---
 src/primitives.jl                | 3 +++
 src/transformations/svpoints.jl  | 5 +++--
 src/transformations/transform.jl | 7 ++++---
 3 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/src/primitives.jl b/src/primitives.jl
index 49e832461..8dcc409fb 100644
--- a/src/primitives.jl
+++ b/src/primitives.jl
@@ -78,6 +78,9 @@ const TupleEdge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T
 struct SVPoint{N,T,Z,M} <: GeometryBasics.StaticArraysCore.StaticVector{N,T}
     vals::NTuple{N,T}  # TODO: Should Z and M be booleans or BoolsAsTypes?
 end
+
+SVPoint(p) = SVPoint(tuples(p))
+
 const SVPoint_2D{T} = SVPoint{2, T, _False, _False} where T <: AbstractFloat
 const SVPoint_3D{T} = SVPoint{3, T, _True, _False} where T <: AbstractFloat
 const SVPoint_4D{T} = SVPoint{4, T, _True, _True} where T <: AbstractFloat
diff --git a/src/transformations/svpoints.jl b/src/transformations/svpoints.jl
index 4150683c5..1c1ad4a9e 100644
--- a/src/transformations/svpoints.jl
+++ b/src/transformations/svpoints.jl
@@ -3,10 +3,11 @@
 """
     svpoints(obj)
 
-Convert all points in `obj` to GI.Points wrapping StaticVectors, wherever the are nested.
+Convert all points in `obj` to SVPoints, which are a subtype of a StaticVector, wherever
+they are nested.
 
 Returns a similar object or collection of objects using GeoInterface.jl geometries wrapping
-GI.Points containing a StaticVector.
+SVPoints.
 
 # Keywords
 
diff --git a/src/transformations/transform.jl b/src/transformations/transform.jl
index 93c9c870a..1d69030e3 100644
--- a/src/transformations/transform.jl
+++ b/src/transformations/transform.jl
@@ -41,15 +41,16 @@ GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearR
 function transform(f, geom, ::Type{T} = Float64; kw...) where T
     if _ismeasured(geom)
         return apply(PointTrait(), geom; kw...) do p
-            f(SVPoint_4D{T}(p))
+            SVPoint_4D{T}(f(SVPoint_4D{T}(p)))
         end
     elseif _is3d(geom)
         return apply(PointTrait(), geom; kw...) do p
-            f(SVPoint_3D{T}(p))
+            SVPoint_3D{T}(f(SVPoint_3D{T}(p)))
         end
     else
         return apply(PointTrait(), geom; kw...) do p
-            f(SVPoint_2D{T}(p))
+            SVPoint_2D{T}(f(SVPoint_2D{T}(p)))
+            # GI.Point(T.(f(SA.SVector{2}((GI.x(p), GI.y(p))))))
         end
     end
 end

From 201da836d902fedcbaa26daca252f09a25e2f431 Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Tue, 11 Jun 2024 18:54:26 -0400
Subject: [PATCH 11/16] Remove ambiguity

---
 src/methods/clipping/clipping_processor.jl |  2 +-
 src/primitives.jl                          | 20 ++++++++++++--------
 src/transformations/flip.jl                |  6 +++---
 src/transformations/svpoints.jl            |  6 +++---
 src/transformations/transform.jl           |  7 +++----
 src/utils.jl                               |  4 ++--
 6 files changed, 24 insertions(+), 21 deletions(-)

diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index a2f067eb9..72b6828ab 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/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::SVPoint_2D{T}
+    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
diff --git a/src/primitives.jl b/src/primitives.jl
index ff0e070fe..44c81cc69 100644
--- a/src/primitives.jl
+++ b/src/primitives.jl
@@ -60,16 +60,20 @@ struct SVPoint{N,T,Z,M} <: GeometryBasics.StaticArraysCore.StaticVector{N,T}
     vals::NTuple{N,T}  # TODO: Should Z and M be booleans or BoolsAsTypes?
 end
 
-SVPoint(p) = SVPoint(tuples(p))
+Base.getindex(p::SVPoint, i::Int64) = p.vals[i]
 
-const SVPoint_2D{T} = SVPoint{2, T, _False, _False} where T <: AbstractFloat
-const SVPoint_3D{T} = SVPoint{3, T, _True, _False} where T <: AbstractFloat
-const SVPoint_4D{T} = SVPoint{4, T, _True, _True} where T <: AbstractFloat
+SVPoint_2D(vals::NTuple{2,T}) where T = SVPoint{2, T, _False, _False}(vals)
+SVPoint_2D(vals) = SVPoint_2D((GI.x(vals), GI.y(vals)))
 
-Base.getindex(p::SVPoint, i::Int64) = p.vals[i]
-SVPoint(vals::NTuple{2, T}) where T = SVPoint_2D{T}(vals)
-SVPoint(vals::NTuple{3, T}) where T = SVPoint_3D{T}(vals)
-SVPoint(vals::NTuple{4, T}) where T = SVPoint_4D{T}(vals)
+SVPoint_3D(vals::NTuple{3,T}) where T = SVPoint{3, T, _True, _False}(vals)
+SVPoint_3D(vals) = SVPoint_3D((GI.x(vals), GI.y(vals), GI.z(vals)))
+
+SVPoint_4D(vals::NTuple{4,T}) where T = SVPoint{4, T, _True, _True}(vals)
+SVPoint_4D(vals) = SVPoint_4D((GI.x(vals), GI.y(vals), GI.z(vals), GI.m(vals)))
+
+SVPoint(vals::NTuple{2, <:AbstractFloat}) = SVPoint_2D(vals)
+SVPoint(vals::NTuple{3, <:AbstractFloat}) = SVPoint_3D(vals)
+SVPoint(vals::NTuple{4, <:AbstractFloat}) = SVPoint_4D(vals)
 
 const SVEdge{T} = Tuple{SVPoint{N,T,Z,M}, SVPoint{N,T,Z,M}} where {N,T,Z,M}
 const Edge{T} = Union{TupleEdge{T}, SVEdge{T}} where T
diff --git a/src/transformations/flip.jl b/src/transformations/flip.jl
index 98a54f8f0..76ab724ef 100644
--- a/src/transformations/flip.jl
+++ b/src/transformations/flip.jl
@@ -17,15 +17,15 @@ $APPLY_KEYWORDS
 function flip(geom; kw...)
     if _ismeasured(geom)
         return apply(PointTrait(), geom; kw...) do p
-            svpoints((GI.y(p), GI.x(p), GI.z(p), GI.m(p)))
+            SVPoint_4D(StaticArrays.SVector{4}((GI.y(p), GI.x(p), GI.z(p), GI.m(p))))
         end
     elseif _is3d(geom)
         return apply(PointTrait(), geom; kw...) do p
-            svpoints((GI.y(p), GI.x(p), GI.z(p)))
+            SVPoint_3D(StaticArrays.SVector{3}((GI.y(p), GI.x(p), GI.z(p))))
         end
     else
         return apply(PointTrait(), geom; kw...) do p
-            svpoints((GI.y(p), GI.x(p)))
+            SVPoint_2D(StaticArrays.SVector{2}((GI.y(p), GI.x(p))))
         end
     end
 end
diff --git a/src/transformations/svpoints.jl b/src/transformations/svpoints.jl
index 1c1ad4a9e..98269dfae 100644
--- a/src/transformations/svpoints.jl
+++ b/src/transformations/svpoints.jl
@@ -16,15 +16,15 @@ $APPLY_KEYWORDS
 function svpoints(geom, ::Type{T} = Float64; kw...) where T
     if _ismeasured(geom)
         return apply(PointTrait(), geom; kw...) do p
-            SVPoint{4, T, _True, _True}((GI.x(p), GI.y(p), GI.z(p), GI.m(p)))
+            SVPoint_4D(p)
         end
     elseif _is3d(geom)
         return apply(PointTrait(), geom; kw...) do p
-            SVPoint{3, T, _True, _False}((GI.x(p), GI.y(p), GI.z(p)))
+            SVPoint_3D(p)
         end
     else
         return apply(PointTrait(), geom; kw...) do p
-            SVPoint{2, T, _False, _False}((GI.x(p), GI.y(p)))
+            SVPoint_2D(p)
         end
     end
 end
diff --git a/src/transformations/transform.jl b/src/transformations/transform.jl
index 1d69030e3..47005e352 100644
--- a/src/transformations/transform.jl
+++ b/src/transformations/transform.jl
@@ -41,16 +41,15 @@ GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearR
 function transform(f, geom, ::Type{T} = Float64; kw...) where T
     if _ismeasured(geom)
         return apply(PointTrait(), geom; kw...) do p
-            SVPoint_4D{T}(f(SVPoint_4D{T}(p)))
+            SVPoint_4D(f(StaticArrays.SVector{4}((GI.x(p), GI.y(p), GI.z(p), GI.m(p)))))
         end
     elseif _is3d(geom)
         return apply(PointTrait(), geom; kw...) do p
-            SVPoint_3D{T}(f(SVPoint_3D{T}(p)))
+            SVPoint_3D(f(StaticArrays.SVector{3}((GI.x(p), GI.y(p), GI.z(p)))))
         end
     else
         return apply(PointTrait(), geom; kw...) do p
-            SVPoint_2D{T}(f(SVPoint_2D{T}(p)))
-            # GI.Point(T.(f(SA.SVector{2}((GI.x(p), GI.y(p))))))
+            SVPoint_2D(f(StaticArrays.SVector{2}((GI.x(p), GI.y(p)))))
         end
     end
 end
diff --git a/src/utils.jl b/src/utils.jl
index 7b8f35824..d7826eff5 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -132,9 +132,9 @@ function _point_in_extent(p, extent::Extents.Extent)
     return x1 ≤ GI.x(p) ≤ x2 && y1 ≤ GI.y(p) ≤ y2
 end
 
-_get_point_type(::Type{T}) where T = SVPoint_2D{T}
+_get_point_type(::Type{T}) where T = SVPoint{2, T, _False, _False}
 
-_sv_point(p, ::Type{T}) where T = SVPoint_2D{T}(_tuple_point(p))
+_sv_point(p, ::Type{T}) where T = SVPoint_2D(_tuple_point(p, T))
 # Get type of polygons that will be made
 # TODO: Increase type options
 _get_poly_type(::Type{T}) where T =

From f0e74fb50e73fcca834074171054768c02cf2884 Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Wed, 12 Jun 2024 13:59:44 -0400
Subject: [PATCH 12/16] Update and consolidate point functions

---
 ext/GeometryOpsLibGEOSExt/simplify.jl         |  1 -
 ext/GeometryOpsProjExt/segmentize.jl          |  6 +--
 src/methods/centroid.jl                       |  8 +--
 src/methods/clipping/clipping_processor.jl    |  6 +--
 src/methods/clipping/coverage.jl              |  4 +-
 src/methods/clipping/cut.jl                   |  2 +-
 src/methods/clipping/intersection.jl          | 50 +++++++++----------
 src/methods/clipping/predicates.jl            |  6 +--
 .../geom_relations/geom_geom_processors.jl    | 26 +++++-----
 src/primitives.jl                             | 26 ----------
 src/transformations/segmentize.jl             |  6 +--
 src/transformations/simplify.jl               |  8 ---
 src/types.jl                                  | 49 ++++++++++++++++++
 src/utils.jl                                  | 30 +++++------
 14 files changed, 121 insertions(+), 107 deletions(-)

diff --git a/ext/GeometryOpsLibGEOSExt/simplify.jl b/ext/GeometryOpsLibGEOSExt/simplify.jl
index c370d7976..df40f9497 100644
--- a/ext/GeometryOpsLibGEOSExt/simplify.jl
+++ b/ext/GeometryOpsLibGEOSExt/simplify.jl
@@ -21,7 +21,6 @@ function GO._simplify(::GI.AbstractGeometryTrait, alg::GO.GEOS, geom; kwargs...)
 end
 
 function GO._simplify(trait::GI.AbstractCurveTrait, alg::GO.GEOS, geom; kw...)
-    # TODO: Not sure what to do about T... LibGEOS only works in Float64 and don't think I can convert while still returning LG object
     Base.invoke(
         GO._simplify,
         Tuple{GI.AbstractGeometryTrait, GO.GEOS, typeof(geom)},
diff --git a/ext/GeometryOpsProjExt/segmentize.jl b/ext/GeometryOpsProjExt/segmentize.jl
index d9d053207..803b9cf29 100644
--- a/ext/GeometryOpsProjExt/segmentize.jl
+++ b/ext/GeometryOpsProjExt/segmentize.jl
@@ -1,6 +1,6 @@
 # This holds the `segmentize` geodesic functionality.
 
-import GeometryOps: GeodesicSegments, _fill_linear_kernel!, _sv_point
+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))
@@ -17,11 +17,11 @@ function GeometryOps._fill_linear_kernel!(::Type{T}, method::GeodesicSegments{Pr
         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, GeometryOps._sv_point((x, y), T))
+            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, GeometryOps._sv_point((x2, y2), T))
+    push!(new_coords, GeometryOps.SVPoint_2D((x2, y2), T))
     return nothing
 end
\ No newline at end of file
diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl
index 8cb01146b..48a84876a 100644
--- a/src/methods/centroid.jl
+++ b/src/methods/centroid.jl
@@ -79,10 +79,10 @@ function centroid_and_length(
     xcentroid = T(0)
     ycentroid = T(0)
     length = T(0)
-    x1, y1 = _tuple_point(GI.getpoint(geom, 1), T)
+    x1, y1 = TuplePoint_2D(GI.getpoint(geom, 1), T)
     # Loop over line segments of line string
     for point₂ in GI.getpoint(geom)
-        x2, y2 = _tuple_point(point₂, T)
+        x2, y2 = TuplePoint_2D(point₂, T)
         # Calculate length of line segment
         length_component = distance((x1, y1), (x2, y2), T)
         # Accumulate the line segment length into `length`
@@ -95,7 +95,7 @@ function centroid_and_length(
     end
     xcentroid /= length
     ycentroid /= length
-    return _sv_point((xcentroid, ycentroid), T), length
+    return SVPoint_2D((xcentroid, ycentroid), T), length
 end
 
 """
@@ -109,7 +109,7 @@ function centroid_and_area(geom, ::Type{T}=Float64; threaded=false) where T
     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
+    return SVPoint_2D(centroid, T), area
 end
 
 function _centroid_and_area(
diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 72b6828ab..7a8ff31e4 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -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 = _sv_point(a_p2, T)
+        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 = _sv_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 = _sv_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
diff --git a/src/methods/clipping/coverage.jl b/src/methods/clipping/coverage.jl
index 597fb470f..f894b8ad3 100644
--- a/src/methods/clipping/coverage.jl
+++ b/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)
diff --git a/src/methods/clipping/cut.jl b/src/methods/clipping/cut.jl
index 9b38ac6a4..8a4d158a7 100644
--- a/src/methods/clipping/cut.jl
+++ b/src/methods/clipping/cut.jl
@@ -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 = [_sv_point((Inf, Inf), T)]
+    cross_backs = [SVPoint_2D((Inf, Inf), T)]
     poly_idx = 1
     n_polys = 1
     # Walk around original polygon to find split polygons
diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl
index 211c34f1a..de397ca95 100644
--- a/src/methods/clipping/intersection.jl
+++ b/src/methods/clipping/intersection.jl
@@ -250,12 +250,12 @@ Calculation derivation can be found here: https://stackoverflow.com/questions/56
 function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge; exact) where T
     # Default answer for no intersection
     line_orient = line_out
-    intr1 = (_sv_point((zero(T), zero(T)), T), (zero(T), zero(T)))
+    intr1 = (SVPoint_2D((zero(T), zero(T)), T), (zero(T), zero(T)))
     intr2 = intr1
     no_intr_result = (line_orient, intr1, intr2)
     # Seperate out line segment points
-    (a1x, a1y), (a2x, a2y) = _tuple_point(a1, T), _tuple_point(a2, T)
-    (b1x, b1y), (b2x, b2y) = _tuple_point(b1, T), _tuple_point(b2, T)
+    (a1x, a1y), (a2x, a2y) = TuplePoint_2D(a1, T), TuplePoint_2D(a2, T)
+    (b1x, b1y), (b2x, b2y) = TuplePoint_2D(b1, T), TuplePoint_2D(b2, T)
     # Check if envelopes of lines intersect
     a_ext = Extent(X = minmax(a1x, a2x), Y = minmax(a1y, a2y))
     b_ext = Extent(X = minmax(b1x, b2x), Y = minmax(b1y, b2y))
@@ -304,18 +304,18 @@ function _find_collinear_intersection(::Type{T}, a1, a2, b1, b2, a_ext, b_ext, n
         line_orient = line_over
         β1 = _clamped_frac(distance(a1, b1, T), b_dist)
         β2 = _clamped_frac(distance(a2, b1, T), b_dist)
-        intr1 = (_sv_point(a1, T), (zero(T), β1))
-        intr2 = (_sv_point(a2, T), (one(T), β2))
+        intr1 = (SVPoint_2D(a1, T), (zero(T), β1))
+        intr2 = (SVPoint_2D(a2, T), (one(T), β2))
     elseif b1_in_a && b2_in_a  # 1st vertex of b and 2nd vertex of b form overlap
         line_orient = line_over
         α1 = _clamped_frac(distance(b1, a1, T), a_dist)
         α2 = _clamped_frac(distance(b2, a1, T), a_dist)
-        intr1 = (_sv_point(b1, T), (α1, zero(T)))
-        intr2 = (_sv_point(b2, T), (α2, one(T)))
+        intr1 = (SVPoint_2D(b1, T), (α1, zero(T)))
+        intr2 = (SVPoint_2D(b2, T), (α2, one(T)))
     elseif a1_in_b && b1_in_a  # 1st vertex of a and 1st vertex of b form overlap
         if equals(a1, b1)
             line_orient = line_hinge
-            intr1 = (_sv_point(a1, T), (zero(T), zero(T)))
+            intr1 = (SVPoint_2D(a1, T), (zero(T), zero(T)))
         else
             line_orient = line_over
             intr1, intr2 = _set_ab_collinear_intrs(T, a1, b1, zero(T), zero(T), a1, b1, a_dist, b_dist)
@@ -323,7 +323,7 @@ function _find_collinear_intersection(::Type{T}, a1, a2, b1, b2, a_ext, b_ext, n
     elseif a1_in_b && b2_in_a  # 1st vertex of a and 2nd vertex of b form overlap
         if equals(a1, b2)
             line_orient = line_hinge
-            intr1 = (_sv_point(a1, T), (zero(T), one(T)))
+            intr1 = (SVPoint_2D(a1, T), (zero(T), one(T)))
         else
             line_orient = line_over
             intr1, intr2 = _set_ab_collinear_intrs(T, a1, b2, zero(T), one(T), a1, b1, a_dist, b_dist) 
@@ -331,7 +331,7 @@ function _find_collinear_intersection(::Type{T}, a1, a2, b1, b2, a_ext, b_ext, n
     elseif a2_in_b && b1_in_a  # 2nd vertex of a and 1st vertex of b form overlap
         if equals(a2, b1)
             line_orient = line_hinge
-            intr1 = (_sv_point(a2, T), (one(T), zero(T)))
+            intr1 = (SVPoint_2D(a2, T), (one(T), zero(T)))
         else
             line_orient = line_over
             intr1, intr2 = _set_ab_collinear_intrs(T, a2, b1, one(T), zero(T), a1, b1, a_dist, b_dist)
@@ -339,7 +339,7 @@ function _find_collinear_intersection(::Type{T}, a1, a2, b1, b2, a_ext, b_ext, n
     elseif a2_in_b && b2_in_a  # 2nd vertex of a and 2nd vertex of b form overlap
         if equals(a2, b2)
             line_orient = line_hinge
-            intr1 = (_sv_point(a2, T), (one(T), one(T)))
+            intr1 = (SVPoint_2D(a2, T), (one(T), one(T)))
         else
             line_orient = line_over
             intr1, intr2 = _set_ab_collinear_intrs(T, a2, b2, one(T), one(T), a1, b1, a_dist, b_dist)
@@ -352,8 +352,8 @@ end
 endpoint of segment (a1, a2) and one endpoint of segment (b1, b2). =#
 _set_ab_collinear_intrs(::Type{T}, a_pt, b_pt, a_pt_α, b_pt_β, a1, b1, a_dist, b_dist) where T =
     (
-        (_sv_point(a_pt, T), (a_pt_α, _clamped_frac(distance(a_pt, b1, T), b_dist))),
-        (_sv_point(b_pt, T), (_clamped_frac(distance(b_pt, a1, T), a_dist), b_pt_β))
+        (SVPoint_2D(a_pt, T), (a_pt_α, _clamped_frac(distance(a_pt, b1, T), b_dist))),
+        (SVPoint_2D(b_pt, T), (_clamped_frac(distance(b_pt, a1, T), a_dist), b_pt_β))
     )
 
 #= If lines defined by (a1, a2) and (b1, b2) are just touching at one of those endpoints and
@@ -363,25 +363,25 @@ to avoid floating point errors. If the points are not equal, we know that the hi
 take place at an endpoint and the fractions must be between 0 or 1 (exclusive). =#
 function _find_hinge_intersection(::Type{T}, a1, a2, b1, b2, a1_orient, a2_orient, b1_orient) where T
     pt, α, β = if equals(a1, b1)
-        _sv_point(a1, T), zero(T), zero(T)
+        SVPoint_2D(a1, T), zero(T), zero(T)
     elseif equals(a1, b2)
-        _sv_point(a1, T), zero(T), one(T)
+        SVPoint_2D(a1, T), zero(T), one(T)
     elseif equals(a2, b1)
-        _sv_point(a2, T), one(T), zero(T)
+        SVPoint_2D(a2, T), one(T), zero(T)
     elseif equals(a2, b2)
-        _sv_point(a2, T), one(T), one(T)
+        SVPoint_2D(a2, T), one(T), one(T)
     elseif a1_orient == 0
         β_val = _clamped_frac(distance(b1, a1, T), distance(b1, b2, T), eps(T))
-        _sv_point(a1, T), zero(T), β_val
+        SVPoint_2D(a1, T), zero(T), β_val
     elseif a2_orient == 0
         β_val = _clamped_frac(distance(b1, a2, T), distance(b1, b2, T), eps(T))
-        _sv_point(a2, T), one(T), β_val
+        SVPoint_2D(a2, T), one(T), β_val
     elseif b1_orient == 0
         α_val = _clamped_frac(distance(a1, b1, T), distance(a1, a2, T), eps(T))
-        _sv_point(b1, T), α_val, zero(T)
+        SVPoint_2D(b1, T), α_val, zero(T)
     else  # b2_orient == 0
         α_val = _clamped_frac(distance(a1, b2, T), distance(a1, a2, T), eps(T))
-        _sv_point(b2, T), α_val, one(T)
+        SVPoint_2D(b2, T), α_val, one(T)
     end
     return pt, (α, β)
 end
@@ -397,10 +397,10 @@ Regardless of point value, we know that it does not actually occur at an endpoin
 fractions must be between 0 or 1 (exclusive). =#
 function _find_cross_intersection(::Type{T}, a1, a2, b1, b2, a_ext, b_ext) where T
     # First line runs from a to a + Δa
-    (a1x, a1y), (a2x, a2y) = _tuple_point(a1, T), _tuple_point(a2, T)
+    (a1x, a1y), (a2x, a2y) = TuplePoint_2D(a1, T), TuplePoint_2D(a2, T)
     Δax, Δay = a2x - a1x, a2y - a1y
     # Second line runs from b to b + Δb 
-    (b1x, b1y), (b2x, b2y) = _tuple_point(b1, T), _tuple_point(b2, T)
+    (b1x, b1y), (b2x, b2y) = TuplePoint_2D(b1, T), TuplePoint_2D(b2, T)
     Δbx, Δby = b2x - b1x, b2y - b1y
     # Differences between starting points
     Δbax = b1x - a1x
@@ -430,7 +430,7 @@ function _find_cross_intersection(::Type{T}, a1, a2, b1, b2, a_ext, b_ext) where
     else
         (a1y + α * Δay + b1y + β * Δby) / 2
     end
-    pt = _sv_point((x, y), T)
+    pt = SVPoint_2D((x, y), T)
     # Check if point is within segment envelopes and adjust to endpoint if not
     if !_point_in_extent(pt, a_ext) || !_point_in_extent(pt, b_ext)
         pt, α, β = _nearest_endpoint(T, a1, a2, b1, b2)
@@ -465,7 +465,7 @@ function _nearest_endpoint(::Type{T}, a1, a2, b1, b2) where T
         α, β = _clamped_frac(distance(min_pt, a2, T), a_dist, eps(T)), one(T) - eps(T)
     end
     # Return point with smallest distance
-    return _sv_point(min_pt, T), α, β
+    return SVPoint_2D(min_pt, T), α, β
 end
 
 # Return value of x/y clamped between ϵ and 1 - ϵ
diff --git a/src/methods/clipping/predicates.jl b/src/methods/clipping/predicates.jl
index 8da6001b1..4b3d8326a 100644
--- a/src/methods/clipping/predicates.jl
+++ b/src/methods/clipping/predicates.jl
@@ -2,7 +2,7 @@ module Predicates
     using ExactPredicates, ExactPredicates.Codegen
     import ExactPredicates: ext
     import ExactPredicates.Codegen: group!, @genpredicate
-    import GeometryOps: _False, _True, _booltype, _tuple_point
+    import GeometryOps: _False, _True, _booltype, TuplePoint_2D
     import GeoInterface as GI
 
     #= Determine the orientation of c with regards to the oriented segment (a, b).
@@ -12,7 +12,7 @@ module Predicates
     orient(a, b, c; exact) = _orient(_booltype(exact), a, b, c)
     
     # If `exact` is `true`, use `ExactPredicates` to calculate the orientation.
-    _orient(::_True, a, b, c) = ExactPredicates.orient(_tuple_point(a, Float64), _tuple_point(b, Float64), _tuple_point(c, Float64))
+    _orient(::_True, a, b, c) = ExactPredicates.orient(TuplePoint_2D(a, Float64), TuplePoint_2D(b, Float64), TuplePoint_2D(c, Float64))
     # If `exact` is `false`, calculate the orientation without using `ExactPredicates`.
     function _orient(exact::_False, a, b, c)
         a = a .- c
@@ -29,7 +29,7 @@ module Predicates
     #= If `exact` is `true`, use exact cross product calculation created using
     `ExactPredicates`generated predicate. Note that as of now `ExactPredicates` requires
     Float64 so we must convert points a and b. =#
-    _cross(::_True, a, b) = _cross_exact(_tuple_point(a, Float64), _tuple_point(b, Float64))
+    _cross(::_True, a, b) = _cross_exact(TuplePoint_2D(a, Float64), TuplePoint_2D(b, Float64))
 
     # Exact cross product calculation using `ExactPredicates`.
     @genpredicate function _cross_exact(a :: 2, b :: 2)
diff --git a/src/methods/geom_relations/geom_geom_processors.jl b/src/methods/geom_relations/geom_geom_processors.jl
index 9df829974..23bdafd84 100644
--- a/src/methods/geom_relations/geom_geom_processors.jl
+++ b/src/methods/geom_relations/geom_geom_processors.jl
@@ -139,14 +139,14 @@ function _inner_line_curve_process(
     closed_line |= first_last_equal_line
     closed_curve |= first_last_equal_curve
     # Loop over each line segment
-    l_start = _tuple_point(GI.getpoint(line, closed_line ? nl : 1))
+    l_start = TuplePoint_2D(GI.getpoint(line, closed_line ? nl : 1))
     i = closed_line ? 1 : 2
     while i ≤ nl
-        l_end = _tuple_point(GI.getpoint(line, i))
-        c_start = _tuple_point(GI.getpoint(curve, closed_curve ? nc : 1))
+        l_end = TuplePoint_2D(GI.getpoint(line, i))
+        c_start = TuplePoint_2D(GI.getpoint(curve, closed_curve ? nc : 1))
         # Loop over each curve segment
         for j in (closed_curve ? 1 : 2):nc
-            c_end = _tuple_point(GI.getpoint(curve, j))
+            c_end = TuplePoint_2D(GI.getpoint(curve, j))
             # Check if line and curve segments meet
             seg_val, intr1, _ = _intersection_point(Float64, (l_start, l_end), (c_start, c_end); exact)
             # If segments are co-linear
@@ -228,12 +228,12 @@ end
 function _find_hinge_next_segments(α, β, ls, le, cs, ce, i, line, j, curve) 
     next_seg = if β == 1
         if α == 1  # hinge at endpoints, so next segment of both is needed
-            ((le, _tuple_point(GI.getpoint(line, i + 1))), (ce, _tuple_point(GI.getpoint(curve, j + 1))))
+            ((le, TuplePoint_2D(GI.getpoint(line, i + 1))), (ce, TuplePoint_2D(GI.getpoint(curve, j + 1))))
         else  # hinge at curve endpoint and line interior point, curve next segment needed 
-            ((ls, le), (ce, _tuple_point(GI.getpoint(curve, j + 1))))
+            ((ls, le), (ce, TuplePoint_2D(GI.getpoint(curve, j + 1))))
         end
     else  # hinge at curve interior point and line endpoint, line next segment needed
-        ((le, _tuple_point(GI.getpoint(line, i + 1))), (cs, ce))
+        ((le, TuplePoint_2D(GI.getpoint(line, i + 1))), (cs, ce))
     end
     return next_seg
 end
@@ -559,7 +559,7 @@ function _line_filled_curve_interactions(
     closed_line |= first_last_equal_line
 
     # See if first point is in an acceptable orientation
-    l_start = _tuple_point(GI.getpoint(line, closed_line ? nl : 1))
+    l_start = TuplePoint_2D(GI.getpoint(line, closed_line ? nl : 1))
     point_val = _point_filled_curve_orientation(l_start, curve; exact)
     if point_val == point_in
         in_curve = true
@@ -571,13 +571,13 @@ function _line_filled_curve_interactions(
 
     # Check for any intersections between line and curve
     for i in (closed_line ? 1 : 2):nl
-        l_end = _tuple_point(GI.getpoint(line, i))
-        c_start = _tuple_point(GI.getpoint(curve, nc))
+        l_end = TuplePoint_2D(GI.getpoint(line, i))
+        c_start = TuplePoint_2D(GI.getpoint(curve, nc))
         # If already interacted with all regions of curve, can stop
         in_curve && on_curve && out_curve && break
         # Check next segment of line against curve
         for j in 1:nc
-            c_end = _tuple_point(GI.getpoint(curve, j))
+            c_end = TuplePoint_2D(GI.getpoint(curve, j))
             # Check if two line and curve segments meet
             seg_val, _, _ = _intersection_point(Float64, (l_start, l_end), (c_start, c_end); exact)
             if seg_val != line_out
@@ -612,9 +612,9 @@ function _line_filled_curve_interactions(
                             x -> _euclid_distance(Float64, x, l_start)
                         end
                         sort!(ipoints, by = dist_from_lstart)
-                        p_start = _tuple_point(l_start)
+                        p_start = TuplePoint_2D(l_start)
                         for i in 1:(npoints + 1)
-                            p_end = i ≤ npoints ? _tuple_point(ipoints[i]) : l_end
+                            p_end = i ≤ npoints ? TuplePoint_2D(ipoints[i]) : l_end
                             mid_val = _point_filled_curve_orientation((p_start .+ p_end) ./ 2, curve; exact)
                             if mid_val == point_in
                                 in_curve = true
diff --git a/src/primitives.jl b/src/primitives.jl
index 44c81cc69..01de73cf8 100644
--- a/src/primitives.jl
+++ b/src/primitives.jl
@@ -52,32 +52,6 @@ TraitTarget
 ## Implementation
 
 =#
-
-const TuplePoint{T} = Tuple{T, T} where T <: AbstractFloat
-const TupleEdge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T
-
-struct SVPoint{N,T,Z,M} <: GeometryBasics.StaticArraysCore.StaticVector{N,T}
-    vals::NTuple{N,T}  # TODO: Should Z and M be booleans or BoolsAsTypes?
-end
-
-Base.getindex(p::SVPoint, i::Int64) = p.vals[i]
-
-SVPoint_2D(vals::NTuple{2,T}) where T = SVPoint{2, T, _False, _False}(vals)
-SVPoint_2D(vals) = SVPoint_2D((GI.x(vals), GI.y(vals)))
-
-SVPoint_3D(vals::NTuple{3,T}) where T = SVPoint{3, T, _True, _False}(vals)
-SVPoint_3D(vals) = SVPoint_3D((GI.x(vals), GI.y(vals), GI.z(vals)))
-
-SVPoint_4D(vals::NTuple{4,T}) where T = SVPoint{4, T, _True, _True}(vals)
-SVPoint_4D(vals) = SVPoint_4D((GI.x(vals), GI.y(vals), GI.z(vals), GI.m(vals)))
-
-SVPoint(vals::NTuple{2, <:AbstractFloat}) = SVPoint_2D(vals)
-SVPoint(vals::NTuple{3, <:AbstractFloat}) = SVPoint_3D(vals)
-SVPoint(vals::NTuple{4, <:AbstractFloat}) = SVPoint_4D(vals)
-
-const SVEdge{T} = Tuple{SVPoint{N,T,Z,M}, SVPoint{N,T,Z,M}} where {N,T,Z,M}
-const Edge{T} = Union{TupleEdge{T}, SVEdge{T}} where T
-
 const THREADED_KEYWORD = "- `threaded`: `true` or `false`. Whether to use multithreading. Defaults to `false`."
 const CRS_KEYWORD = "- `crs`: The CRS to attach to geometries. Defaults to `nothing`."
 const CALC_EXTENT_KEYWORD = "- `calc_extent`: `true` or `false`. Whether to calculate the extent. Defaults to `false`."
diff --git a/src/transformations/segmentize.jl b/src/transformations/segmentize.jl
index 2a6d09828..1473fb982 100644
--- a/src/transformations/segmentize.jl
+++ b/src/transformations/segmentize.jl
@@ -187,7 +187,7 @@ function _segmentize(::Type{T}, method::Union{LinearSegments, GeodesicSegments},
     x1, y1 = GI.x(first_coord), GI.y(first_coord)
     new_coords = _get_point_type(T)[] #NTuple{2, Float64}[]
     sizehint!(new_coords, GI.npoint(geom))
-    push!(new_coords, _sv_point((x1, y1), T))
+    push!(new_coords, SVPoint_2D((x1, y1), T))
     for coord in Iterators.drop(GI.getpoint(geom), 1)
         x2, y2 = GI.x(coord), GI.y(coord)
         _fill_linear_kernel!(T, method, new_coords, x1, y1, x2, y2)
@@ -203,12 +203,12 @@ function _fill_linear_kernel!(::Type{T}, method::LinearSegments, new_coords::Vec
         n_segments = ceil(Int, distance / method.max_distance)
         for i in 1:(n_segments - 1)
             t = i / n_segments
-            push!(new_coords, _sv_point((x1 + t * dx, y1 + t * dy), T))
+            push!(new_coords, SVPoint_2D((x1 + t * dx, y1 + t * dy), T))
         end
     end
     # End the line with the original coordinate,
     # to avoid any multiplication errors.
-    push!(new_coords, _sv_point((x2, y2), T))
+    push!(new_coords, SVPoint_2D((x2, y2), T))
     return nothing
 end
 #=
diff --git a/src/transformations/simplify.jl b/src/transformations/simplify.jl
index 65848e644..0263e90b3 100644
--- a/src/transformations/simplify.jl
+++ b/src/transformations/simplify.jl
@@ -508,14 +508,6 @@ function _build_tolerances(f, points)
     return real_tolerances
 end
 
-function _sv_points(::Type{T}, geom) where T
-    points = Array{_get_point_type(T)}(undef, GI.npoint(geom))
-    for (i, p) in enumerate(GI.getpoint(geom))
-        points[i] = _sv_point(p, T)
-    end
-    return points
-end
-
 function _get_points(alg, points, tolerances)
     ## This assumes that `alg` has the properties
     ## `tol`, `number`, and `ratio` available...
diff --git a/src/types.jl b/src/types.jl
index 4e6a0dec1..75a97d787 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -68,6 +68,55 @@ struct _False <: BoolsAsTypes end
 @inline _booltype(x::Bool)::BoolsAsTypes = x ? _True() : _False()
 @inline _booltype(x::BoolsAsTypes)::BoolsAsTypes = x
 
+
+const TuplePoint{N, T} = Tuple{T, T} where T <: AbstractFloat
+const TupleEdge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T
+
+# Should we just have default Float64??
+TuplePoint_2D(vals) = (GI.x(vals), GI.y(vals))
+TuplePoint_2D(vals, ::Type{T}) where T <: AbstractFloat = T.((GI.x(vals), GI.y(vals)))
+
+TuplePoint_3D(vals) = (GI.x(vals), GI.y(vals), GI.z(vals))
+TuplePoint_3D(vals, ::Type{T}) where T <: AbstractFloat = T.((GI.x(vals), GI.y(vals), GI.z(vals)))
+
+TuplePoint_4D(vals) = (GI.x(vals), GI.y(vals), GI.z(vals), GI.m(vals))
+TuplePoint_4D(vals, ::Type{T}) where T <: AbstractFloat = T.((GI.x(vals), GI.y(vals), GI.z(vals), GI.m(vals)))
+
+#=
+## `SVPoint`
+
+
+=#
+struct SVPoint{N, T, Z <: BoolsAsTypes, M <: BoolsAsTypes} <: GeometryBasics.StaticArraysCore.StaticVector{N,T}
+    vals::NTuple{N,T}  # TODO: Should Z and M be booleans or BoolsAsTypes?
+end
+
+Base.getindex(p::SVPoint, i::Int64) = p.vals[i]
+
+# Should we just have default Float64??
+SVPoint_2D(vals) = _SVPoint_2D(TuplePoint_2D(vals))
+SVPoint_2D(vals, ::Type{T}) where T <: AbstractFloat = _SVPoint_2D(TuplePoint_2D(vals, T))
+_SVPoint_2D(vals::NTuple{2,T}) where T = SVPoint{2, T, _False, _False}(vals)
+
+SVPoint_3D(vals) = _SVPoint_3D(TuplePoint_3D(vals))
+SVPoint_3D(vals, ::Type{T}) where T <: AbstractFloat = _SVPoint_3D(TuplePoint_3D(vals, T))
+_SVPoint_3D(vals::NTuple{3,T}) where T = SVPoint{3, T, _True, _False}(vals)
+
+SVPoint_4D(vals) = _SVPoint_4D(TuplePoint_4D(vals))
+SVPoint_4D(vals, ::Type{T}) where T <: AbstractFloat = _SVPoint_4D(TuplePoint_4D(vals, T))
+_SVPoint_4D(vals::NTuple{4,T}) where T = SVPoint{4, T, _True, _True}(vals)
+
+const SVEdge{T} = Tuple{SVPoint{N,T,Z,M}, SVPoint{N,T,Z,M}} where {N,T,Z,M}
+
+#=
+Get type of points and polygons made through library functionality (e.g. clipping)
+TODO: Increase type options as library expands capabilities
+=#
+_get_point_type(::Type{T}) where T = SVPoint{2, T, _False, _False}
+_get_poly_type(::Type{T}) where T =
+    GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{_get_point_type(T)}, Nothing, Nothing}}, Nothing, Nothing}
+
+const Edge{T} = Union{TupleEdge{T}, SVEdge{T}} where T
 #=
 
 ## `GEOS`
diff --git a/src/utils.jl b/src/utils.jl
index d7826eff5..e9154f583 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -54,6 +54,10 @@ function polygon_to_line(poly)
     return GI.LineString(collect(GI.getgeom(GI.getgeom(poly, 1))))
 end
 
+#= TODO: Should the `_to_edges` and `_to_points` functions return tuples or svpoints and
+should they be dimensionally specific? If the dimension is found with if/else like in
+the transform functions, is that going to effect performance negatively? Note the existance
+of the _sv_points function that was within `simplify.jl`=#
 
 """
     to_edges()
@@ -81,10 +85,10 @@ function _to_edges!(edges::Vector, ::GI.AbstractGeometryTrait, fc, n)
 end
 function _to_edges!(edges::Vector, ::GI.AbstractCurveTrait, geom, n)
     p1 = GI.getpoint(geom, 1) 
-    p1x, p1y = GI.x(p1), GI.y(p1)
+    p1x, p1y = TuplePoint_2D(p1)
     for i in 2:GI.npoint(geom)
         p2 = GI.getpoint(geom, i)
-        p2x, p2y = GI.x(p2), GI.y(p2)
+        p2x, p2y = TuplePoint_2D(p2)
         edges[n] = (p1x, p1y), (p2x, p2y)
         p1x, p1y = p2x, p2y
         n += 1
@@ -92,9 +96,6 @@ function _to_edges!(edges::Vector, ::GI.AbstractCurveTrait, geom, n)
     return n
 end
 
-_tuple_point(p) = GI.x(p), GI.y(p)
-_tuple_point(p, ::Type{T}) where T = T(GI.x(p)), T(GI.y(p))
-
 function to_extent(edges::Vector{<:Edge})
     x, y = extrema(first, edges)
     Extents.Extent(X=x, Y=y)
@@ -122,21 +123,20 @@ function _to_points!(points::Vector, ::Union{AbstractCurveTrait,MultiPointTrait}
     n = 0
     for p in GI.getpoint(geom)
         n += 1
-        points[n] = _tuple_point(p)
+        points[n] = TuplePoint_2D(p)
     end
     return n
 end
 
+function _sv_points(::Type{T}, geom) where T
+    points = Array{_get_point_type(T)}(undef, GI.npoint(geom))
+    for (i, p) in enumerate(GI.getpoint(geom))
+        points[i] = SVPoint_2D(p, T)
+    end
+    return points
+end
+
 function _point_in_extent(p, extent::Extents.Extent)
     (x1, x2), (y1, y2) = extent.X, extent.Y
     return x1 ≤ GI.x(p) ≤ x2 && y1 ≤ GI.y(p) ≤ y2
 end
-
-_get_point_type(::Type{T}) where T = SVPoint{2, T, _False, _False}
-
-_sv_point(p, ::Type{T}) where T = SVPoint_2D(_tuple_point(p, T))
-# 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{_get_point_type(T)}, Nothing, Nothing}}, Nothing, Nothing}
-

From ebf83be5dd99368797e19526351529bbaddb9546 Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Wed, 12 Jun 2024 15:52:26 -0400
Subject: [PATCH 13/16] Change comment labels

---
 src/types.jl | 2 +-
 src/utils.jl | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/types.jl b/src/types.jl
index 75a97d787..a6b97b8ab 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -88,7 +88,7 @@ TuplePoint_4D(vals, ::Type{T}) where T <: AbstractFloat = T.((GI.x(vals), GI.y(v
 
 =#
 struct SVPoint{N, T, Z <: BoolsAsTypes, M <: BoolsAsTypes} <: GeometryBasics.StaticArraysCore.StaticVector{N,T}
-    vals::NTuple{N,T}  # TODO: Should Z and M be booleans or BoolsAsTypes?
+    vals::NTuple{N,T}
 end
 
 Base.getindex(p::SVPoint, i::Int64) = p.vals[i]
diff --git a/src/utils.jl b/src/utils.jl
index e9154f583..c4904b001 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -54,7 +54,7 @@ function polygon_to_line(poly)
     return GI.LineString(collect(GI.getgeom(GI.getgeom(poly, 1))))
 end
 
-#= TODO: Should the `_to_edges` and `_to_points` functions return tuples or svpoints and
+#= Question: Should the `_to_edges` and `_to_points` functions return tuples or svpoints and
 should they be dimensionally specific? If the dimension is found with if/else like in
 the transform functions, is that going to effect performance negatively? Note the existance
 of the _sv_points function that was within `simplify.jl`=#

From a20b8b5010f8d78d1ae3cc59f202ea23d5316201 Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Thu, 13 Jun 2024 13:41:53 -0400
Subject: [PATCH 14/16] Add general point constructor

---
 src/methods/clipping/clipping_processor.jl |  2 +-
 src/types.jl                               | 32 +++++++++++++---------
 2 files changed, 20 insertions(+), 14 deletions(-)

diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 7a8ff31e4..0f8bfef06 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/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::SVPoint{2, T, _False, _False}
+    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
diff --git a/src/types.jl b/src/types.jl
index a6b97b8ab..ac87b549e 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -87,24 +87,30 @@ TuplePoint_4D(vals, ::Type{T}) where T <: AbstractFloat = T.((GI.x(vals), GI.y(v
 
 
 =#
-struct SVPoint{N, T, Z <: BoolsAsTypes, M <: BoolsAsTypes} <: GeometryBasics.StaticArraysCore.StaticVector{N,T}
+struct SVPoint{N, T, Z, M} <: GeometryBasics.StaticArraysCore.StaticVector{N,T}
     vals::NTuple{N,T}
 end
-
 Base.getindex(p::SVPoint, i::Int64) = p.vals[i]
 
-# Should we just have default Float64??
-SVPoint_2D(vals) = _SVPoint_2D(TuplePoint_2D(vals))
-SVPoint_2D(vals, ::Type{T}) where T <: AbstractFloat = _SVPoint_2D(TuplePoint_2D(vals, T))
-_SVPoint_2D(vals::NTuple{2,T}) where T = SVPoint{2, T, _False, _False}(vals)
+# Syntactic sugar for type stability within functions with known point types
+SVPoint_2D(vals, ::Type{T} = Float64) where T <: AbstractFloat = _SVPoint_2D(TuplePoint_2D(vals, T))
+_SVPoint_2D(vals::NTuple{2,T}) where T = SVPoint{2, T, false, false}(vals)
+
+SVPoint_3D(vals, ::Type{T} = Float64) where T <: AbstractFloat = _SVPoint_3D(TuplePoint_3D(vals, T))
+_SVPoint_3D(vals::NTuple{3,T}) where T = SVPoint{3, T, true, false}(vals)
 
-SVPoint_3D(vals) = _SVPoint_3D(TuplePoint_3D(vals))
-SVPoint_3D(vals, ::Type{T}) where T <: AbstractFloat = _SVPoint_3D(TuplePoint_3D(vals, T))
-_SVPoint_3D(vals::NTuple{3,T}) where T = SVPoint{3, T, _True, _False}(vals)
+SVPoint_4D(vals, ::Type{T} = Float64) where T <: AbstractFloat = _SVPoint_4D(TuplePoint_4D(vals, T))
+_SVPoint_4D(vals::NTuple{4,T}) where T = SVPoint{4, T, true, true}(vals)
 
-SVPoint_4D(vals) = _SVPoint_4D(TuplePoint_4D(vals))
-SVPoint_4D(vals, ::Type{T}) where T <: AbstractFloat = _SVPoint_4D(TuplePoint_4D(vals, T))
-_SVPoint_4D(vals::NTuple{4,T}) where T = SVPoint{4, T, _True, _True}(vals)
+# General constructor when point type/size isn't known
+SVPoint(geom) = _SVPoint(GI.trait(geom), geom)
+# Make sure geometry is a point type
+_SVPoint(::GI.PointTrait, geom) = _SVPoint(tuples(geom))
+_SVPoint(trait::GI.AbstractTrait, _) = throw(ArgumentError("Geometry with trait $trait cannot be made into a point."))
+# Dispatch off of NTuple length to make point of needed dimension
+_SVPoint(geom::NTuple{2, T}) where T = _SVPoint_2D(geom)
+_SVPoint(geom::NTuple{3, T}) where T = _SVPoint_3D(geom)
+_SVPoint(geom::NTuple{4, T}) where T = _SVPoint_4D(geom)
 
 const SVEdge{T} = Tuple{SVPoint{N,T,Z,M}, SVPoint{N,T,Z,M}} where {N,T,Z,M}
 
@@ -112,7 +118,7 @@ const SVEdge{T} = Tuple{SVPoint{N,T,Z,M}, SVPoint{N,T,Z,M}} where {N,T,Z,M}
 Get type of points and polygons made through library functionality (e.g. clipping)
 TODO: Increase type options as library expands capabilities
 =#
-_get_point_type(::Type{T}) where T = SVPoint{2, T, _False, _False}
+_get_point_type(::Type{T}) where T = SVPoint{2, T, false, false}
 _get_poly_type(::Type{T}) where T =
     GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{_get_point_type(T)}, Nothing, Nothing}}, Nothing, Nothing}
 

From 5bddba73db0c02bf6a91a75678daa7260a393c24 Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Tue, 18 Jun 2024 15:09:53 -0400
Subject: [PATCH 15/16] Update point constructors

---
 src/types.jl | 17 ++++++++---------
 src/utils.jl | 12 ++++++++++++
 2 files changed, 20 insertions(+), 9 deletions(-)

diff --git a/src/types.jl b/src/types.jl
index ac87b549e..5775ad6f5 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -69,18 +69,13 @@ struct _False <: BoolsAsTypes end
 @inline _booltype(x::BoolsAsTypes)::BoolsAsTypes = x
 
 
-const TuplePoint{N, T} = Tuple{T, T} where T <: AbstractFloat
+const TuplePoint{T} = Tuple{T, T} where T <: AbstractFloat
 const TupleEdge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T
 
-# Should we just have default Float64??
-TuplePoint_2D(vals) = (GI.x(vals), GI.y(vals))
-TuplePoint_2D(vals, ::Type{T}) where T <: AbstractFloat = T.((GI.x(vals), GI.y(vals)))
 
-TuplePoint_3D(vals) = (GI.x(vals), GI.y(vals), GI.z(vals))
-TuplePoint_3D(vals, ::Type{T}) where T <: AbstractFloat = T.((GI.x(vals), GI.y(vals), GI.z(vals)))
-
-TuplePoint_4D(vals) = (GI.x(vals), GI.y(vals), GI.z(vals), GI.m(vals))
-TuplePoint_4D(vals, ::Type{T}) where T <: AbstractFloat = T.((GI.x(vals), GI.y(vals), GI.z(vals), GI.m(vals)))
+TuplePoint_2D(vals, ::Type{T} = Float64) where T <: AbstractFloat = T.((GI.x(vals), GI.y(vals)))
+TuplePoint_3D(vals, ::Type{T} = Float64) where T <: AbstractFloat = T.((GI.x(vals), GI.y(vals), GI.z(vals)))
+TuplePoint_4D(vals, ::Type{T} = Float64) where T <: AbstractFloat = T.((GI.x(vals), GI.y(vals), GI.z(vals), GI.m(vals)))
 
 #=
 ## `SVPoint`
@@ -92,6 +87,10 @@ struct SVPoint{N, T, Z, M} <: GeometryBasics.StaticArraysCore.StaticVector{N,T}
 end
 Base.getindex(p::SVPoint, i::Int64) = p.vals[i]
 
+const Point_2D = SVPoint{2, T, false, false}
+const Point_3D = SVPoint{3, T, true, false}
+const Point_4D = SVPoint{4, T, true, true}
+
 # Syntactic sugar for type stability within functions with known point types
 SVPoint_2D(vals, ::Type{T} = Float64) where T <: AbstractFloat = _SVPoint_2D(TuplePoint_2D(vals, T))
 _SVPoint_2D(vals::NTuple{2,T}) where T = SVPoint{2, T, false, false}(vals)
diff --git a/src/utils.jl b/src/utils.jl
index c4904b001..db9340a6b 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -66,11 +66,23 @@ Convert any geometry or collection of geometries into a flat
 vector of `Tuple{Tuple{Float64,Float64},Tuple{Float64,Float64}}` edges.
 """
 function to_edges(x, ::Type{T} = Float64) where T
+    
     edges = Vector{TupleEdge{T}}(undef, _nedge(x))
     _to_edges!(edges, x, 1)
     return edges
 end
 
+function find_point_constructor(x, ::Type{T})
+    if _ismeasured(x)
+        SVPoint_4D
+    elseif _is3d
+        SVPoint_3D
+    else
+        SVPoint_2D
+    end
+
+end
+
 _to_edges!(edges::Vector, x, n) = _to_edges!(edges, trait(x), x, n)
 function _to_edges!(edges::Vector, ::GI.FeatureCollectionTrait, fc, n)
     for f in GI.getfeature(fc)

From 3624dc670d0b43c5d94db400422203e6895fde1f Mon Sep 17 00:00:00 2001
From: Skylar Gering <sgering@g.hmc.edu>
Date: Thu, 20 Jun 2024 14:27:18 -0400
Subject: [PATCH 16/16] Update new point types for measured

---
 src/methods/clipping/clipping_processor.jl |  2 +-
 src/methods/clipping/difference.jl         |  2 +-
 src/methods/clipping/intersection.jl       |  6 +--
 src/transformations/flip.jl                |  6 +--
 src/transformations/segmentize.jl          |  2 +-
 src/transformations/transform.jl           |  6 +--
 src/transformations/tuples.jl              |  6 +--
 src/types.jl                               | 59 +++++++++++-----------
 src/utils.jl                               | 25 +++++----
 9 files changed, 57 insertions(+), 57 deletions(-)

diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 0f8bfef06..d75bb1456 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -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
diff --git a/src/methods/clipping/difference.jl b/src/methods/clipping/difference.jl
index 4928c3917..96db1c2a0 100644
--- a/src/methods/clipping/difference.jl
+++ b/src/methods/clipping/difference.jl
@@ -136,7 +136,7 @@ function _difference(
     if !isnothing(fix_multipoly) # Fix multipoly_a to prevent returning an invalid multipolygon
         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, T; target))
diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl
index de397ca95..40f8556dd 100644
--- a/src/methods/clipping/intersection.jl
+++ b/src/methods/clipping/intersection.jl
@@ -120,7 +120,7 @@ function _intersection(
     if !isnothing(fix_multipoly) # Fix multipoly_b to prevent duplicated intersection regions
         multipoly_b = fix_multipoly(multipoly_b, T)
     end
-    polys = Vector{_get_poly_type(T)}()
+    polys = Vector{PolyType2D{T}}()
     for poly_b in GI.getpolygon(multipoly_b)
         append!(polys, intersection(poly_a, poly_b, T; target))
     end
@@ -153,7 +153,7 @@ function _intersection(
         multipoly_b = fix_multipoly(multipoly_b, T)
         fix_multipoly = nothing
     end
-    polys = Vector{_get_poly_type(T)}()
+    polys = Vector{PolyType2D{T}}()
     for poly_a in GI.getpolygon(multipoly_a)
         append!(polys, intersection(poly_a, multipoly_b, T; target, fix_multipoly))
     end
@@ -197,7 +197,7 @@ were possible given geometry extents or if none are found, return an empty list
 GI.Points. =#
 function _intersection_points(::Type{T}, ::GI.AbstractTrait, a, ::GI.AbstractTrait, b; exact = _False()) where T
     # Initialize an empty list of points
-    result = Vector{_get_point_type(T)}()
+    result = Vector{PointType2D{T}}()
     # Check if the geometries extents even overlap
     Extents.intersects(GI.extent(a), GI.extent(b)) || return result
     # Create a list of edges from the two input geometries
diff --git a/src/transformations/flip.jl b/src/transformations/flip.jl
index 76ab724ef..2d5a184d3 100644
--- a/src/transformations/flip.jl
+++ b/src/transformations/flip.jl
@@ -17,15 +17,15 @@ $APPLY_KEYWORDS
 function flip(geom; kw...)
     if _ismeasured(geom)
         return apply(PointTrait(), geom; kw...) do p
-            SVPoint_4D(StaticArrays.SVector{4}((GI.y(p), GI.x(p), GI.z(p), GI.m(p))))
+            SVPoint_4D((GI.y(p), GI.x(p), GI.z(p), GI.m(p)))
         end
     elseif _is3d(geom)
         return apply(PointTrait(), geom; kw...) do p
-            SVPoint_3D(StaticArrays.SVector{3}((GI.y(p), GI.x(p), GI.z(p))))
+            SVPoint_3D((GI.y(p), GI.x(p), GI.z(p)))
         end
     else
         return apply(PointTrait(), geom; kw...) do p
-            SVPoint_2D(StaticArrays.SVector{2}((GI.y(p), GI.x(p))))
+            SVPoint_2D((GI.y(p), GI.x(p)))
         end
     end
 end
diff --git a/src/transformations/segmentize.jl b/src/transformations/segmentize.jl
index 1473fb982..8da2bc5fd 100644
--- a/src/transformations/segmentize.jl
+++ b/src/transformations/segmentize.jl
@@ -185,7 +185,7 @@ and calls out to the "kernel" function which we've defined per linesegment.
 function _segmentize(::Type{T}, method::Union{LinearSegments, GeodesicSegments}, geom, ::Union{GI.LineStringTrait, GI.LinearRingTrait}) where T
     first_coord = GI.getpoint(geom, 1)
     x1, y1 = GI.x(first_coord), GI.y(first_coord)
-    new_coords = _get_point_type(T)[] #NTuple{2, Float64}[]
+    new_coords = PointType2D{T}[]
     sizehint!(new_coords, GI.npoint(geom))
     push!(new_coords, SVPoint_2D((x1, y1), T))
     for coord in Iterators.drop(GI.getpoint(geom), 1)
diff --git a/src/transformations/transform.jl b/src/transformations/transform.jl
index 47005e352..7faa8dc54 100644
--- a/src/transformations/transform.jl
+++ b/src/transformations/transform.jl
@@ -41,15 +41,15 @@ GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearR
 function transform(f, geom, ::Type{T} = Float64; kw...) where T
     if _ismeasured(geom)
         return apply(PointTrait(), geom; kw...) do p
-            SVPoint_4D(f(StaticArrays.SVector{4}((GI.x(p), GI.y(p), GI.z(p), GI.m(p)))))
+            SVPoint_4D(f(SVPoint_4D((GI.x(p), GI.y(p), GI.z(p), GI.m(p)))))
         end
     elseif _is3d(geom)
         return apply(PointTrait(), geom; kw...) do p
-            SVPoint_3D(f(StaticArrays.SVector{3}((GI.x(p), GI.y(p), GI.z(p)))))
+            SVPoint_3D(f(SVPoint_3D((GI.x(p), GI.y(p), GI.z(p)))))
         end
     else
         return apply(PointTrait(), geom; kw...) do p
-            SVPoint_2D(f(StaticArrays.SVector{2}((GI.x(p), GI.y(p)))))
+            SVPoint_2D(f(SVPoint_2D((GI.x(p), GI.y(p)))))
         end
     end
 end
diff --git a/src/transformations/tuples.jl b/src/transformations/tuples.jl
index 7ac511857..c3c6bf2ed 100644
--- a/src/transformations/tuples.jl
+++ b/src/transformations/tuples.jl
@@ -15,15 +15,15 @@ $APPLY_KEYWORDS
 function tuples(geom, ::Type{T} = Float64; kw...) where T
     if _ismeasured(geom)
         return apply(PointTrait(), geom; kw...) do p
-            (T(GI.x(p)), T(GI.y(p)), T(GI.z(p)), T(GI.m(p)))
+            TuplePoint_4D(p, T)
         end
     elseif _is3d(geom)
         return apply(PointTrait(), geom; kw...) do p
-            (T(GI.x(p)), T(GI.y(p)), T(GI.z(p)))
+            TuplePoint_3D(p, T)
         end
     else
         return apply(PointTrait(), geom; kw...) do p
-            (T(GI.x(p)), T(GI.y(p)))
+            TuplePoint_2D(p, T)
         end
     end
 end
diff --git a/src/types.jl b/src/types.jl
index 5775ad6f5..e2bc65d60 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -69,13 +69,20 @@ struct _False <: BoolsAsTypes end
 @inline _booltype(x::BoolsAsTypes)::BoolsAsTypes = x
 
 
-const TuplePoint{T} = Tuple{T, T} where T <: AbstractFloat
-const TupleEdge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T
+const TuplePoint{N, T} = NTuple{N, T} where {N, T}
+const TupleEdge{N, T} = Tuple{TuplePoint{N, T}, TuplePoint{N, T}} where {N, T}
 
+TuplePoint(geom, ::Type{T} = Float64) where T = _TuplePoint(T, GI.trait(geom), geom)
+_TuplePoint(::Type{T}, ::GI.PointTrait, geom) where T = tuples(geom, T)
+_TuplePoint(::Type{T}, trait::GI.AbstractTrait, _) where T = throw(ArgumentError("Geometry with trait $trait cannot be made into a point."))
 
-TuplePoint_2D(vals, ::Type{T} = Float64) where T <: AbstractFloat = T.((GI.x(vals), GI.y(vals)))
-TuplePoint_3D(vals, ::Type{T} = Float64) where T <: AbstractFloat = T.((GI.x(vals), GI.y(vals), GI.z(vals)))
-TuplePoint_4D(vals, ::Type{T} = Float64) where T <: AbstractFloat = T.((GI.x(vals), GI.y(vals), GI.z(vals), GI.m(vals)))
+TuplePoint_2D(vals, ::Type{T} = Float64) where T = TuplePoint{2, T}((GI.x(vals), GI.y(vals)))
+
+TuplePoint_3D(vals, ::Type{T} = Float64, M = _False()) where T = _TuplePoint_3D(T, _booltype(M), vals)
+_TuplePoint_3D(::Type{T}, ::_False, vals) where T = TuplePoint{3, T}((GI.x(vals), GI.y(vals), GI.z(vals)))
+_TuplePoint_3D(::Type{T}, ::_True, vals) where T = TuplePoint{3, T}((GI.x(vals), GI.y(vals), GI.m(vals)))
+
+TuplePoint_4D(vals, ::Type{T} = Float64) where T = TuplePoint{4, T}((GI.x(vals), GI.y(vals), GI.z(vals), GI.m(vals)))
 
 #=
 ## `SVPoint`
@@ -86,42 +93,36 @@ struct SVPoint{N, T, Z, M} <: GeometryBasics.StaticArraysCore.StaticVector{N,T}
     vals::NTuple{N,T}
 end
 Base.getindex(p::SVPoint, i::Int64) = p.vals[i]
+# TODO: overload `similar_type``
 
-const Point_2D = SVPoint{2, T, false, false}
-const Point_3D = SVPoint{3, T, true, false}
-const Point_4D = SVPoint{4, T, true, true}
+const SVEdge{N, T, Z, M} = Tuple{SVPoint{N,T,Z,M}, SVPoint{N,T,Z,M}} where {N,T,Z,M}
 
-# Syntactic sugar for type stability within functions with known point types
-SVPoint_2D(vals, ::Type{T} = Float64) where T <: AbstractFloat = _SVPoint_2D(TuplePoint_2D(vals, T))
-_SVPoint_2D(vals::NTuple{2,T}) where T = SVPoint{2, T, false, false}(vals)
+# General SVPoint constructor when point type/size isn't known
+SVPoint(geom, ::Type{T} = Float64) where T = _SVPoint(T, GI.trait(geom), geom)
+_SVPoint(::Type{T}, ::GI.PointTrait, geom) where T = svpoints(geom, T)
+_SVPoint(::Type{T}, trait::GI.AbstractTrait, _) where T = throw(ArgumentError("Geometry with trait $trait cannot be made into a point."))
 
-SVPoint_3D(vals, ::Type{T} = Float64) where T <: AbstractFloat = _SVPoint_3D(TuplePoint_3D(vals, T))
-_SVPoint_3D(vals::NTuple{3,T}) where T = SVPoint{3, T, true, false}(vals)
+# Syntactic sugar for type stability within functions with known point types
+const PointType2D{T} = SVPoint{2, T, false, false} where T
+const PointType3D{T} = SVPoint{3, T, true, false} where T
+const PointType3DM{T} = SVPoint{3, T, false, true} where T
+const PointType4D{T} = SVPoint{4, T, true, true} where T
 
-SVPoint_4D(vals, ::Type{T} = Float64) where T <: AbstractFloat = _SVPoint_4D(TuplePoint_4D(vals, T))
-_SVPoint_4D(vals::NTuple{4,T}) where T = SVPoint{4, T, true, true}(vals)
+SVPoint_2D(vals, ::Type{T} = Float64) where T <: AbstractFloat = PointType2D{T}(TuplePoint_2D(vals, T))
 
-# General constructor when point type/size isn't known
-SVPoint(geom) = _SVPoint(GI.trait(geom), geom)
-# Make sure geometry is a point type
-_SVPoint(::GI.PointTrait, geom) = _SVPoint(tuples(geom))
-_SVPoint(trait::GI.AbstractTrait, _) = throw(ArgumentError("Geometry with trait $trait cannot be made into a point."))
-# Dispatch off of NTuple length to make point of needed dimension
-_SVPoint(geom::NTuple{2, T}) where T = _SVPoint_2D(geom)
-_SVPoint(geom::NTuple{3, T}) where T = _SVPoint_3D(geom)
-_SVPoint(geom::NTuple{4, T}) where T = _SVPoint_4D(geom)
+SVPoint_3D(vals, ::Type{T} = Float64, M = _False()) where {T <: AbstractFloat} = _SVPoint_3D(T, _booltype(M), vals)
+_SVPoint_3D(T, M::_False, vals) = PointType3D{T}(TuplePoint_3D(vals, T, M))
+_SVPoint_3D(T, M::_True, vals) = PointType3DM{T}(TuplePoint_3D(vals, T, M))
 
-const SVEdge{T} = Tuple{SVPoint{N,T,Z,M}, SVPoint{N,T,Z,M}} where {N,T,Z,M}
+SVPoint_4D(vals, ::Type{T} = Float64) where T <: AbstractFloat = PointType4D{T}(TuplePoint_4D(vals, T))
 
 #=
 Get type of points and polygons made through library functionality (e.g. clipping)
 TODO: Increase type options as library expands capabilities
 =#
-_get_point_type(::Type{T}) where T = SVPoint{2, T, false, false}
-_get_poly_type(::Type{T}) where T =
-    GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{_get_point_type(T)}, Nothing, Nothing}}, Nothing, Nothing}
+const PolyType2D{T} = GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{PointType2D{T}}, Nothing, Nothing}}, Nothing, Nothing} where T
 
-const Edge{T} = Union{TupleEdge{T}, SVEdge{T}} where T
+const Edge{N, T} = Union{TupleEdge{N, T}, SVEdge{N, T, Z, M}} where {N, T, Z, M}
 #=
 
 ## `GEOS`
diff --git a/src/utils.jl b/src/utils.jl
index db9340a6b..b6534b818 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -65,23 +65,22 @@ of the _sv_points function that was within `simplify.jl`=#
 Convert any geometry or collection of geometries into a flat 
 vector of `Tuple{Tuple{Float64,Float64},Tuple{Float64,Float64}}` edges.
 """
-function to_edges(x, ::Type{T} = Float64) where T
-    
-    edges = Vector{TupleEdge{T}}(undef, _nedge(x))
+function to_edges(x, ::Type{T} = Float64, ::Val{N} = Val(2)) where {T, N}
+    edges = Vector{TupleEdge{N, T}}(undef, _nedge(x))
     _to_edges!(edges, x, 1)
     return edges
 end
 
-function find_point_constructor(x, ::Type{T})
-    if _ismeasured(x)
-        SVPoint_4D
-    elseif _is3d
-        SVPoint_3D
-    else
-        SVPoint_2D
-    end
+# function find_point_constructor(x, ::Type{T})
+#     if _ismeasured(x)
+#         SVPoint_4D
+#     elseif _is3d(x)
+#         SVPoint_3D
+#     else
+#         SVPoint_2D
+#     end
 
-end
+# end
 
 _to_edges!(edges::Vector, x, n) = _to_edges!(edges, trait(x), x, n)
 function _to_edges!(edges::Vector, ::GI.FeatureCollectionTrait, fc, n)
@@ -141,7 +140,7 @@ function _to_points!(points::Vector, ::Union{AbstractCurveTrait,MultiPointTrait}
 end
 
 function _sv_points(::Type{T}, geom) where T
-    points = Array{_get_point_type(T)}(undef, GI.npoint(geom))
+    points = Array{PointType2D{T}}(undef, GI.npoint(geom))
     for (i, p) in enumerate(GI.getpoint(geom))
         points[i] = SVPoint_2D(p, T)
     end