diff --git a/ext/GeometryOpsProjExt/segmentize.jl b/ext/GeometryOpsProjExt/segmentize.jl
index 308c97188..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!
+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))
@@ -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.SVPoint_2D((x, y), T))
         end
     end
     # End the line with the original coordinate,
     # to avoid any multiplication errors.
-    push!(new_coords, (x2, y2))
+    push!(new_coords, GeometryOps.SVPoint_2D((x2, y2), T))
     return nothing
 end
\ No newline at end of file
diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index 1aced0e82..788cde8df 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -14,9 +14,7 @@ using GeoInterface.Extents: Extents
 
 const GI = GeoInterface
 const GB = GeometryBasics
-
-const TuplePoint{T} = Tuple{T, T} where T <: AbstractFloat
-const Edge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T
+const SA = GeometryBasics.StaticArrays
 
 include("types.jl")
 include("primitives.jl")
@@ -55,6 +53,7 @@ include("transformations/flip.jl")
 include("transformations/reproject.jl")
 include("transformations/segmentize.jl")
 include("transformations/simplify.jl")
+include("transformations/svpoints.jl")
 include("transformations/tuples.jl")
 include("transformations/transform.jl")
 include("transformations/correction/geometry_correction.jl")
diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl
index 9dd0e4167..48a84876a 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 = TuplePoint_2D(GI.getpoint(geom, 1), T)
     # Loop over line segments of line string
     for point₂ in GI.getpoint(geom)
+        x2, y2 = TuplePoint_2D(point₂, T)
         # Calculate length of line segment
-        length_component = sqrt(
-            (GI.x(point₂) - GI.x(point₁))^2 +
-            (GI.y(point₂) - GI.y(point₁))^2
-        )
+        length_component = distance((x1, y1), (x2, y2), T)
         # Accumulate the line segment length into `length`
         length += length_component
         # Weighted average of line segment centroids
-        xcentroid += (GI.x(point₁) + GI.x(point₂)) * (length_component / 2)
-        ycentroid += (GI.y(point₁) + GI.y(point₂)) * (length_component / 2)
-        #centroid = centroid .+ ((point₁ .+ point₂) .* (length_component / 2))
+        xcentroid += (x1 + x2) * (length_component / 2)
+        ycentroid += (y1 + y2) * (length_component / 2)
         # Advance the point buffer by 1 point to move to next line segment
-        point₁ = point₂
+        x1, y1 = x2, y2
     end
     xcentroid /= length
     ycentroid /= length
