diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index d626e7bc2..b1e78fb70 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -28,6 +28,10 @@ include("methods/bools.jl")
 include("methods/centroid.jl")
 include("methods/distance.jl")
 include("methods/equals.jl")
+include("methods/clipping/clipping_processor.jl")
+include("methods/clipping/intersection.jl")
+include("methods/clipping/difference.jl")
+include("methods/clipping/union.jl")
 include("methods/geom_relations/contains.jl")
 include("methods/geom_relations/coveredby.jl")
 include("methods/geom_relations/covers.jl")
diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
new file mode 100644
index 000000000..3e787880e
--- /dev/null
+++ b/src/methods/clipping/clipping_processor.jl
@@ -0,0 +1,291 @@
+# # This file contains the shared helper functions forlyNode the polygon clipping functionalities.
+
+#= This is the struct that makes up a_list and b_list. Many values are only used if point is
+an intersection point (ipt). =#
+struct PolyNode{T <: AbstractFloat}
+    idx::Int           # If ipt, index of point in a_idx_list, else 0
+    point::Tuple{T,T}  # (x, y) values of given point
+    inter::Bool        # If ipt, true, else 0
+    neighbor::Int      # If ipt, index of equivalent point in a_list or b_list, else 0
+    ent_exit::Bool     # If ipt, true if enter and false if exit, else false
+    fracs::Tuple{T,T}  # If ipt, fractions along edges to ipt (a_frac, b_frac), else (0, 0)
+end
+
+#=
+    _build_ab_list(::Type{T}, poly_a, poly_b) -> (a_list, b_list, a_idx_list)
+
+This function takes in two polygon rings and calls '_build_a_list', '_build_b_list', and
+'_flag_ent_exit' in order to fully form a_list and b_list. The 'a_list' and 'b_list' that it
+returns are the fully updated vectors of PolyNodes that represent the rings 'poly_a' and
+'poly_b', respectively. This function also returns 'a_idx_list', which at its "ith" index
+stores the index in 'a_list' at which the "ith" intersection point lies.
+=#
+function _build_ab_list(::Type{T}, poly_a, poly_b) where T
+    # Make a list for nodes of each polygon
+    a_list, a_idx_list = _build_a_list(T, poly_a, poly_b)
+    b_list = _build_b_list(T, a_idx_list, a_list, poly_b)
+
+    # Flag the entry and exits
+    _flag_ent_exit!(poly_b, a_list)
+    _flag_ent_exit!(poly_a, b_list)
+
+    return a_list, b_list, a_idx_list
+end
+
+#=
+    _build_a_list(::Type{T}, poly_a, poly_b) -> (a_list, a_idx_list)
+
+This function take in two polygon rings and creates a vector of PolyNodes to represent
+poly_a, including its intersection points with poly_b. The information stored in each
+PolyNode is needed for clipping using the Greiner-Hormann clipping algorithm.
+    
+Note: After calling this function, a_list is not fully formed because the neighboring
+indicies of the intersection points in b_list still need to be updated. Also we still have
+not update the entry and exit flags for a_list.
+    
+The a_idx_list is a list of the indicies of intersection points in a_list. The value at
+index i of a_idx_list is the location in a_list where the ith intersection point lies.
+=#
+function _build_a_list(::Type{T}, poly_a, poly_b) where T
+    n_a_edges = _nedge(poly_a)
+    a_list = Vector{PolyNode{T}}(undef, n_a_edges)  # list of points in poly_a
+    a_idx_list = Vector{Int}()  # finds indices of intersection points in a_list
+    intr_count = 0  # number of intersection points found
+    a_count = 0  # number of points added to a_list
+    # 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)))
+        if i <= 1
+            a_pt1 = a_pt2
+            continue
+        end
+        # Add the first point of the edge to the list of points in a_list
+        new_point = PolyNode(0, a_pt1, false, 0, false, (zero(T), zero(T)))
+        a_count += 1
+        _add!(a_list, a_count, new_point, n_a_edges)
+        # Find intersections with edges of poly_b
+        local b_pt1
+        prev_counter = intr_count
+        for (j, b_p2) in enumerate(GI.getpoint(poly_b))
+            b_pt2 = _tuple_point(b_p2)
+            if j <=1
+                b_pt1 = b_pt2
+                continue
+            end
+            int_pt, fracs = _intersection_point(T, (a_pt1, a_pt2), (b_pt1, b_pt2))
+            # if no intersection point, skip this edge
+            if !isnothing(int_pt) && all(0 .≤ fracs .≤ 1)
+                # Set neighbor field to b edge (j-1) to keep track of intersection
+                new_intr = PolyNode(intr_count, int_pt, true, j - 1, false, fracs)
+                a_count += 1
+                intr_count += 1
+                _add!(a_list, a_count, new_intr, n_a_edges)
+                push!(a_idx_list, a_count)
+            end
+            b_pt1 = b_pt2
+        end
+
+        # Order intersection points by placement along edge using fracs value
+        if prev_counter < intr_count
+            Δintrs = intr_count - prev_counter
+            inter_points = @view a_list[(a_count - Δintrs + 1):a_count]
+            sort!(inter_points, by = x -> x.fracs[1])
+            for (i, p) in enumerate(inter_points)
+                inter_points[i] = PolyNode(prev_counter + i, p.point, p.inter, p.neighbor, p.ent_exit, p.fracs)
+            end
+        end
+    
+        a_pt1 = a_pt2
+    end
+    return a_list, a_idx_list
+end
+
+# Add value x at index i to given array - if list isn't long enough, push value to array
+function _add!(arr, i, x, l = length(arr))
+    if i <= l
+        arr[i] = x
+    else
+        push!(arr, x)
+    end
+    return
+end
+
+#=
+    _build_b_list(::Type{T}, a_idx_list, a_list, poly_b) -> b_list
+
+This function takes in the a_list and a_idx_list build in _build_a_list and poly_b and
+creates a vector of PolyNodes to represent poly_b. The information stored in each PolyNode
+is needed for clipping using the Greiner-Hormann clipping algorithm.
+    
+Note: after calling this function, b_list is not fully updated. The entry/exit flags still
+need to be updated. However, the neightbor value in a_list is now updated.
+=#
+function _build_b_list(::Type{T}, a_idx_list, a_list, poly_b) where T
+    # Sort intersection points by insertion order in b_list
+    sort!(a_idx_list, by = x-> a_list[x].neighbor + a_list[x].fracs[2])
+    # Initialize needed values and lists
+    n_b_edges = _nedge(poly_b)
+    n_intr_pts = length(a_idx_list)
+    b_list = Vector{PolyNode{T}}(undef, n_b_edges + n_intr_pts)
+    intr_curr = 1
+    b_count = 0
+    # Loop over points in poly_b and add each point and intersection point
+    for (i, p) in enumerate(GI.getpoint(poly_b))
+        (i == n_b_edges + 1) && break
+        b_count += 1
+        pt = (T(GI.x(p)), T(GI.y(p)))
+        b_list[b_count] = PolyNode(0, pt, false, 0, false, (zero(T), zero(T)))
+        if intr_curr ≤ n_intr_pts
+            curr_idx = a_idx_list[intr_curr]
+            curr_node = a_list[curr_idx]
+            while curr_node.neighbor == i  # Add all intersection points in current edge
+                b_count += 1
+                b_list[b_count] = PolyNode(curr_node.idx, curr_node.point, true, curr_idx, false, curr_node.fracs)
+                a_list[curr_idx] = PolyNode(curr_node.idx, curr_node.point, curr_node.inter, b_count, curr_node.ent_exit, curr_node.fracs)
+                curr_node = a_list[curr_idx]
+                intr_curr += 1
+                intr_curr > n_intr_pts && break
+                curr_idx = a_idx_list[intr_curr]
+                curr_node = a_list[curr_idx]
+            end
+        end
+    end
+    sort!(a_idx_list)  # return a_idx_list to order of points in a_list
+    return b_list
+end
+
+#=
+    _flag_ent_exit(poly_b, a_list)
+
+This function flags all the intersection points as either an 'entry' or 'exit' point in
+relation to the given polygon.
+=#
+function _flag_ent_exit!(poly, pt_list)
+    local status
+    for ii in eachindex(pt_list)
+        if ii == 1
+            status = !_point_filled_curve_orientation(
+                pt_list[ii].point, poly;
+                in = true, on = false, out = false
+            )
+        elseif pt_list[ii].inter
+            pt_list[ii] = PolyNode(pt_list[ii].idx, pt_list[ii].point, pt_list[ii].inter, pt_list[ii].neighbor, status, pt_list[ii].fracs)
+            status = !status
+        end
+    end
+    return
+end
+
+#=
+    _trace_polynodes(::Type{T}, a_list, b_list, a_idx_list, f_step)::Vector{GI.Polygon}
+
+This function takes the outputs of _build_ab_list and traces the lists to determine which
+polygons are formed as described in Greiner and Hormann. The function f_step determines in
+which direction the lists are traced.  This function is different for intersection,
+difference, and union. f_step must take in two arguments: the most recent intersection
+node's entry/exit status and a boolean that is true if we are currently tracing a_list and
+false if we are tracing b_list. The functions used for each clipping operation are follows:
+    - Intersection: (x, y) -> x ? 1 : (-1)
+    - Difference: (x, y) -> (x ⊻ y) ? 1 : (-1)
+    - Union: (x, y) -> (x ⊻ y) ? 1 : (-1)
+
+A list of GeoInterface polygons is returned from this function. 
+=#
+function _trace_polynodes(::Type{T}, a_list, b_list, a_idx_list, f_step) where T
+    n_a_pts, n_b_pts = length(a_list), length(b_list)
+    n_intr_pts = length(a_idx_list)
+    return_polys = Vector{_get_poly_type(T)}(undef, 0)
+    # Keep track of number of processed intersection points
+    processed_pts = 0
+    while processed_pts < n_intr_pts
+        curr_list, curr_npoints = a_list, n_a_pts
+        on_a_list = true
+        # Find first unprocessed intersecting point in subject polygon
+        processed_pts += 1
+        first_idx = findnext(x -> x != 0, a_idx_list, processed_pts)
+        idx = a_idx_list[first_idx]
+        a_idx_list[first_idx] = 0
+        start_pt = a_list[idx]
+
+        # Set first point in polygon
+        curr = curr_list[idx]
+        pt_list = [curr.point]
+
+        curr_not_start = true
+        while curr_not_start
+            step = f_step(curr.ent_exit, on_a_list)
+            curr_not_intr = true
+            while curr_not_intr
+                # Traverse polygon either forwards or backwards
+                idx += step
+                idx = (idx > curr_npoints) ? mod(idx, curr_npoints) : idx
+                idx = (idx == 0) ? curr_npoints : idx
+
+                # Get current node and add to pt_list
+                curr = curr_list[idx]
+                push!(pt_list, curr.point)
+                if curr.inter 
+                    # Keep track of processed intersection points
+                    curr_not_start = curr != start_pt && curr != b_list[start_pt.neighbor]
+                    if curr_not_start
+                        processed_pts += 1
+                        a_idx_list[curr.idx] = 0
+                    end
+                    curr_not_intr = false
+                end
+            end
+
+            # Switch to next list and next point
+            curr_list, curr_npoints = on_a_list ? (b_list, n_b_pts) : (a_list, n_a_pts)
+            on_a_list = !on_a_list
+            idx = curr.neighbor
+            curr = curr_list[idx]
+        end
+        push!(return_polys, GI.Polygon([pt_list]))
+    end
+    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}
+
+#=
+    _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator)
+
+The holes specified by the hole iterator are added to the polygons in the return_polys list.
+If this creates more polygon, they are added to the end of the list. If this removes
+polygons, they are removed from the list
+=#
+function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator) where T
+    n_polys = length(return_polys)
+    # Remove set of holes from all polygons
+    for i in 1:n_polys
+        n_new_per_poly = 0
+        for hole in hole_iterator  # loop through all holes
+            hole_poly = GI.Polygon([hole])
+            # 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)))
+                if !isnothing(return_polys[j])
+                    new_polys = difference(return_polys[j], hole_poly, T; target = GI.PolygonTrait)
+                    n_new_polys = length(new_polys)
+                    if n_new_polys == 0  # hole covered whole polygon
+                        return_polys[j] = nothing
+                    else
+                        return_polys[j] = new_polys[1]  # replace original
+                        if n_new_polys > 1  # add any extra pieces
+                            append!(return_polys, @view new_polys[2:end])
+                            n_new_per_poly += n_new_polys - 1
+                        end
+                    end
+                end
+            end
+        end
+        n_polys += n_new_per_poly
+    end
+    # Remove all polygon that were marked for removal
+    filter!(!isnothing, return_polys)
+    return
+end
\ No newline at end of file
diff --git a/src/methods/clipping/difference.jl b/src/methods/clipping/difference.jl
new file mode 100644
index 000000000..4e290d7d3
--- /dev/null
+++ b/src/methods/clipping/difference.jl
@@ -0,0 +1,83 @@
+# #  Difference Polygon Clipping
+export difference
+
+
+"""
+    difference(geom_a, geom_b, [T::Type]; target::Type)
+
+Return the difference between two geometries as a list of geometries. Return an empty list
+if none are found. The type of the list will be constrained as much as possible given the
+input geometries. Furthermore, the user can provide a `taget` type as a keyword argument and
+a list of target geometries found in the difference will be returned. The user can also
+provide a float type that they would like the points of returned geometries to be. 
+
+## Example 
+
+```jldoctest
+import GeoInterface as GI, GeometryOps as GO
+
+poly1 = GI.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]])
+poly2 = GI.Polygon([[[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]])
+diff_poly = GO.difference(poly1, poly2; target = GI.PolygonTrait)
+GI.coordinates.(diff_poly)
+
+# output
+1-element Vector{Vector{Vector{Vector{Float64}}}}:
+ [[[6.5, 3.5], [5.0, 5.0], [0.0, 0.0], [5.0, -5.0], [6.5, -3.5], [3.0, 0.0], [6.5, 3.5]]]
+```
+"""
+function difference(
+    geom_a, geom_b, ::Type{T} = Float64; target::Type{Target} = Nothing,
+) where {T <: AbstractFloat, Target <: GI.AbstractTrait}
+    return _difference(Target, T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
+end
+
+#= The 'difference' function returns the difference of two polygons as a list of polygons.
+The algorithm to determine the difference was adapted from "Efficient clipping of efficient
+polygons," by Greiner and Hormann (1998). DOI: https://doi.org/10.1145/274363.274364 =#
+function _difference(
+    ::Type{GI.PolygonTrait}, ::Type{T},
+    ::GI.PolygonTrait, poly_a,
+    ::GI.PolygonTrait, poly_b,
+) where T
+    # Get the exterior of the polygons
+    ext_poly_a = GI.getexterior(poly_a)
+    ext_poly_b = GI.getexterior(poly_b)
+    # Find the difference of the exterior of the polygons
+    a_list, b_list, a_idx_list = _build_ab_list(T, ext_poly_a, ext_poly_b)
+    polys = _trace_polynodes(T, a_list, b_list, a_idx_list, (x, y) -> (x ⊻ y) ? 1 : (-1))
+    if isempty(polys)
+        if _point_filled_curve_orientation(b_list[1].point, ext_poly_a) == point_in
+            poly_a_b_hole = GI.Polygon([ext_poly_a, ext_poly_b])
+            push!(polys, poly_a_b_hole)
+        elseif _point_filled_curve_orientation(a_list[1].point, ext_poly_b) != point_in
+            # Two polygons don't intersect and are not contained in one another
+            push!(polys, GI.Polygon([ext_poly_a]))
+        end
+    end
+
+    # If the original polygons had holes, take that into account.
+    if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0
+        _add_holes_to_polys!(T, polys, GI.gethole(poly_a))
+        for hole in GI.gethole(poly_b)
+            new_polys = intersection(GI.Polygon([hole]), poly_a, T; target = GI.PolygonTrait)
+            if length(new_polys) > 0
+                append!(polys, new_polys)
+            end
+        end
+    end
+    return polys
+end
+
+# Many type and target combos aren't implemented
+function _difference(
+    ::Type{Target}, ::Type{T},
+    trait_a::GI.AbstractTrait, geom_a,
+    trait_b::GI.AbstractTrait, geom_b,
+) where {Target, T}
+    @assert(
+        false,
+        "Intersection between $trait_a and $trait_b with target $Target isn't implemented yet.",
+    )
+    return nothing
+end
\ No newline at end of file
diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl
new file mode 100644
index 000000000..8f81a699a
--- /dev/null
+++ b/src/methods/clipping/intersection.jl
@@ -0,0 +1,174 @@
+# #  Intersection
+export intersection, intersection_points
+
+
+"""
+    intersection(geom_a, geom_b, [T::Type]; target::Type)
+
+Return the intersection between two geometries as a list of geometries. Return an empty list
+if none are found. The type of the list will be constrained as much as possible given the
+input geometries. Furthermore, the user can provide a `taget` type as a keyword argument and
+a list of target geometries found in the intersection will be returned. The user can also
+provide a float type that they would like the points of returned geometries to be. 
+
+## Example
+
+```jldoctest
+import GeoInterface as GI, GeometryOps as GO
+
+line1 = GI.Line([(124.584961,-12.768946), (126.738281,-17.224758)])
+line2 = GI.Line([(123.354492,-15.961329), (127.22168,-14.008696)])
+inter_points = GO.intersection(line1, line2; target = GI.PointTrait)
+GI.coordinates.(inter_points)
+
+# output
+1-element Vector{Vector{Float64}}:
+ [125.58375366067547, -14.83572303404496]
+```
+"""
+function intersection(
+    geom_a, geom_b, ::Type{T} = Float64; target::Type{Target} = Nothing,
+) where {T <: AbstractFloat, Target <: GI.AbstractTrait}
+    return _intersection(Target, T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
+end
+
+# Curve-Curve Intersections with target Point
+_intersection(
+    ::Type{GI.PointTrait}, ::Type{T},
+    trait_a::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, geom_a,
+    trait_b::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, geom_b,
+) where T = _intersection_points(T, trait_a, geom_a, trait_b, geom_b)
+
+
+#= Polygon-Polygon Intersections with target Polygon
+The algorithm to determine the intersection was adapted from "Efficient clipping
+of efficient polygons," by Greiner and Hormann (1998).
+DOI: https://doi.org/10.1145/274363.274364 =#
+function _intersection(
+    ::Type{GI.PolygonTrait}, ::Type{T},
+    ::GI.PolygonTrait, poly_a,
+    ::GI.PolygonTrait, poly_b,
+) where {T}
+    # First we get the exteriors of 'poly_a' and 'poly_b'
+    ext_poly_a = GI.getexterior(poly_a)
+    ext_poly_b = GI.getexterior(poly_b)
+    # Then we find the intersection of the exteriors
+    a_list, b_list, a_idx_list = _build_ab_list(T, ext_poly_a, ext_poly_b)
+    polys = _trace_polynodes(T, a_list, b_list, a_idx_list, (x, y) -> x ? 1 : (-1))
+
+    if isempty(polys)
+        if _point_filled_curve_orientation(a_list[1].point, ext_poly_b) == point_in
+            push!(polys, GI.Polygon([ext_poly_a]))
+        elseif _point_filled_curve_orientation(b_list[1].point, ext_poly_a) == point_in
+            push!(polys, GI.Polygon([ext_poly_b]))
+        end
+    end
+    # If the original polygons had holes, take that into account.
+    if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0
+        hole_iterator = Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b)))
+        _add_holes_to_polys!(T, polys, hole_iterator)
+    end    
+    return polys
+end
+
+# Many type and target combos aren't implemented
+function _intersection(
+    ::Type{Target}, ::Type{T},
+    trait_a::GI.AbstractTrait, geom_a,
+    trait_b::GI.AbstractTrait, geom_b,
+) where {Target, T}
+    @assert(
+        false,
+        "Intersection between $trait_a and $trait_b with target $Target isn't implemented yet.",
+    )
+    return nothing
+end
+
+"""
+    intersection_points(
+        geom_a,
+        geom_b,
+    )::Union{
+        ::Vector{::Tuple{::Real, ::Real}},
+        ::Nothing,
+    }
+
+Return a list of intersection points between two geometries of type GI.Point.
+If no intersection point was possible given geometry extents, returns an empty
+list.
+"""
+intersection_points(geom_a, geom_b, ::Type{T} = Float64) where T <: AbstractFloat =
+    _intersection_points(T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
+
+
+#= Calculates the list of intersection points between two geometries, inlcuding line
+segments, line strings, linear rings, polygons, and multipolygons. If no intersection points
+were possible given geometry extents or if none are found, return an empty list of
+GI.Points. =#
+function _intersection_points(::Type{T}, ::GI.AbstractTrait, a, ::GI.AbstractTrait, b) where T
+    # Initialize an empty list of points
+    result = GI.Point[]
+    # 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
+    edges_a, edges_b = map(sort! ∘ to_edges, (a, b))
+    npoints_a, npoints_b  = length(edges_a), length(edges_b)
+    a_closed = npoints_a > 1 && edges_a[1][1] == edges_a[end][1]
+    b_closed = npoints_b > 1 && edges_b[1][1] == edges_b[end][1]
+    if npoints_a > 0 && npoints_b > 0
+        # Loop over pairs of edges and add any intersection points to results
+        for i in eachindex(edges_a), j in eachindex(edges_b)
+            point, fracs = _intersection_point(T, edges_a[i], edges_b[j])
+            if !isnothing(point)
+                #=
+                Determine if point is on edge (all edge endpoints excluded
+                except for the last edge for an open geometry)
+                =#
+                α, β = fracs
+                on_a_edge = (!a_closed && i == npoints_a && 0 <= α <= 1) ||
+                    (0 <= α < 1)
+                on_b_edge = (!b_closed && j == npoints_b && 0 <= β <= 1) ||
+                    (0 <= β < 1)
+                if on_a_edge && on_b_edge
+                    push!(result, GI.Point(point))
+                end
+            end
+        end
+    end
+    return result
+end
+
+#= Calculates the intersection point between two lines if it exists, and as if the line
+extended to infinity, and the fractional component of each line from the initial end point
+to the intersection point.
+Inputs:
+    (a1, a2)::Tuple{Tuple{::Real, ::Real}, Tuple{::Real, ::Real}} first line
+    (b1, b2)::Tuple{Tuple{::Real, ::Real}, Tuple{::Real, ::Real}} second line
+Outputs:
+    (x, y)::Tuple{::Real, ::Real} intersection point
+    (t, u)::Tuple{::Real, ::Real} fractional length of lines to intersection
+    Both are ::Nothing if point doesn't exist!
+
+Calculation derivation can be found here:
+    https://stackoverflow.com/questions/563198/
+=#
+function _intersection_point(::Type{T}, (a1, a2)::Tuple, (b1, b2)::Tuple) where T
+    # First line runs from p to p + r
+    px, py = GI.x(a1), GI.y(a1)
+    rx, ry = GI.x(a2) - px, GI.y(a2) - py
+    # Second line runs from q to q + s 
+    qx, qy = GI.x(b1), GI.y(b1)
+    sx, sy = GI.x(b2) - qx, GI.y(b2) - qy
+    # Intersection will be where p + tr = q + us where 0 < t, u < 1 and
+    r_cross_s = rx * sy - ry * sx
+    if r_cross_s != 0
+        Δqp_x = qx - px
+        Δqp_y = qy - py
+        t = (Δqp_x * sy - Δqp_y * sx) / r_cross_s
+        u = (Δqp_x * ry - Δqp_y * rx) / r_cross_s
+        x = px + t * rx
+        y = py + t * ry
+        return (T(x), T(y)), (T(t), T(u))
+    end
+    return nothing, nothing
+end
\ No newline at end of file
diff --git a/src/methods/clipping/union.jl b/src/methods/clipping/union.jl
new file mode 100644
index 000000000..c76b84ddf
--- /dev/null
+++ b/src/methods/clipping/union.jl
@@ -0,0 +1,92 @@
+# #  Union Polygon Clipping
+export union
+
+"""
+    union(geom_a, geom_b, [::Type{T}]; target::Type)
+
+Return the union between two geometries as a list of geometries. Return an empty list if
+none are found. The type of the list will be constrained as much as possible given the input
+geometries. Furthermore, the user can provide a `taget` type as a keyword argument and a
+list of target geometries found in the difference will be returned. The user can also
+provide a float type 'T' that they would like the points of returned geometries to be. 
+    
+Calculates the union between two polygons.
+## Example
+
+```jldoctest
+import GeoInterface as GI, GeometryOps as GO
+
+p1 = GI.Polygon([[(0.0, 0.0), (5.0, 5.0), (10.0, 0.0), (5.0, -5.0), (0.0, 0.0)]])
+p2 = GI.Polygon([[(3.0, 0.0), (8.0, 5.0), (13.0, 0.0), (8.0, -5.0), (3.0, 0.0)]])
+union_poly = GO.union(p1, p2; target = GI.PolygonTrait)
+GI.coordinates.(union_poly)
+
+# output
+1-element Vector{Vector{Vector{Vector{Float64}}}}:
+ [[[6.5, 3.5], [5.0, 5.0], [0.0, 0.0], [5.0, -5.0], [6.5, -3.5], [8.0, -5.0], [13.0, 0.0], [8.0, 5.0], [6.5, 3.5]]]
+```
+"""
+function union(
+    geom_a, geom_b, ::Type{T} = Float64; target::Type{Target} = Nothing,
+) where {T <: AbstractFloat, Target <: GI.AbstractTrait}
+    _union(Target, T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
+end
+
+#= This 'union' implementation returns the union of two polygons. The algorithm to determine
+the union was adapted from "Efficient clipping of efficient polygons," by Greiner and
+Hormann (1998). DOI: https://doi.org/10.1145/274363.274364 =#
+function _union(
+    ::Type{GI.PolygonTrait}, ::Type{T},
+    ::GI.PolygonTrait, poly_a,
+    ::GI.PolygonTrait, poly_b,
+) where T
+    # First, I get the exteriors of the two polygons
+    ext_poly_a = GI.getexterior(poly_a)
+    ext_poly_b = GI.getexterior(poly_b)
+    # Then, I get the union of the exteriors
+    a_list, b_list, a_idx_list = _build_ab_list(T, ext_poly_a, ext_poly_b)
+    polys = _trace_polynodes(T, a_list, b_list, a_idx_list, (x, y) -> x ? (-1) : 1)
+    # Check if one polygon totally within other and if so, return the larger polygon.
+    if isempty(polys)
+        if _point_filled_curve_orientation(a_list[1].point, ext_poly_b) == point_in
+            push!(polys, GI.Polygon([ext_poly_b]))
+        elseif _point_filled_curve_orientation(b_list[1].point, ext_poly_a) == point_in
+            push!(polys,  GI.Polygon([ext_poly_a]))
+        else
+            push!(polys, poly_a)
+            push!(polys, poly_b)
+            return polys
+        end
+    else
+        sort!(polys, by = area, rev = true)
+        polys = [GI.Polygon([GI.getexterior(p) for p in polys])]
+    end
+
+    n_b_holes = GI.nhole(poly_b)
+    if GI.nhole(poly_a) != 0 || n_b_holes != 0
+        new_poly = [GI.getexterior(polys[1]); collect(GI.gethole(polys[1]))]
+        current_poly = GI.Polygon([ext_poly_b])
+        for (i, hole) in enumerate(Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b))))
+            # Use ext_poly_b to not overcount overlapping holes in poly_a and in poly_b
+            new_hole = difference(GI.Polygon([hole]), current_poly, T; target = GI.PolygonTrait)
+            for h in new_hole
+                push!(new_poly, GI.getexterior(h))
+            end
+            if i == n_b_holes
+                current_poly = poly_a
+            end
+        end
+        polys[1] = GI.Polygon(new_poly)
+    end
+    return polys
+end
+
+# Many type and target combos aren't implemented
+function _union(
+    ::Type{Target}, ::Type{T},
+    trait_a::GI.AbstractTrait, geom_a,
+    trait_b::GI.AbstractTrait, geom_b,
+) where {Target, T}
+    throw(ArgumentError("Intersection between $trait_a and $trait_b with target $Target isn't implemented yet."))
+    return nothing
+end
\ No newline at end of file
diff --git a/src/methods/geom_relations/crosses.jl b/src/methods/geom_relations/crosses.jl
index b9f0cbcab..0e122337a 100644
--- a/src/methods/geom_relations/crosses.jl
+++ b/src/methods/geom_relations/crosses.jl
@@ -43,7 +43,6 @@ function multipoint_crosses_line(geom1, geom2)
         end
         i += 1
     end