-    return (xcentroid, ycentroid), length
+    return SVPoint_2D((xcentroid, ycentroid), T), length
 end
 
 """
@@ -109,9 +106,10 @@ Returns the centroid and area of a given geometry.
 function centroid_and_area(geom, ::Type{T}=Float64; threaded=false) where T
     target = TraitTarget{Union{GI.PolygonTrait,GI.LineStringTrait,GI.LinearRingTrait}}()
     init = (zero(T), zero(T)), zero(T)
-    applyreduce(_combine_centroid_and_area, target, geom; threaded, init) do g
+    centroid, area = applyreduce(_combine_centroid_and_area, target, geom; threaded, init) do g
         _centroid_and_area(GI.trait(g), g, T)
     end
+    return SVPoint_2D(centroid, T), area
 end
 
 function _centroid_and_area(
@@ -146,14 +144,14 @@ function _centroid_and_area(
 end
 function _centroid_and_area(::GI.PolygonTrait, geom, ::Type{T}) where T
     # Exterior ring's centroid and area
-    (xcentroid, ycentroid), area = centroid_and_area(GI.getexterior(geom), T)
+    (xcentroid, ycentroid), area = _centroid_and_area(GI.LinearRingTrait(), GI.getexterior(geom), T)
     # Weight exterior centroid by area
     xcentroid *= area
     ycentroid *= area
     # Loop over any holes within the polygon
     for hole in GI.gethole(geom)
         # Hole polygon's centroid and area
-        (xinterior, yinterior), interior_area = centroid_and_area(hole, T)
+        (xinterior, yinterior), interior_area = _centroid_and_area(GI.LinearRingTrait(), hole, T)
         # Accumulate the area component into `area`
         area -= interior_area
         # Weighted average of centroid components
diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index aaf938d95..d75bb1456 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::SVPoint{2, T, false, false}
     inter::Bool = false        # If ipt, true, else 0
     neighbor::Int = 0          # If ipt, index of equivalent point in a_list or b_list, else 0
     idx::Int = 0               # If crossing point, index within sorted a_idx_list
@@ -90,7 +90,7 @@ function _build_a_list(::Type{T}, poly_a, poly_b; exact) where T
     # Loop through points of poly_a
     local a_pt1
     for (i, a_p2) in enumerate(GI.getpoint(poly_a))
-        a_pt2 = (T(GI.x(a_p2)), T(GI.y(a_p2)))
+        a_pt2 = SVPoint_2D(a_p2, T)
         if i <= 1 || (a_pt1 == a_pt2)  # don't repeat points
             a_pt1 = a_pt2
             continue
@@ -103,7 +103,7 @@ function _build_a_list(::Type{T}, poly_a, poly_b; exact) where T
         local b_pt1
         prev_counter = a_count
         for (j, b_p2) in enumerate(GI.getpoint(poly_b))
-            b_pt2 = _tuple_point(b_p2, T)
+            b_pt2 = SVPoint_2D(b_p2, T)
             if j <= 1 || (b_pt1 == b_pt2)  # don't repeat points
                 b_pt1 = b_pt2
                 continue
@@ -194,7 +194,7 @@ function _build_b_list(::Type{T}, a_idx_list, a_list, n_b_intrs, poly_b) where T
     # Loop over points in poly_b and add each point and intersection point
     local b_pt1
     for (i, b_p2) in enumerate(GI.getpoint(poly_b))
-        b_pt2 = _tuple_point(b_p2, T)
+        b_pt2 = SVPoint_2D(b_p2, T)
         if i ≤ 1 || (b_pt1 == b_pt2)  # don't repeat points
             b_pt1 = b_pt2
             continue
@@ -543,7 +543,7 @@ function _trace_polynodes(::Type{T}, a_list, b_list, a_idx_list, f_step, poly_a,
     n_a_pts, n_b_pts = length(a_list), length(b_list)
     total_pts = n_a_pts + n_b_pts
     n_cross_pts = length(a_idx_list)
-    return_polys = Vector{_get_poly_type(T)}(undef, 0)
+    return_polys = Vector{PolyType2D{T}}(undef, 0)
     # Keep track of number of processed intersection points
     visited_pts = 0
     processed_pts = 0
@@ -601,11 +601,6 @@ function _trace_polynodes(::Type{T}, a_list, b_list, a_idx_list, f_step, poly_a,
     return return_polys
 end
 
-# Get type of polygons that will be made
-# TODO: Increase type options
-_get_poly_type(::Type{T}) where T =
-    GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{Tuple{T, T}}, Nothing, Nothing}}, Nothing, Nothing}
-
 #=
     _find_non_cross_orientation(a_list, b_list, a_poly, b_poly; exact)
 
@@ -642,12 +637,12 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol
     # Remove set of holes from all polygons
     for i in 1:n_polys
         n_new_per_poly = 0
-        for curr_hole in Iterators.map(tuples, hole_iterator) # loop through all holes
+        for curr_hole in Iterators.map(Base.Fix2(svpoints, T), hole_iterator) # loop through all holes
             # loop through all pieces of original polygon (new pieces added to end of list)
             for j in Iterators.flatten((i:i, (n_polys + 1):(n_polys + n_new_per_poly)))
                 curr_poly = return_polys[j]
                 remove_poly_idx[j] && continue
-                curr_poly_ext = GI.nhole(curr_poly) > 0 ? GI.Polygon(StaticArrays.SVector(GI.getexterior(curr_poly))) : curr_poly
+                curr_poly_ext = GI.nhole(curr_poly) > 0 ? GI.Polygon(SA.SVector(GI.getexterior(curr_poly))) : curr_poly
                 in_ext, on_ext, out_ext = _line_polygon_interactions(curr_hole, curr_poly_ext; exact, closed_line = true)
                 if in_ext  # hole is at least partially within the polygon's exterior
                     new_hole, new_hole_poly, n_new_pieces = _combine_holes!(T, curr_hole, curr_poly, return_polys, remove_hole_idx)
@@ -670,7 +665,7 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol
                         end
                     end
                 # polygon is completly within hole
-                elseif coveredby(curr_poly_ext, GI.Polygon(StaticArrays.SVector(curr_hole)))
+                elseif coveredby(curr_poly_ext, GI.Polygon(SA.SVector(curr_hole)))
                     remove_poly_idx[j] = true
                 end
             end
@@ -698,16 +693,16 @@ changes.
 function _combine_holes!(::Type{T}, new_hole, curr_poly, return_polys, remove_hole_idx) where T
     n_new_polys = 0
     empty!(remove_hole_idx)
-    new_hole_poly = GI.Polygon(StaticArrays.SVector(new_hole))
+    new_hole_poly = GI.Polygon(SA.SVector(new_hole))
     # Combine any existing holes in curr_poly with new hole
     for (k, old_hole) in enumerate(GI.gethole(curr_poly))
-        old_hole_poly = GI.Polygon(StaticArrays.SVector(old_hole))
+        old_hole_poly = GI.Polygon(SA.SVector(old_hole))
         if intersects(new_hole_poly, old_hole_poly)
             # If the holes intersect, combine them into a bigger hole
             hole_union = union(new_hole_poly, old_hole_poly, T; target = GI.PolygonTrait())[1]
             push!(remove_hole_idx, k + 1)
             new_hole = GI.getexterior(hole_union)
-            new_hole_poly = GI.Polygon(StaticArrays.SVector(new_hole))
+            new_hole_poly = GI.Polygon(SA.SVector(new_hole))
             n_pieces = GI.nhole(hole_union)
             if n_pieces > 0  # if the hole has a hole, then this is a new polygon piece! 
                 append!(return_polys, [GI.Polygon([h]) for h in GI.gethole(hole_union)])
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 1303ceb7e..8a4d158a7 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; exact) w
     n_intr_pts = length(intr_list)
     # If an impossible number of intersection points, return original polygon
     if n_intr_pts < 2 || isodd(n_intr_pts)
-        return [tuples(poly)]
+        return [svpoints(poly, T)]
     end
     # Cut polygon by line
     cut_coords = _cut(T, ext_poly, line, poly_list, intr_list, n_intr_pts; exact)
@@ -103,7 +103,7 @@ function _cut(::Type{T}, geom, line, geom_list, intr_list, n_intr_pts; exact) wh
     _flag_ent_exit!(GI.LineTrait(), line, geom_list; exact)
     # Add first point to output list
     return_coords = [[geom_list[1].point]]
-    cross_backs = [(T(Inf),T(Inf))]
+    cross_backs = [SVPoint_2D((Inf, Inf), T)]
     poly_idx = 1
     n_polys = 1
     # Walk around original polygon to find split polygons
diff --git a/src/methods/clipping/difference.jl b/src/methods/clipping/difference.jl
index 604e4606b..96db1c2a0 100644
--- a/src/methods/clipping/difference.jl
+++ b/src/methods/clipping/difference.jl
@@ -62,10 +62,10 @@ function _difference(
         # add case for if they polygons are the same (all intersection points!)
         # add a find_first check to find first non-inter poly!
         if b_in_a && !a_in_b  # b in a and can't be the same polygon
-            poly_a_b_hole = GI.Polygon([tuples(ext_a), tuples(ext_b)])
+            poly_a_b_hole = GI.Polygon([svpoints(ext_a, T), svpoints(ext_b, T)])
             push!(polys, poly_a_b_hole)
         elseif !b_in_a && !a_in_b # polygons don't intersect
-            push!(polys, tuples(poly_a))
+            push!(polys, svpoints(poly_a, T))
             return polys
         end
     end
@@ -76,7 +76,7 @@ function _difference(
     end
     if GI.nhole(poly_b) != 0
         for hole in GI.gethole(poly_b)
-            hole_poly = GI.Polygon(StaticArrays.SVector(hole))
+            hole_poly = GI.Polygon(SA.SVector(hole))
             new_polys = intersection(hole_poly, poly_a, T; target = GI.PolygonTrait)
             if length(new_polys) > 0
                 append!(polys, new_polys)
@@ -115,10 +115,10 @@ function _difference(
     ::GI.MultiPolygonTrait, multipoly_b;
     kwargs...,
 ) where T
-    polys = [tuples(poly_a, T)]
+    polys = [svpoints(poly_a, T)]
     for poly_b in GI.getpolygon(multipoly_b)
         isempty(polys) && break
-        polys = mapreduce(p -> difference(p, poly_b; target), append!, polys)
+        polys = mapreduce(p -> difference(p, poly_b, T; target), append!, polys)
     end
     return polys
 end
@@ -134,12 +134,12 @@ function _difference(
     fix_multipoly = UnionIntersectingPolygons(), kwargs...,
 ) where T
     if !isnothing(fix_multipoly) # Fix multipoly_a to prevent returning an invalid multipolygon
-        multipoly_a = fix_multipoly(multipoly_a)
+        multipoly_a = fix_multipoly(multipoly_a, T)
     end
-    polys = Vector{_get_poly_type(T)}()
+    polys = Vector{PolyType2D{T}}()
     sizehint!(polys, GI.npolygon(multipoly_a))
     for poly_a in GI.getpolygon(multipoly_a)
-        append!(polys, difference(poly_a, poly_b; target))
+        append!(polys, difference(poly_a, poly_b, T; target))
     end
     return polys
 end
@@ -156,7 +156,7 @@ function _difference(
     fix_multipoly = UnionIntersectingPolygons(), kwargs...,
 ) where T
     if !isnothing(fix_multipoly) # Fix multipoly_a to prevent returning an invalid multipolygon
-        multipoly_a = fix_multipoly(multipoly_a)
+        multipoly_a = fix_multipoly(multipoly_a, T)
         fix_multipoly = nothing
     end
     local polys
diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl
index 5b4444ea9..40f8556dd 100644
--- a/src/methods/clipping/intersection.jl
+++ b/src/methods/clipping/intersection.jl
@@ -75,9 +75,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; exact)
         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))
@@ -118,11 +118,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)}()
+    polys = Vector{PolyType2D{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
@@ -135,7 +135,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
@@ -149,13 +149,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)}()
+    polys = Vector{PolyType2D{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
@@ -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.
 """
@@ -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 = GI.Point[]
+    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
@@ -221,7 +221,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
@@ -242,20 +242,20 @@ 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 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 = ((zero(T), zero(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 = (_tuple_point(a1, T), (zero(T), β1))
-        intr2 = (_tuple_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 = (_tuple_point(b1, T), (α1, zero(T)))
-        intr2 = (_tuple_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 = (_tuple_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 = (_tuple_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 = (_tuple_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 = (_tuple_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 =
     (
-        (_tuple_point(a_pt, T), (a_pt_α, _clamped_frac(distance(a_pt, b1, T), b_dist))),
-        (_tuple_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)
-        _tuple_point(a1, T), zero(T), zero(T)
+        SVPoint_2D(a1, T), zero(T), zero(T)
     elseif equals(a1, b2)
-        _tuple_point(a1, T), zero(T), one(T)
+        SVPoint_2D(a1, T), zero(T), one(T)
     elseif equals(a2, b1)
-        _tuple_point(a2, T), one(T), zero(T)
+        SVPoint_2D(a2, T), one(T), zero(T)
     elseif equals(a2, b2)
-        _tuple_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))
-        _tuple_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))
-        _tuple_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))
-        _tuple_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))
-        _tuple_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 = (x, y)
+    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 _tuple_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/clipping/union.jl b/src/methods/clipping/union.jl
index 8e1d785e4..2563ca6f1 100644
--- a/src/methods/clipping/union.jl
+++ b/src/methods/clipping/union.jl
@@ -62,12 +62,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; exact)
         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
@@ -81,7 +81,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; exact)
+        _add_union_holes!(T, polys, a_in_b, b_in_a, poly_a, poly_b; exact)
     end
     # Remove uneeded collinear points on same edge
     _remove_collinear_points!(polys, [false], poly_a, poly_b)