-
     return int_point && ext_point
 end
 
diff --git a/src/methods/geom_relations/geom_geom_processors.jl b/src/methods/geom_relations/geom_geom_processors.jl
index d26818083..bc631be71 100644
--- a/src/methods/geom_relations/geom_geom_processors.jl
+++ b/src/methods/geom_relations/geom_geom_processors.jl
@@ -123,14 +123,14 @@ function _line_curve_process(
     closed_line |= first_last_equal_line
     closed_curve |= first_last_equal_curve
     # Loop over each line segment
-    l_start = GI.getpoint(line, closed_line ? nl : 1)
+    l_start = _tuple_point(GI.getpoint(line, closed_line ? nl : 1))
     i = closed_line ? 1 : 2
     while i ≤ nl
-        l_end = GI.getpoint(line, i)
-        c_start = GI.getpoint(curve, closed_curve ? nc : 1)
+        l_end = _tuple_point(GI.getpoint(line, i))
+        c_start = _tuple_point(GI.getpoint(curve, closed_curve ? nc : 1))
         # Loop over each curve segment
         for j in (closed_curve ? 1 : 2):nc
-            c_end = GI.getpoint(curve, j)
+            c_end = _tuple_point(GI.getpoint(curve, j))
             # Check if line and curve segments meet
             seg_val = _segment_segment_orientation((l_start, l_end), (c_start, c_end))
             # If segments are co-linear
@@ -208,11 +208,12 @@ function _find_new_seg(i, ls, le, cs, ce)
 end
 
 #= Find where line and curve segments intersect by fraction of length. α is the fraction of
-the line (ls to le) and β is the traction of the curve (cs to ce). =#
+the line (ls to le) and β is the traction of the curve (cs to ce). All inputs are tuples. =#
 function _find_intersect_fracs(ls, le, cs, ce)
     _, fracs = _intersection_point(
-        (_tuple_point(ls), _tuple_point(le)),
-        (_tuple_point(cs), _tuple_point(ce))
+        Float64,
+        (ls, le),
+        (cs, ce)
     )
     (α, β) = if !isnothing(fracs)
         fracs
@@ -231,17 +232,18 @@ function _find_intersect_fracs(ls, le, cs, ce)
 end
 
 #= Find next set of segments needed to determine if given hinge segments cross or not.=#
-_find_hinge_next_segments(α, β, ls, le, cs, ce, i, line, j, curve) = 
-    if β == 1
+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, GI.getpoint(line, i + 1)), (ce, GI.getpoint(curve, j + 1)))
+            ((le, _tuple_point(GI.getpoint(line, i + 1))), (ce, _tuple_point(GI.getpoint(curve, j + 1))))
         else  # hinge at curve endpoint and line interior point, curve next segment needed 