@@ -108,15 +108,15 @@ _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; exact)
+function _add_union_holes!(::Type{T}, polys, a_in_b, b_in_a, poly_a, poly_b; exact) where T
     if a_in_b
-        _add_union_holes_contained_polys!(polys, poly_a, poly_b; exact)
+        _add_union_holes_contained_polys!(T, polys, poly_a, poly_b; exact)
     elseif b_in_a
-        _add_union_holes_contained_polys!(polys, poly_b, poly_a; exact)
+        _add_union_holes_contained_polys!(T, polys, poly_b, poly_a; exact)
     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 =#
@@ -129,15 +129,15 @@ function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b; exact)
                 #= 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
                 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))
-                new_holes = difference(h_poly, current_poly; target = GI.PolygonTrait())
+                h_poly = GI.Polygon(SA.SVector(ih))
+                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
@@ -152,51 +152,51 @@ 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; exact)
+function _add_union_holes_contained_polys!(::Type{T}, polys, interior_poly, exterior_poly; exact) where T
     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; exact, closed_line = true)
         if in_ih  # at least part of interior polygon exterior is within the ith hole
             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
                 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; exact, 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())
+                            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
@@ -215,21 +215,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
@@ -242,7 +242,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
@@ -255,13 +255,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/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/methods/polygonize.jl b/src/methods/polygonize.jl
index b238f4c63..6c669e46a 100644
--- a/src/methods/polygonize.jl
+++ b/src/methods/polygonize.jl
@@ -109,6 +109,7 @@ function polygonize(f::Base.Callable, xs::AbstractRange, ys::AbstractRange, A::A
     else
         _polygonize_featurecollection(f, xs, ys, A; kw...) 
     end
+    return [GI.Polygon([c]) for c in contours]
 end
 function _polygonize(f::Base.Callable, xs::AbstractRange, ys::AbstractRange, A::AbstractMatrix; 
     kw...
diff --git a/src/primitives.jl b/src/primitives.jl
index 01014c1cf..01de73cf8 100644
--- a/src/primitives.jl
+++ b/src/primitives.jl
@@ -52,7 +52,6 @@ TraitTarget
 ## Implementation
 
 =#
-
 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/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..14b01e63f 100644
--- a/src/transformations/correction/geometry_correction.jl
+++ b/src/transformations/correction/geometry_correction.jl
@@ -40,18 +40,18 @@ 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])
+        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..2668b88c1 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 = 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)(::GI.MultiPolygonTrait, multipoly)
                 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