-            ((ls, le), (ce, GI.getpoint(curve, j + 1)))
+            ((ls, le), (ce, _tuple_point(GI.getpoint(curve, j + 1))))
         end
     else  # hinge at curve interior point and line endpoint, line next segment needed
-        ((le, GI.getpoint(line, i + 1)), (cs, ce))
+        ((le, _tuple_point(GI.getpoint(line, i + 1))), (cs, ce))
     end
-
+    return next_seg
+end
 #=
 Determines if a line meets the given checks with respect to a polygon.
 
@@ -540,10 +542,10 @@ function _segment_segment_orientation(
     cross::T = line_cross, hinge::T = line_hinge,
     over::T = line_over, out::T = line_out,
 ) where T
-    (ax, ay) = Float64.(_tuple_point(a_point))
-    (bx, by) = Float64.(_tuple_point(b_point))
-    (cx, cy) = Float64.(_tuple_point(c_point))
-    (dx, dy) = Float64.(_tuple_point(d_point))
+    (ax, ay) = Float64.(a_point)
+    (bx, by) = Float64.(b_point)
+    (cx, cy) = Float64.(c_point)
+    (dx, dy) = Float64.(d_point)
     meet_type = ExactPredicates.meet((ax, ay), (bx, by), (cx, cy), (dx, dy))
     # Lines meet at one point within open segments 
     meet_type == 1 && return cross
@@ -600,7 +602,7 @@ function _line_filled_curve_interactions(
     closed_line |= first_last_equal_line
 
     # See if first point is in an acceptable orientation
-    l_start = GI.getpoint(line, closed_line ? nl : 1)
+    l_start = _tuple_point(GI.getpoint(line, closed_line ? nl : 1))
     point_val = _point_filled_curve_orientation(l_start, curve)
     if point_val == point_in
         in_curve = true
@@ -612,13 +614,13 @@ function _line_filled_curve_interactions(
 
     # Check for any intersections between line and curve
     for i in (closed_line ? 1 : 2):nl
-        l_end = GI.getpoint(line, i)
-        c_start = GI.getpoint(curve, nc)
+        l_end = _tuple_point(GI.getpoint(line, i))
+        c_start = _tuple_point(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 = GI.getpoint(curve, j)
+            c_end = _tuple_point(GI.getpoint(curve, j))
             # Check if two line and curve segments meet
             seg_val = _segment_segment_orientation(
                 (l_start, l_end),
@@ -659,8 +661,8 @@ function _line_filled_curve_interactions(
                         p_start = _tuple_point(l_start)
                         for i in 1:(npoints + 1)
                             p_end = i ≤ npoints ?
-                                ipoints[i] :
-                                _tuple_point(l_end)
+                                _tuple_point(ipoints[i]) :
+                                l_end
                             mid_val = _point_filled_curve_orientation(
                                 (p_start .+ p_end) ./ 2,
                                 curve,
diff --git a/src/methods/geom_relations/intersects.jl b/src/methods/geom_relations/intersects.jl
index 8e0a71027..83bf4a98c 100644
--- a/src/methods/geom_relations/intersects.jl
+++ b/src/methods/geom_relations/intersects.jl
@@ -1,25 +1,14 @@
 # # Intersection checks
 
-export intersects, intersection, intersection_points
+export intersects
 
 #=
-## What is `intersects` vs `intersection` vs `intersection_points`?
+## What is `intersects`?
 
 The intersects function checks if a given geometry intersects with another
 geometry, or in other words, the either the interiors or boundaries of the two
 geometries intersect.
 
-The `intersection` methods return the geometry intersection between the two
-input geometries. The `intersection_points` method returns a list of
-intersection points between two geometries.
-
-The `intersects` methods will always return a Boolean. However, note that the
-`intersection` methods will not all return the same type. For example, the
-intersection of two lines will be a point in most cases, unless the lines are
-parallel. On the other hand, the intersection of two polygons will be another
-polygon in most cases. Finally, the `intersection_points` method returns a list
-of tuple points.
-
 To provide an example, consider these two lines:
 ```@example intersects_intersection
 import GeometryOps as GO
@@ -46,7 +35,6 @@ This is the GeoInterface-compatible implementation.
 Given that intersects is the exact opposite of disjoint, we simply pass the two
 inputs variables, swapped in order, to disjoint.
 =#
-
 """
     intersects(geom1, geom2)::Bool
 
@@ -70,207 +58,23 @@ true
 intersects(geom1, geom2) = !disjoint(geom1, geom2)
 
 
-"""
-    intersection(geom_a, geom_b)::Union{Tuple{::Real, ::Real}, ::Nothing}
-
-Return an intersection point between two geometries. Return nothing if none are
-found. Else, the return type depends on the input. It will be a union between:
-a point, a line, a linear ring, a polygon, or a multipolygon
-
-## Example
-
-```jldoctest
-import GeoInterface as GI, GeometryOps as GO
-
-line1 = GI.Line([(124.584961,-12.768946), (126.738281,-17.224758)])
-line2 = GI.Line([(123.354492,-15.961329), (127.22168,-14.008696)])
-GO.intersection(line1, line2)
-
-# output
-(125.58375366067547, -14.83572303404496)
-```
-"""
-intersection(geom_a, geom_b) =
-    intersection(GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
-
-"""
-    intersection(
-        ::GI.LineTrait, line_a,
-        ::GI.LineTrait, line_b,
-    )::Union{
-        ::Tuple{::Real, ::Real},
-        ::Nothing
-    }
-
-Calculates the intersection between two line segments. Return nothing if
-there isn't one.
-"""
-function intersection(::GI.LineTrait, line_a, ::GI.LineTrait, line_b)
-    # Get start and end points for both lines
-    a1 = GI.getpoint(line_a, 1)
-    a2 = GI.getpoint(line_a, 2)
-    b1 = GI.getpoint(line_b, 1)
-    b2 = GI.getpoint(line_b, 2)
-    # Determine the intersection point
-    point, fracs = _intersection_point((a1, a2), (b1, b2))
-    # Determine if intersection point is on line segments
-    if !isnothing(point) && 0 <= fracs[1] <= 1 && 0 <= fracs[2] <= 1
-        return point
-    end
-    return nothing
-end
-
-intersection(
-    trait_a::Union{GI.LineStringTrait, GI.LinearRingTrait},
-    geom_a,
-    trait_b::Union{GI.LineStringTrait, GI.LinearRingTrait},
-    geom_b,
-) = intersection_points(trait_a, geom_a, trait_b, geom_b)
-
-"""
-    intersection(
-        ::GI.PolygonTrait, poly_a,
-        ::GI.PolygonTrait, poly_b,
-    )::Union{
-        ::Vector{Vector{Tuple{::Real, ::Real}}}, # is this a good return type?
-        ::Nothing
-    }
-
-Calculates the intersection between two line segments. Return nothing if
-there isn't one.
-"""
-function intersection(::GI.PolygonTrait, poly_a, ::GI.PolygonTrait, poly_b)
-    @assert false "Polygon intersection isn't implemented yet."
-    return nothing
-end
-
-"""
-    intersection(
-        ::GI.AbstractTrait, geom_a,
-        ::GI.AbstractTrait, geom_b,
-    )::Union{
-        ::Vector{Vector{Tuple{::Real, ::Real}}}, # is this a good return type?
-        ::Nothing
-    }
-
-Calculates the intersection between two line segments. Return nothing if
-there isn't one.
-"""
-function intersection(
-    trait_a::GI.AbstractGeometryTrait, geom_a,
-    trait_b::GI.AbstractGeometryTrait, geom_b,
+#= Returns true if there is at least one intersection between edges within the
+two lists of edges. =#
+function _line_intersects(
+    edges_a::Vector{Edge},
+    edges_b::Vector{Edge}
 )
-    @assert(
-        false,
-        "Intersection between $trait_a and $trait_b isn't implemented yet.",
-    )
-    return nothing
-end
-
-"""
-    intersection_points(
-        geom_a,
-        geom_b,
-    )::Union{
-        ::Vector{::Tuple{::Real, ::Real}},
-        ::Nothing,
-    }
-
-Return a list of intersection points between two geometries. If no intersection
-point was possible given geometry extents, return nothing. If none are found,
-return an empty list.
-"""
-intersection_points(geom_a, geom_b) =
-    intersection_points(GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
-
-"""
-    intersection_points(
-        ::GI.AbstractTrait, geom_a,
-        ::GI.AbstractTrait, geom_b,
-    )::Union{
-        ::Vector{::Tuple{::Real, ::Real}},
-        ::Nothing,
-    }
-
-Calculates the list of intersection points between two geometries, inlcuding
-line segments, line strings, linear rings, polygons, and multipolygons. If no
-intersection points were possible given geometry extents, return nothing. If
-none are found, return an empty list.
-"""
-function intersection_points(::GI.AbstractTrait, a, ::GI.AbstractTrait, b)
-    # Check if the geometries extents even overlap
-    Extents.intersects(GI.extent(a), GI.extent(b)) || return nothing
-    # Create a list of edges from the two input geometries
-    edges_a, edges_b = map(sort! ∘ to_edges, (a, b))
-    npoints_a, npoints_b  = length(edges_a), length(edges_b)
-    a_closed = npoints_a > 1 && edges_a[1][1] == edges_a[end][1]
-    b_closed = npoints_b > 1 && edges_b[1][1] == edges_b[end][1]
-    if npoints_a > 0 && npoints_b > 0
-        # Initialize an empty list of points
-        T = typeof(edges_a[1][1][1]) # x-coordinate of first point in first edge
-        result = Tuple{T,T}[]
-        # Loop over pairs of edges and add any intersection points to results
-        for i in eachindex(edges_a)
-            for j in eachindex(edges_b)
-                point, fracs = _intersection_point(edges_a[i], edges_b[j])
-                if !isnothing(point)
-                    #=
-                    Determine if point is on edge (all edge endpoints excluded
-                    except for the last edge for an open geometry)
-                    =#
-                    α, β = fracs
-                    on_a_edge = (!a_closed && i == npoints_a && 0 <= α <= 1) ||
-                        (0 <= α < 1)
-                    on_b_edge = (!b_closed && j == npoints_b && 0 <= β <= 1) ||
-                        (0 <= β < 1)
-                    if on_a_edge && on_b_edge
-                        push!(result, point)
-                    end
-                end
-            end
+    # Extents.intersects(to_extent(edges_a), to_extent(edges_b)) || return false
+    for edge_a in edges_a
+        for edge_b in edges_b
+            _line_intersects(edge_a, edge_b) && return true 
         end
-        return result
     end
-    return nothing
+    return false
 end
 
-"""
-    _intersection_point(
-        (a1, a2)::Tuple,
-        (b1, b2)::Tuple,
-    )
-
-Calculates the intersection point between two lines if it exists, and as if the
-line extended to infinity, and the fractional component of each line from the
-initial end point to the intersection point.
-Inputs:
-    (a1, a2)::Tuple{Tuple{::Real, ::Real}, Tuple{::Real, ::Real}} first line
-    (b1, b2)::Tuple{Tuple{::Real, ::Real}, Tuple{::Real, ::Real}} second line
-Outputs:
-    (x, y)::Tuple{::Real, ::Real} intersection point
-    (t, u)::Tuple{::Real, ::Real} fractional length of lines to intersection
-    Both are ::Nothing if point doesn't exist!
-
-Calculation derivation can be found here:
-    https://stackoverflow.com/questions/563198/
-"""
-function _intersection_point((a1, a2)::Tuple, (b1, b2)::Tuple)
-    # First line runs from p to p + r
-    px, py = GI.x(a1), GI.y(a1)
-    rx, ry = GI.x(a2) - px, GI.y(a2) - py
-    # Second line runs from q to q + s 
-    qx, qy = GI.x(b1), GI.y(b1)
-    sx, sy = GI.x(b2) - qx, GI.y(b2) - qy
-    # Intersection will be where p + tr = q + us where 0 < t, u < 1 and
-    r_cross_s = rx * sy - ry * sx
-    if r_cross_s != 0
-        Δqp_x = qx - px
-        Δqp_y = qy - py
-        t = (Δqp_x * sy - Δqp_y * sx) / r_cross_s
-        u = (Δqp_x * ry - Δqp_y * rx) / r_cross_s
-        x = px + t * rx
-        y = py + t * ry
-        return (x, y), (t, u)
-    end
-    return nothing, nothing
-end
+# Returns true if there is at least one intersection between two edges.
+function _line_intersects(edge_a::Edge, edge_b::Edge)
+    meet_type = ExactPredicates.meet(edge_a..., edge_b...)
+    return meet_type == 0 || meet_type == 1
+end
\ No newline at end of file
diff --git a/src/utils.jl b/src/utils.jl
index b95b78172..c934ed6fc 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -94,7 +94,7 @@ function to_extent(edges::Vector{Edge})
     Extents.Extent(X=x, Y=y)
 end
 
-function to_points(xs)
+function to_points(x)
     points = Vector{TuplePoint}(undef, _npoint(x))
     _to_points!(points, x, 1)
     return points
@@ -123,4 +123,4 @@ function _to_points!(points::Vector, ::Union{AbstractCurveTrait,MultiPointTrait}
         n += 1
     end
     return n
-end
+end
\ No newline at end of file
diff --git a/test/methods/bools.jl b/test/methods/bools.jl
index 35706e2bf..134a81f7b 100644
--- a/test/methods/bools.jl
+++ b/test/methods/bools.jl
@@ -32,18 +32,4 @@ import GeometryOps as GO
 
 	@test isclockwise(l1) == true
 	@test isclockwise(l2) == false
-
-	l3 = GI.LineString([[0, 0], [3, 3], [4, 4]])
-	p1 = GI.Point([1,1])
-
-	l4 = GI.LineString([[0, 0], [3, 3]])
-	p2 = GI.Point([0, 0])
-
-	p3 = GI.Point([20, 20])
-	l5 = GI.LineString([[0, 0], [3, 3], [38.32, 5.96]])
-
-	pt = (-77, 44)
-	poly = GI.Polygon([[[-81, 41], [-81, 47], [-72, 47], [-72, 41], [-81, 41]]])
-
-
 end
diff --git a/test/methods/clipping/difference.jl b/test/methods/clipping/difference.jl
new file mode 100644
index 000000000..c113e1d64
--- /dev/null
+++ b/test/methods/clipping/difference.jl
@@ -0,0 +1,41 @@
+#=
+    compare_GO_LG_difference(p1, p2, ϵ)::Bool
+
+    Returns true if the 'difference' function from LibGEOS and 
+    GeometryOps return similar enough polygons (determined by ϵ).
+=#
+function compare_GO_LG_difference(p1, p2, ϵ)
+    GO_difference = GO.difference(p1,p2; target = GI.PolygonTrait)
+    LG_difference = LG.difference(p1,p2)
+    if isempty(GO_difference) && LG.isEmpty(LG_difference)
+        return true
+    end
+
+    if length(GO_difference)==1
+        GO_difference_poly = GO_difference[1]
+    else
+        GO_difference_poly = GI.MultiPolygon(GO_difference)
+    end
+    return LG.area(LG.difference(GO_difference_poly, LG_difference)) < ϵ
+end
+
+@testset "Difference_polygons" begin
+    # Two "regular" polygons that intersect
+    p1 = [[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]
+    p2 = [[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]
+    @test compare_GO_LG_difference(GI.Polygon([p1]), GI.Polygon([p2]), 1e-5)
+
+    # Two ugly polygons with 2 holes each
+    p3 = [[(0.0, 0.0), (5.0, 0.0), (5.0, 8.0), (0.0, 8.0), (0.0, 0.0)], [(4.0, 0.5), (4.5, 0.5), (4.5, 3.5), (4.0, 3.5), (4.0, 0.5)], [(2.0, 4.0), (4.0, 4.0), (4.0, 6.0), (2.0, 6.0), (2.0, 4.0)]]
+    p4 = [[(3.0, 1.0), (8.0, 1.0), (8.0, 7.0), (3.0, 7.0), (3.0, 5.0), (6.0, 5.0), (6.0, 3.0), (3.0, 3.0), (3.0, 1.0)], [(3.5, 5.5), (6.0, 5.5), (6.0, 6.5), (3.5, 6.5), (3.5, 5.5)], [(5.5, 1.5), (5.5, 2.5), (3.5, 2.5), (3.5, 1.5), (5.5, 1.5)]]
+    @test compare_GO_LG_difference(GI.Polygon(p3), GI.Polygon(p4), 1e-5)
+
+    # # The two polygons that intersect from the Greiner paper
+    # greiner_1 = [(0.0, 0.0), (0.0, 4.0), (7.0, 4.0), (7.0, 0.0), (0.0, 0.0)]
+    # greiner_2 = [(1.0, -3.0), (1.0, 1.0), (3.5, -1.5), (6.0, 1.0), (6.0, -3.0), (1.0, -3.0)]
+    # @test compare_GO_LG_difference(GI.Polygon([greiner_1]), GI.Polygon([greiner_2]), 1e-5)
+
+    # ugly difference test
+    pa = [[(0.0, 0.0), (8.0, 0.0), (8.0, 8.0), (0.0, 8.0), (0.0, 0.0)], [(5.0, 5.0), (7.0, 5.0), (7.0, 7.0), (5.0, 7.0), (5.0, 5.0)]]
+    pb = [[(3.0, -1.0), (10.0, -1.0), (10.0, 9.0), (3.0, 9.0), (3.0, -1.0)], [(4.0, 3.0), (5.0, 3.0), (5.0, 4.0), (4.0, 4.0), (4.0, 3.0)], [(6.0, 1.0), (7.0, 1.0), (7.0, 2.0), (6.0, 2.0), (6.0, 1.0)]]
+end
\ No newline at end of file
diff --git a/test/methods/clipping/intersection.jl b/test/methods/clipping/intersection.jl
new file mode 100644
index 000000000..37ced0560
--- /dev/null
+++ b/test/methods/clipping/intersection.jl
@@ -0,0 +1,279 @@
+#=
+    compare_GO_LG_intersection(p1, p2, ϵ)::Bool
+
+    Returns true if the 'intersection' function from LibGEOS and 
+    GeometryOps return similar enough polygons (determined by ϵ).
+=#
+function compare_GO_LG_intersection(p1, p2, ϵ)
+    GO_intersection = GO.intersection(p1,p2; target = GI.PolygonTrait)
+    LG_intersection = LG.intersection(p1,p2)
+    if isempty(GO_intersection) && LG.isEmpty(LG_intersection)
+        return true
+    end
+
+    if length(GO_intersection)==1
+        GO_intersection_poly = GO_intersection[1]
+    else
+        GO_intersection_poly = GI.MultiPolygon(GO_intersection)
+    end
+    return LG.area(LG.difference(GO_intersection_poly, LG_intersection)) < ϵ
+end
+
+@testset "Line-Line Intersection" begin
+    # Parallel lines
+    @test isempty(GO.intersection(
+        GI.Line([(0.0, 0.0), (2.5, 0.0)]),
+        GI.Line([(0.0, 1.0), (2.5, 1.0)]);
+        target = GI.PointTrait
+    ))
+    # Non-parallel lines that don't intersect
+    @test isempty(GO.intersection(
+        GI.Line([(0.0, 0.0), (2.5, 0.0)]),
+        GI.Line([(2.0, -3.0), (3.0, 0.0)]);
+        target = GI.PointTrait
+    ))
+    # Test for lines only touching at endpoint
+    l1 = GI.Line([(0.0, 0.0), (2.5, 0.0)])
+    l2 = GI.Line([(2.0, -3.0), (2.5, 0.0)])
+    @test all(GO.equals(
+        GO.intersection(
+            GI.Line([(0.0, 0.0), (2.5, 0.0)]),
+            GI.Line([(2.0, -3.0), (2.5, 0.0)]);
+            target = GI.PointTrait
+        )[1], (2.5, 0.0),
+    ))
+    # Test for lines that intersect in the middle
+    @test all(GO.equals(
+        GO.intersection(
+            GI.Line([(0.0, 0.0), (5.0, 5.0)]),
+            GI.Line([(0.0, 5.0), (5.0, 0.0)]);
+            target = GI.PointTrait
+        )[1], GI.Point((2.5, 2.5)),
+    ))
+    # Single element line strings crossing over each other
+    l1 = LG.LineString([[5.5, 7.2], [11.2, 12.7]])
+    l2 = LG.LineString([[4.3, 13.3], [9.6, 8.1]])
+    go_inter = GO.intersection(l1, l2; target = GI.PointTrait)
+    lg_inter = LG.intersection(l1, l2)
+    @test GI.x(go_inter[1]) .≈ GI.x(lg_inter)
+    @test GI.y(go_inter[1]) .≈ GI.y(lg_inter)
+    # Multi-element line strings crossing over on vertex
+    l1 = LG.LineString([[0.0, 0.0], [2.5, 0.0], [5.0, 0.0]])
+    l2 = LG.LineString([[2.0, -3.0], [3.0, 0.0], [4.0, 3.0]])
+    go_inter = GO.intersection(l1, l2; target = GI.PointTrait)
+    lg_inter = LG.intersection(l1, l2)
+    @test GI.x(go_inter[1]) .≈ GI.x(lg_inter)
+    @test GI.y(go_inter[1]) .≈ GI.y(lg_inter)
+    # Multi-element line strings crossing over with multiple intersections
+    l1 = LG.LineString([[0.0, -1.0], [1.0, 1.0], [2.0, -1.0], [3.0, 1.0]])
+    l2 = LG.LineString([[0.0, 0.0], [1.0, 0.0], [3.0, 0.0]])
+    go_inter = GO.intersection(l1, l2; target = GI.PointTrait)
+    lg_inter = LG.intersection(l1, l2)
+    @test length(go_inter) == 3
+    @test issetequal(Set(GO._tuple_point.(go_inter)), Set(GO._tuple_point.(GI.getpoint(lg_inter))))
+    # Line strings far apart so extents don't overlap
+    @test isempty(GO.intersection(
+        LG.LineString([[100.0, 0.0], [101.0, 0.0], [103.0, 0.0]]),
+        LG.LineString([[0.0, 0.0], [1.0, 0.0], [3.0, 0.0]]);
+        target = GI.PointTrait
+    ))
+    # Line strings close together that don't overlap
+    @test isempty(GO.intersection(
+        LG.LineString([[3.0, 0.25], [5.0, 0.25], [7.0, 0.25]]),
+        LG.LineString([[0.0, 0.0], [5.0, 10.0], [10.0, 0.0]]);
+        target = GI.PointTrait
+    ))
+    # Closed linear ring with open line string
+    r1 = LG.LinearRing([[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]])
+    l2 = LG.LineString([[0.0, -2.0], [12.0, 10.0],])
+    go_inter = GO.intersection(r1, l2; target = GI.PointTrait)
+    lg_inter = LG.intersection(r1, l2)
+    @test length(go_inter) == 2
+    @test issetequal(Set(GO._tuple_point.(go_inter)), Set(GO._tuple_point.(GI.getpoint(lg_inter))))
+    # Closed linear ring with closed linear ring
+    r1 = LG.LinearRing([[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]])
+    r2 = LG.LineString([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]])
+    go_inter = GO.intersection(r1, r2; target = GI.PointTrait)
+    lg_inter = LG.intersection(r1, r2)
+    @test length(go_inter) == 2
+    @test issetequal(Set(GO._tuple_point.(go_inter)), Set(GO._tuple_point.(GI.getpoint(lg_inter))))
+end
+
+@testset "Intersection Points" begin
+    # Two polygons that intersect
+    @test all(GO.equals.(
+        GO.intersection_points(
+            LG.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]]),
+            LG.Polygon([[[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]]),
+        ), [(6.5, 3.5), (6.5, -3.5)]))
+    # Two polygons that don't intersect
+    @test isempty(GO.intersection_points(
+        LG.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]]),
+        LG.Polygon([[[13.0, 0.0], [18.0, 5.0], [23.0, 0.0], [18.0, -5.0], [13.0, 0.0]]]),
+    ))
+    # Polygon that intersects with linestring
+    @test all(GO.equals.(
+        GO.intersection_points(
+            LG.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]]),
+            LG.LineString([[0.0, 0.0], [10.0, 0.0]]),
+        ),[(0.0, 0.0), (10.0, 0.0)]))
+
+    # Polygon with a hole, line through polygon and hole
+    @test all(GO.equals.(
+        GO.intersection_points(
+            LG.Polygon([
+                [[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]],
+                [[2.0, -1.0], [2.0, 1.0], [3.0, 1.0], [3.0, -1.0], [2.0, -1.0]]
+            ]),
+            LG.LineString([[0.0, 0.0], [10.0, 0.0]]),
+        ), [(0.0, 0.0), (2.0, 0.0), (3.0, 0.0), (10.0, 0.0)]))
+
+    # Polygon with a hole, line only within the hole
+    @test isempty(GO.intersection_points(
+        LG.Polygon([
+            [[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]],
+            [[2.0, -1.0], [2.0, 1.0], [3.0, 1.0], [3.0, -1.0], [2.0, -1.0]]
+        ]),
+        LG.LineString([[2.25, 0.0], [2.75, 0.0]]),
+    ))
+end
+
+@testset "Polygon-Polygon Intersection" begin
+    # Two "regular" polygons that intersect
+    p1 = [[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]
+    p2 = [[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]
+    @test compare_GO_LG_intersection(GI.Polygon([p1]), GI.Polygon([p2]), 1e-5)
+
+    # Polygons made with low spikiniess so that they are convex that intersect
+    poly_1 = [(4.526700198111509, 3.4853728532584696), (2.630732683726619, -4.126134282323841),
+     (-0.7638522032421201, -4.418734350277446), (-4.367920073785058, -0.2962672719707883),
+     (4.526700198111509, 3.4853728532584696)]
+
+    poly_2 = [(5.895141140952208, -0.8095078714426418), (2.8634927670695283, -4.625511746720306),
+     (-1.790623183259246, -4.138092164660989), (-3.9856656502985843, -0.5275687876429914),
+     (-2.554809853598822, 3.553455552936806), (1.1865909598835922, 4.984203644564732),
+     (5.895141140952208, -0.8095078714426418)]
+     
+    @test compare_GO_LG_intersection(GI.Polygon([poly_1]), GI.Polygon([poly_2]), 1e-5)
+
+    poly_3 = [(0.5901491691115478, 5.417479903013745),
+            (5.410537020963228, 1.6848837893332593),
+            (3.0666309839349886, -3.9084006463341616),
+            (-2.514245379142719, -3.513692563749087),
+            (0.5901491691115478, 5.417479903013745)]
+
+    poly_4 = [(-0.3877415677821795, 4.820886659632285),
+            (2.4402316126286077, 4.815948978468927),
+            (5.619171316140831, 1.820304652282779),
+            (5.4834492257940335, -2.0751598463740715),
+            (0.6301127609188103, -4.9788233300733635),
+            (-3.577552839568708, -2.039090724825537),
+            (-0.3877415677821795, 4.820886659632285)]
+
+    @test compare_GO_LG_intersection(GI.Polygon([poly_3]), GI.Polygon([poly_4]), 1e-5)
+
+    # Highly irregular convex polygons that intersect
+    poly_5 = [(2.5404227081738795, 0.5995497066446837),
+    (0.7215593458353178, -1.5811392990170074),
+    (-0.6792714151561866, -1.0909218208298457),
+    (-0.5721092724334685, 2.0387826195734795),
+    (0.0011462224659918308, 2.3273077404755487),
+    (2.5404227081738795, 0.5995497066446837)]
+
+    poly_6 = [(3.2022522653586183, -4.4613815131276615),
+    (-1.0482425878695998, -4.579816661708281),
+    (-3.630239248625253, 2.0443933767558677),
+    (-2.6164940041615927, 3.4708149011067224),
+    (1.725945294696213, 4.954192017601067),
+    (3.2022522653586183, -4.4613815131276615)]
+
+    @test compare_GO_LG_intersection(GI.Polygon([poly_5]), GI.Polygon([poly_6]), 1e-5)
+
+     # Concave polygons that intersect
+    poly_7 = [(1.2938349167338743, -3.175128530227131),
+    (-2.073885870841754, -1.6247711001754137),
+    (-5.787437985975053, 0.06570713422599561),
+    (-2.1308128111898093, 5.426689675486368),
+    (2.3058074184797244, 6.926652158268195),
+    (1.2938349167338743, -3.175128530227131)]
+
+    poly_8 = [(-2.1902469793743924, -1.9576242117579579),
+    (-4.726006206053999, 1.3907098941556428),
+    (-3.165301985923147, 2.847612825874245),
+    (-2.5529280962099428, 4.395492123980911),
+    (0.5677700216973937, 6.344638314896882),
+    (3.982554842356183, 4.853519613487035),
+    (5.251193948893394, 0.9343031382106848),
+    (5.53045582244555, -3.0101433691361734),
+    (-2.1902469793743924, -1.9576242117579579)]
+
+    @test compare_GO_LG_intersection(GI.Polygon([poly_7]), GI.Polygon([poly_8]), 1e-5)
+
+    poly_9 = [(5.249356078602074, -2.8345817731726015),
+    (2.3352302336587907, -3.8552311330323303),
+    (-2.5682755307722944, -3.2242349917570725),
+    (-5.090361026588791, 3.1787748367040436),
+    (2.510187688737385, 6.80232133791085),
+    (5.249356078602074, -2.8345817731726015)]
+
+    poly_10 = [(0.9429816692520585, -4.059373565182292),
+    (-1.9709301763568616, -2.2232906176621063),
+    (-3.916179758100803, 0.11584395040497442),
+    (-4.029910114454336, 2.544556062019178),
+    (0.03076291665013231, 5.265930727801515),
+    (3.9227716264243835, 6.050009375494719),
+    (5.423174791323876, 4.332978069727759),
+    (6.111111791385024, 0.660511002979265),
+    (0.9429816692520585, -4.059373565182292)]
+
+    @test compare_GO_LG_intersection(GI.Polygon([poly_9]), GI.Polygon([poly_10]), 1e-5)
+
+    # The two polygons that intersect from the Greiner paper
+    greiner_1 = [(0.0, 0.0), (0.0, 4.0), (7.0, 4.0), (7.0, 0.0), (0.0, 0.0)]
+    greiner_2 = [(1.0, -3.0), (1.0, 1.0), (3.5, -1.5), (6.0, 1.0), (6.0, -3.0), (1.0, -3.0)]
+    @test compare_GO_LG_intersection(GI.Polygon([greiner_1]), GI.Polygon([greiner_2]), 1e-5)
+
+    # Two polygons with two separate polygons as their intersection
+    poly_11 = [(1.0, 1.0), (4.0, 1.0), (4.0, 2.0), (2.0, 2.0), (2.0, 3.0), (4.0, 3.0), (4.0, 4.0), (1.0, 4.0), (1.0, 1.0)]
+    poly_12 = [(3.0, 0.0), (5.0, 0.0), (5.0, 5.0), (3.0, 5.0), (3.0, 0.0)]
+    @test compare_GO_LG_intersection(GI.Polygon([poly_11]), GI.Polygon([poly_12]), 1e-5)
+
+    # Two polygons with four separate polygons as their intersection
+    poly_13 = [(1.0, 1.0), (4.0, 1.0), (4.0, 2.0), (2.0, 2.0), (2.0, 3.0), (4.0, 3.0), (4.0, 4.0), (2.0, 4.0), (2.0, 5.0),
+               (4.0, 5.0), (4.0, 6.0), (2.0, 6.0), (2.0, 7.0), (4.0, 7.0), (4.0, 8.0), (1.0, 8.0), (1.0, 1.0)]
+    poly_14 = [(3.0, 0.0), (5.0, 0.0), (5.0, 9.0), (3.0, 9.0), (3.0, 0.0)]
+    @test compare_GO_LG_intersection(GI.Polygon([poly_13]), GI.Polygon([poly_14]), 1e-5)
+
+    # Polygon completely inside other polygon (no holes)
+    poly_in = [(1.0, 1.0), (2.0, 1.0), (2.0, 2.0), (1.0, 2.0), (1.0, 1.0)]
+    poly_out = [(0.0, 0.0), (3.0, 0.0), (3.0, 3.0), (0.0, 3.0), (0.0, 0.0)]
+    @test compare_GO_LG_intersection(GI.Polygon([poly_in]), GI.Polygon([poly_out]), 1e-5)
+    
+    
+    # Intersection of convex and concave polygon
+    @test compare_GO_LG_intersection(GI.Polygon([poly_1]), GI.Polygon([poly_10]), 1e-5)
+    @test compare_GO_LG_intersection(GI.Polygon([poly_5]), GI.Polygon([poly_7]), 1e-5)
+    @test compare_GO_LG_intersection(GI.Polygon([poly_4]), GI.Polygon([poly_9]), 1e-5)
+
+
+    # Intersection polygon with a hole, but other polygon no hole. Hole completely contained in other polygon
+    p_hole = [[(0.0, 0.0), (4.0, 0.0), (4.0, 3.0), (0.0, 3.0), (0.0, 0.0)], [(2.0, 1.0), (3.0, 1.0), (3.0, 2.0), (2.0, 2.0), (2.0, 1.0)]]
+    p_nohole = [[(1.0, -1.0), (1.0, 4.0), (5.0, 4.0), (5.0, -1.0), (1.0, -1.0)]]
+
+    # testing union of poly with hole (union relies on difference)
+    p_hole = [[(0.0, 0.0), (4.0, 0.0), (4.0, 4.0), (0.0, 4.0), (0.0, 0.0)], [(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]]
+    p_nohole = [[(2.0, -1.0), (5.0, -1.0), (5.0, 5.0), (2.0, 5.0), (2.0, -1.0)]]
+
+
+    # splitting a polygon dif test
+    pbig = [[(0.0, 0.0), (3.0, 0.0), (3.0, 3.0), (0.0, 3.0), (0.0, 0.0)]]
+    psplit = [[(1.0, -1.0), (2.0, -1.0), (2.0, 4.0), (1.0, 4.0), (1.0, -1.0)]]
+
+    # Two donut shaped polygons with intersecting holes
+    tall_donut = [[(0.0, 0.0), (6.0, 0.0), (6.0, 7.0), (0.0, 7.0), (0.0, 0.0)], [(1.0, 1.0), (5.0, 1.0), (5.0, 6.0), (1.0, 6.0), (1.0, 1.0)]]
+    wide_donut = [[(2.0, 2.0), (8.0, 2.0), (8.0, 5.0), (2.0, 5.0), (2.0, 2.0)], [(3.0, 3.0), (7.0, 3.0), (7.0, 4.0), (3.0, 4.0), (3.0, 3.0)]]
+
+    # Two ugly polygons with 2 holes each
+    p1 = [[(0.0, 0.0), (5.0, 0.0), (5.0, 8.0), (0.0, 8.0), (0.0, 0.0)], [(4.0, 0.5), (4.5, 0.5), (4.5, 3.5), (4.0, 3.5), (4.0, 0.5)], [(2.0, 4.0), (4.0, 4.0), (4.0, 6.0), (2.0, 6.0), (2.0, 4.0)]]
+    p2 = [[(3.0, 1.0), (8.0, 1.0), (8.0, 7.0), (3.0, 7.0), (3.0, 5.0), (6.0, 5.0), (6.0, 3.0), (3.0, 3.0), (3.0, 1.0)], [(3.5, 5.5), (6.0, 5.5), (6.0, 6.5), (3.5, 6.5), (3.5, 5.5)], [(5.5, 1.5), (5.5, 2.5), (3.5, 2.5), (3.5, 1.5), (5.5, 1.5)]]
+end
\ No newline at end of file
diff --git a/test/methods/clipping/union.jl b/test/methods/clipping/union.jl
new file mode 100644
index 000000000..1b7f96036
--- /dev/null
+++ b/test/methods/clipping/union.jl
@@ -0,0 +1,43 @@
+#=
+    compare_GO_LG_union(p1, p2, ϵ)::Bool
+
+    Returns true if the 'union' function from LibGEOS and 
+    GeometryOps return similar enough polygons (determined by ϵ).
+=#
+function compare_GO_LG_union(p1, p2, ϵ)
+    GO_union = GO.union(p1,p2; target = GI.PolygonTrait)
+    LG_union = LG.union(p1,p2)
+    if isempty(GO_union) && LG.isEmpty(LG_union)
+        return true
+    end
+
+    if length(GO_union)==1
+        GO_union_poly = GO_union[1]
+    else
+        GO_union_poly = GI.MultiPolygon(GO_union)
+    end
+
+    return LG.area(LG.difference(GO_union_poly, LG_union)) < ϵ
+end
+
+@testset "Union_polygons" begin
+    # Two "regular" polygons that intersect
+    p1 = [[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]]
+    p2 = [[[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]]
+    @test compare_GO_LG_union(GI.Polygon(p1), GI.Polygon(p2), 1e-5)
+
+    # Two ugly polygons with 2 holes each
+    p1 = [[(0.0, 0.0), (5.0, 0.0), (5.0, 8.0), (0.0, 8.0), (0.0, 0.0)], [(4.0, 0.5), (4.5, 0.5), (4.5, 3.5), (4.0, 3.5), (4.0, 0.5)], [(2.0, 4.0), (4.0, 4.0), (4.0, 6.0), (2.0, 6.0), (2.0, 4.0)]]
+    p2 = [[(3.0, 1.0), (8.0, 1.0), (8.0, 7.0), (3.0, 7.0), (3.0, 5.0), (6.0, 5.0), (6.0, 3.0), (3.0, 3.0), (3.0, 1.0)], [(3.5, 5.5), (6.0, 5.5), (6.0, 6.5), (3.5, 6.5), (3.5, 5.5)], [(5.5, 1.5), (5.5, 2.5), (3.5, 2.5), (3.5, 1.5), (5.5, 1.5)]]
+    @test compare_GO_LG_union(GI.Polygon(p1), GI.Polygon(p2), 1e-5)
+
+    # Union test when the two polygons are disjoint and each have one hole (two disjoint square donuts)
+    p1 = [[(0.0, 0.0), (3.0, 0.0), (3.0, 3.0), (0.0, 3.0), (0.0, 0.0)], [(1.0, 1.0), (2.0, 1.0), (2.0, 2.0), (1.0, 2.0), (1.0, 1.0)]]
+    p2 = [[(5.0, 0.0), (8.0, 0.0), (8.0, 3.0), (5.0, 3.0), (5.0, 0.0)], [(6.0, 1.0), (7.0, 1.0), (7.0, 2.0), (7.0, 1.0), (6.0, 1.0)]]
+    @test compare_GO_LG_union(GI.Polygon(p1), GI.Polygon(p2), 1e-5)
+
+    # The two polygons that intersect from the Greiner paper
+    greiner_1 = [(0.0, 0.0), (0.0, 4.0), (7.0, 4.0), (7.0, 0.0), (0.0, 0.0)]
+    greiner_2 = [(1.0, -3.0), (1.0, 1.0), (3.5, -1.5), (6.0, 1.0), (6.0, -3.0), (1.0, -3.0)]
+    @test compare_GO_LG_union(GI.Polygon([greiner_1]), GI.Polygon([greiner_2]), 1e-5)
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 1bb337d63..017239492 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -20,9 +20,13 @@ const GO = GeometryOps
     @testset "Barycentric coordinate operations" begin include("methods/barycentric.jl") end
     @testset "Bools" begin include("methods/bools.jl") end
     @testset "Centroid" begin include("methods/centroid.jl") end
+    @testset "DE-9IM Geom Relations" begin include("methods/geom_relations.jl") end
     @testset "Distance" begin include("methods/distance.jl") end
     @testset "Equals" begin include("methods/equals.jl") end
-    @testset "DE-9IM Geom Relations" begin include("methods/geom_relations.jl") end
+    # Clipping
+    @testset "Difference" begin include("methods/clipping/difference.jl") end
+    @testset "Intersection" begin include("methods/clipping/intersection.jl") end
+    @testset "Union" begin include("methods/clipping/union.jl") end
     # # Transformations
     @testset "Embed Extent" begin include("transformations/extent.jl") end
     @testset "Reproject" begin include("transformations/reproject.jl") end