@@ -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 = svpoints(multipoly, T)
     n_starting_polys = GI.npolygon(multipoly)
     n_polys = n_starting_polys
     if n_polys > 1
@@ -119,7 +119,7 @@ function (::DiffIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly)
                     !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/transformations/flip.jl b/src/transformations/flip.jl
index a8b73766f..2d5a184d3 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))
+            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((GI.y(p), GI.x(p), GI.z(p)))
         end
     else
         return apply(PointTrait(), geom; kw...) do p
-            (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 623637bc5..8da2bc5fd 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, TraitTarget(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 = PointType2D{T}[]
     sizehint!(new_coords, GI.npoint(geom))
-    push!(new_coords, (x1, y1))
+    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!(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, 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, (x2, y2))
+    push!(new_coords, SVPoint_2D((x2, y2), T))
     return nothing
 end
 #=
diff --git a/src/transformations/simplify.jl b/src/transformations/simplify.jl
index 35ccabdb1..0263e90b3 100644
--- a/src/transformations/simplify.jl
+++ b/src/transformations/simplify.jl
@@ -189,45 +189,50 @@ 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...)
 simplify(alg::GEOS, data; kw...) = _simplify(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::Union{SimplifyAlg, GEOS}, 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
+# GEOS only works with Float64 so you cannot pass in a float type argument
+function _simplify(alg::GEOS, data; prefilter_alg=nothing, kw...)
+    simplifier(geom) = _simplify(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,
     )
     lrs = map(simplifier, GI.getgeom(geom))
@@ -445,7 +450,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
 
@@ -503,14 +508,6 @@ function _build_tolerances(f, points)
     return real_tolerances
 end
 
-function tuple_points(geom)
-    points = Array{Tuple{Float64,Float64}}(undef, GI.npoint(geom))
-    for (i, p) in enumerate(GI.getpoint(geom))
-        points[i] = (GI.x(p), GI.y(p))
-    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/transformations/svpoints.jl b/src/transformations/svpoints.jl
new file mode 100644
index 000000000..98269dfae
--- /dev/null
+++ b/src/transformations/svpoints.jl
@@ -0,0 +1,30 @@
+# # SVector Points conversion
+
+"""
+    svpoints(obj)
+
+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
+SVPoints.
+
+# Keywords
+
+$APPLY_KEYWORDS
+"""
+function svpoints(geom, ::Type{T} = Float64; kw...) where T
+    if _ismeasured(geom)
+        return apply(PointTrait(), geom; kw...) do p
+            SVPoint_4D(p)
+        end
+    elseif _is3d(geom)
+        return apply(PointTrait(), geom; kw...) do p
+            SVPoint_3D(p)
+        end
+    else
+        return apply(PointTrait(), geom; kw...) do p
+            SVPoint_2D(p)
+        end
+    end
+end
diff --git a/src/transformations/transform.jl b/src/transformations/transform.jl
index b09e85482..7faa8dc54 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,20 +35,21 @@ 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; 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))))
+            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(SVPoint_3D((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))))
+            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 4bef4f91d..c3c6bf2ed 100644
--- a/src/transformations/tuples.jl
+++ b/src/transformations/tuples.jl
@@ -13,13 +13,17 @@ 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)))
+            TuplePoint_4D(p, T)
+        end
+    elseif _is3d(geom)
+        return apply(PointTrait(), geom; kw...) do 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 4e6a0dec1..e2bc65d60 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -68,6 +68,61 @@ struct _False <: BoolsAsTypes end
 @inline _booltype(x::Bool)::BoolsAsTypes = x ? _True() : _False()
 @inline _booltype(x::BoolsAsTypes)::BoolsAsTypes = x
 
+
+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 = 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`
+
+
+=#
+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 SVEdge{N, T, Z, M} = Tuple{SVPoint{N,T,Z,M}, SVPoint{N,T,Z,M}} where {N,T,Z,M}
+
+# 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."))
+
+# 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_2D(vals, ::Type{T} = Float64) where T <: AbstractFloat = PointType2D{T}(TuplePoint_2D(vals, T))
+
+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))
+
+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
+=#
+const PolyType2D{T} = GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{PointType2D{T}}, Nothing, Nothing}}, Nothing, Nothing} 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 82a1e664d..b6534b818 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)
@@ -49,6 +54,10 @@ function polygon_to_line(poly)
     return GI.LineString(collect(GI.getgeom(GI.getgeom(poly, 1))))
 end
 
+#= 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`=#
 
 """
     to_edges()
@@ -56,12 +65,23 @@ end
 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))
+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(x)
+#         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)
@@ -76,10 +96,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
@@ -87,10 +107,7 @@ 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})
+function to_extent(edges::Vector{<:Edge})
     x, y = extrema(first, edges)
     Extents.Extent(X=x, Y=y)
 end
@@ -117,11 +134,19 @@ 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{PointType2D{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
diff --git a/test/methods/centroid.jl b/test/methods/centroid.jl
index 431f4732b..d7f7c453f 100644
--- a/test/methods/centroid.jl
+++ b/test/methods/centroid.jl
@@ -9,16 +9,16 @@ import GeoInterface as GI,
     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
@@ -28,15 +28,14 @@ import GeoInterface as GI,
     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
 
@@ -47,8 +46,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
@@ -61,8 +60,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
@@ -72,8 +71,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
@@ -84,8 +83,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([
@@ -94,7 +93,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
@@ -110,7 +110,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/flip.jl b/test/transformations/flip.jl
index e977c7504..16259d2b6 100644
--- a/test/transformations/flip.jl
+++ b/test/transformations/flip.jl
@@ -7,6 +7,6 @@ import GeoInterface as GI, 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/segmentize.jl b/test/transformations/segmentize.jl
index 3a5d18333..6f8713c42 100644
--- a/test/transformations/segmentize.jl
+++ b/test/transformations/segmentize.jl
@@ -54,7 +54,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
diff --git a/test/transformations/transform.jl b/test/transformations/transform.jl
index f3d5b50ef..f75ca1a3b 100644
--- a/test/transformations/transform.jl
+++ b/test/transformations/transform.jl
@@ -9,5 +9,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