diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index 588c404ed..dc52c7d17 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -54,6 +54,7 @@ include("transformations/tuples.jl")
 include("transformations/transform.jl")
 include("transformations/correction/geometry_correction.jl")
 include("transformations/correction/closed_ring.jl")
+include("transformations/correction/intersecting_polygons.jl")
 
 # Import all names from GeoInterface and Extents, so users can do `GO.extent` or `GO.trait`.
 for name in names(GeoInterface)
diff --git a/src/methods/clipping/difference.jl b/src/methods/clipping/difference.jl
index 83f945c92..f7477e0d2 100644
--- a/src/methods/clipping/difference.jl
+++ b/src/methods/clipping/difference.jl
@@ -3,13 +3,19 @@ export difference
 
 
 """
-    difference(geom_a, geom_b, [T::Type]; target::Type)
+    difference(geom_a, geom_b, [T::Type]; target::Type, fix_multipoly = UnionIntersectingPolygons())
 
 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. 
+provide a float type that they would like the points of returned geometries to be. If the
+user is taking a intersection involving one or more multipolygons, and the multipolygon
+might be comprised of polygons that intersect, if `fix_multipoly` is set to an
+`IntersectingPolygons` correction (the default is `UnionIntersectingPolygons()`), then the
+needed multipolygons will be fixed to be valid before performing the intersection to ensure
+a correct answer. Only set `fix_multipoly` to false if you know that the multipolygons are
+valid, as it will avoid unneeded computation. 
 
 ## Example 
 
@@ -27,9 +33,9 @@ GI.coordinates.(diff_poly)
 ```
 """
 function difference(
-    geom_a, geom_b, ::Type{T} = Float64; target=nothing,
+    geom_a, geom_b, ::Type{T} = Float64; target=nothing, kwargs...,
 ) where {T<:AbstractFloat}
-    return _difference(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
+    return _difference(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b; kwargs...)
 end
 
 #= The 'difference' function returns the difference of two polygons as a list of polygons.
@@ -38,7 +44,8 @@ polygons," by Greiner and Hormann (1998). DOI: https://doi.org/10.1145/274363.27
 function _difference(
     ::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.PolygonTrait, poly_a,
-    ::GI.PolygonTrait, poly_b,
+    ::GI.PolygonTrait, poly_b;
+    kwargs...
 ) where T
     # Get the exterior of the polygons
     ext_a = GI.getexterior(poly_a)
@@ -99,6 +106,65 @@ tracing a_list, else step backwards, where x is the entry/exit status and y is a
 that is true if we are on a_list and false if we are on b_list. =#
 _diff_step(x, y) = (x ⊻ y) ? 1 : (-1)
 
+#= Polygon with multipolygon difference - note that all intersection regions between
+`poly_a` and any of the sub-polygons of `multipoly_b` are removed from `poly_a`. =#
+function _difference(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.PolygonTrait, poly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    kwargs...,
+) where T
+    polys = [tuples(poly_a, T)]
+    for poly_b in GI.getpolygon(multipoly_b)
+        isempty(polys) && break
+        polys = mapreduce(p -> difference(p, poly_b; target), append!, polys)
+    end
+    return polys
+end
+
+#= Multipolygon with polygon difference - note that all intersection regions between
+sub-polygons of `multipoly_a` and `poly_b` will be removed from the corresponding
+sub-polygon. Unless specified with `fix_multipoly = nothing`, `multipolygon_a` will be
+validated using the given (default is `UnionIntersectingPolygons()`) correction. =#
+function _difference(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.PolygonTrait, poly_b;
+    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)
+    end
+    polys = Vector{_get_poly_type(T)}()
+    sizehint!(polys, GI.npolygon(multipoly_a))
+    for poly_a in GI.getpolygon(multipoly_a)
+        append!(polys, difference(poly_a, poly_b; target))
+    end
+    return polys
+end
+
+#= Multipolygon with multipolygon difference - note that all intersection regions between
+sub-polygons of `multipoly_a` and sub-polygons of `multipoly_b` will be removed from the
+corresponding sub-polygon of `multipoly_a`. Unless specified with `fix_multipoly = nothing`,
+`multipolygon_a` will be validated using the given (defauly is `UnionIntersectingPolygons()`)
+correction. =#
+function _difference(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    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)
+        fix_multipoly = nothing
+    end
+    local polys
+    for (i, poly_b) in enumerate(GI.getpolygon(multipoly_b))
+        polys = difference(i == 1 ? multipoly_a : GI.MultiPolygon(polys), poly_b; target, fix_multipoly)
+    end
+    return polys
+end
+
 # Many type and target combos aren't implemented
 function _difference(
     ::TraitTarget{Target}, ::Type{T},
diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl
index 2387af166..832e3f787 100644
--- a/src/methods/clipping/intersection.jl
+++ b/src/methods/clipping/intersection.jl
@@ -10,13 +10,19 @@ Enum for the orientation of a line with respect to a curve. A line can be
 @enum LineOrientation line_cross=1 line_hinge=2 line_over=3 line_out=4
 
 """
-    intersection(geom_a, geom_b, [T::Type]; target::Type)
+    intersection(geom_a, geom_b, [T::Type]; target::Type, fix_multipoly = UnionIntersectingPolygons())
 
 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 `target` 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. 
+provide a float type that they would like the points of returned geometries to be. If the
+user is taking a intersection involving one or more multipolygons, and the multipolygon
+might be comprised of polygons that intersect, if `fix_multipoly` is set to an
+`IntersectingPolygons` correction (the default is `UnionIntersectingPolygons()`), then the
+needed multipolygons will be fixed to be valid before performing the intersection to ensure
+a correct answer. Only set `fix_multipoly` to nothing if you know that the multipolygons are
+valid, as it will avoid unneeded computation. 
 
 ## Example
 
@@ -34,16 +40,17 @@ GI.coordinates.(inter_points)
 ```
 """
 function intersection(
-    geom_a, geom_b, ::Type{T}=Float64; target=nothing,
+    geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...,
 ) where {T<:AbstractFloat}
-    return _intersection(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
+    return _intersection(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b; kwargs...)
 end
 
 # Curve-Curve Intersections with target Point
 _intersection(
     ::TraitTarget{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,
+    trait_b::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, geom_b;
+    kwargs...,
 ) where T = _intersection_points(T, trait_a, geom_a, trait_b, geom_b)
 
 #= Polygon-Polygon Intersections with target Polygon
@@ -53,7 +60,8 @@ DOI: https://doi.org/10.1145/274363.274364 =#
 function _intersection(
     ::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.PolygonTrait, poly_a,
-    ::GI.PolygonTrait, poly_b,
+    ::GI.PolygonTrait, poly_b;
+    kwargs...,
 ) where {T}
     # First we get the exteriors of 'poly_a' and 'poly_b'
     ext_a = GI.getexterior(poly_a)
@@ -98,11 +106,65 @@ _inter_delay_bounce_f(x, _) = x
 point, else step backwards where x is the entry/exit status. =#
 _inter_step(x, _) =  x ? 1 : (-1)
 
+#= Polygon with multipolygon intersection - note that all intersection regions between
+`poly_a` and any of the sub-polygons of `multipoly_b` are counted as intersection polygons.
+Unless specified with `fix_multipoly = nothing`, `multipolygon_b` will be validated using
+the given (default is `UnionIntersectingPolygons()`) correction. =#
+function _intersection(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.PolygonTrait, poly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    fix_multipoly = UnionIntersectingPolygons(), kwargs...,
+) where T
+    if !isnothing(fix_multipoly) # Fix multipoly_b to prevent duplicated intersection regions
+        multipoly_b = fix_multipoly(multipoly_b)
+    end
+    polys = Vector{_get_poly_type(T)}()
+    for poly_b in GI.getpolygon(multipoly_b)
+        append!(polys, intersection(poly_a, poly_b; target))
+    end
+    return polys
+end
+
+#= Multipolygon with polygon intersection is equivalent to taking the intersection of the
+poylgon with the multipolygon and thus simply switches the order of operations and calls the
+above method. =#
+_intersection(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.PolygonTrait, poly_b;
+    kwargs...,
+) where T = intersection(poly_b, multipoly_a; 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
+as intersection polygons. Unless specified with `fix_multipoly = nothing`, both 
+`multipolygon_a` and `multipolygon_b` will be validated using the given (default is
+`UnionIntersectingPolygons()`) correction. =#
+function _intersection(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    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)
+        fix_multipoly = nothing
+    end
+    polys = Vector{_get_poly_type(T)}()
+    for poly_a in GI.getpolygon(multipoly_a)
+        append!(polys, intersection(poly_a, multipoly_b; target, fix_multipoly))
+    end
+    return polys
+end
+
 # Many type and target combos aren't implemented
 function _intersection(
     ::TraitTarget{Target}, ::Type{T},
     trait_a::GI.AbstractTrait, geom_a,
-    trait_b::GI.AbstractTrait, geom_b,
+    trait_b::GI.AbstractTrait, geom_b;
+    kwargs...,
 ) where {Target, T}
     @assert(
         false,
diff --git a/src/methods/clipping/union.jl b/src/methods/clipping/union.jl
index b09749442..e390ef2f1 100644
--- a/src/methods/clipping/union.jl
+++ b/src/methods/clipping/union.jl
@@ -2,13 +2,19 @@
 export union
 
 """
-    union(geom_a, geom_b, [::Type{T}]; target::Type)
+    union(geom_a, geom_b, [::Type{T}]; target::Type, fix_multipoly = UnionIntersectingPolygons())
 
 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. 
+provide a float type 'T' that they would like the points of returned geometries to be. If
+the user is taking a intersection involving one or more multipolygons, and the multipolygon
+might be comprised of polygons that intersect, if `fix_multipoly` is set to an
+`IntersectingPolygons` correction (the default is `UnionIntersectingPolygons()`), then the
+needed multipolygons will be fixed to be valid before performing the intersection to ensure
+a correct answer. Only set `fix_multipoly` to false if you know that the multipolygons are
+valid, as it will avoid unneeded computation. 
     
 Calculates the union between two polygons.
 ## Example
@@ -27,9 +33,9 @@ GI.coordinates.(union_poly)
 ```
 """
 function union(
-    geom_a, geom_b, ::Type{T}=Float64; target=nothing,
+    geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...
 ) where {T<:AbstractFloat}
-    _union(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
+    _union(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b; kwargs...)
 end
 
 #= This 'union' implementation returns the union of two polygons. The algorithm to determine
@@ -38,7 +44,8 @@ Hormann (1998). DOI: https://doi.org/10.1145/274363.274364 =#
 function _union(
     ::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.PolygonTrait, poly_a,
-    ::GI.PolygonTrait, poly_b,
+    ::GI.PolygonTrait, poly_b;
+    kwargs...,
 ) where T
     # First, I get the exteriors of the two polygons
     ext_a = GI.getexterior(poly_a)
@@ -196,12 +203,75 @@ function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly)
     return
 end
 
+#= Polygon with multipolygon union - note that all sub-polygons of `multipoly_b` will be
+included, unioning these sub-polygons with `poly_a` where they intersect. Unless specified
+with `fix_multipoly = nothing`, `multipolygon_b` will be validated using the given (default
+is `UnionIntersectingPolygons()`) correction. =#
+function _union(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.PolygonTrait, poly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    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)
+    end
+    polys = [tuples(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)
+            if length(new_polys) > 1 # case where they intersect by just one point
+                push!(polys, tuples(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))
+        end
+    end
+    return polys
+end
+
+#= Multipolygon with polygon union is equivalent to taking the union of the poylgon with the
+multipolygon and thus simply switches the order of operations and calls the above method. =#
+_union(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.PolygonTrait, poly_b;
+    kwargs...,
+) where T = union(poly_b, multipoly_a; 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
+intersections. Unless specified with `fix_multipoly = nothing`, `multipolygon_b` will be
+validated using the given (default is `UnionIntersectingPolygons()`) correction. =#
+function _union(
+    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    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)
+        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)
+        multipolys = GI.MultiPolygon(polys)
+    end
+    return polys
+end
 
 # Many type and target combos aren't implemented
 function _union(
     ::TraitTarget{Target}, ::Type{T},
     trait_a::GI.AbstractTrait, geom_a,
-    trait_b::GI.AbstractTrait, geom_b,
+    trait_b::GI.AbstractTrait, geom_b;
+    kwargs...
 ) where {Target,T}
     throw(ArgumentError("Union between $trait_a and $trait_b with target $Target isn't implemented yet."))
     return nothing
diff --git a/src/transformations/correction/intersecting_polygons.jl b/src/transformations/correction/intersecting_polygons.jl
new file mode 100644
index 000000000..dc5888a7a
--- /dev/null
+++ b/src/transformations/correction/intersecting_polygons.jl
@@ -0,0 +1,71 @@
+# # Intersecting Polygons
+
+export UnionIntersectingPolygons
+
+#=
+If the sub-polygons of a multipolygon are intersecting, this makes them invalid according to
+specification. Each sub-polygon of a multipolygon being disjoint (other than by a single
+point) is a requirment for a valid multipolygon. However, different libraries may achieve
+this in different ways. 
+
+For example, taking the union of all sub-polygons of a multipolygon will create a new
+multipolygon where each sub-polygon is disjoint. This can be done with the
+`UnionIntersectingPolygons` correction.
+
+The reason this operates on a multipolygon level is that it is easy for users to mistakenly
+create multipolygon's that overlap, which can then be detrimental to polygon clipping
+performance and even create wrong answers.
+=#
+
+# ## Example
+#=
+
+Multipolygon providers may not check that the polygons making up their multipolygons do not
+intersect, which makes them invalid according to the specification.
+
+For example, the following multipolygon is not valid:
+
+```@example union-multipoly
+import GeoInterface as GI
+polygon = GI.Polygon([[(0.0, 0.0), (3.0, 0.0), (3.0, 3.0), (0.0, 3.0), (0.0, 0.0)]])
+multipolygon = GI.MultiPolygon([polygon, polygon])
+```
+
+given that the two sub-polygons are the exact same shape.
+
+```@example union-multipoly
+import GeometryOps as GO
+GO.fix(multipolygon, corrections = [GO.UnionIntersectingPolygons()])
+```
+
+You can see that the the multipolygon now only contains one sub-polygon, rather than the two
+identical ones provided.
+=#
+
+# ## Implementation
+
+"""
+    UnionIntersectingPolygons() <: GeometryCorrection
+
+This correction ensures that the polygon's included in a multipolygon aren't intersecting.
+If any polygon's are intersecting, they will be combined through the union operation to
+create a unique set of disjoint (other than potentially connections by a single point)
+polygons covering the same area.
+
+See also [`GeometryCorrection`](@ref).
+"""
+struct UnionIntersectingPolygons <: GeometryCorrection end
+
+application_level(::UnionIntersectingPolygons) = GI.MultiPolygonTrait
+
+function (::UnionIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly)
+    union_multipoly = if GI.npolygon(multipoly) > 1
+        # Combine any sub-polygons that intersect
+        first_poly = GI.getpolygon(multipoly, 1)
+        exclude_first_poly = GI.MultiPolygon(collect(Iterators.drop(GI.getpolygon(multipoly), 1)))
+        GI.MultiPolygon(union(first_poly, exclude_first_poly; target = GI.PolygonTrait(), fix_multipoly = nothing))
+    else
+        tuples(multipoly)
+    end
+    return union_multipoly
+end
\ No newline at end of file
diff --git a/test/methods/clipping/polygon_clipping.jl b/test/methods/clipping/polygon_clipping.jl
index 85ac47f90..afc18c538 100644
--- a/test/methods/clipping/polygon_clipping.jl
+++ b/test/methods/clipping/polygon_clipping.jl
@@ -103,6 +103,11 @@ p54 = GI.Polygon([[(2.5, 2.5), (2.5, 7.5), (7.5, 7.5), (7.5, 2.5), (2.5, 2.5)],
 p55 = GI.Polygon([[(5.0, 0.25), (5.0, 5.0), (9.5, 5.0), (9.5, 0.25), (5.0, 0.25)],
     [(6.0, 3.0), (6.0, 4.0), (7.0, 4.0), (7.0, 3.0), (6.0, 3.0)],
     [(7.5, 0.5), (7.5, 2.5), (9.25, 2.5), (9.25, 0.5), (7.5, 0.5)]])
+
+mp1 = GI.MultiPolygon([p1])
+mp2 = GI.MultiPolygon([p2])
+mp3 = GI.MultiPolygon([p34, p36])
+mp4 = GI.MultiPolygon([p33, p35])
     
 test_pairs = [
     (p1, p1, "p1", "p1", "Same polygon"),
@@ -143,6 +148,16 @@ test_pairs = [
     (p52, p53, "p52", "p53", "Polygon with two holes completely inside of one (of three) holes of other polygon"),
     (p52, p54, "p52", "p54", "Polygon with two holes has exterior equal to one (of three) holes of other polygon"),
     (p52, p55, "p52", "p55", "Polygon within another polygon, with intersecting and disjoint holes"),
+    (p1, mp2, "p1", "mp2", "Polygon overlapping with mutlipolygon with one sub-polygon"),
+    (mp1, p2, "mp1", "p2", "Polygon overlapping with mutlipolygon with one sub-polygon"),
+    (mp1, mp2, "mp1", "mp2", "Two overlapping multipolygons, which each have one sub-polygon"),
+    (mp3, p33, "mp3", "p33", "Mutipolygon, with two sub-polygon, both of which overlap with other polygon"),
+    (p33, mp3, "p33", "mp3", "Polygon overlaps with both sub-polygons of a multipolygon"),
+    (mp3, p34, "mp3", "p34", "Multipolygon, with two sub-polygon, overlaps with polygon equal to sub-polygon"),
+    (p34, mp3, "p34", "mp3", "Polygon overlaps with multipolygon, where on of the sub-polygons is equivalent"),
+    (mp3, p35, "mp3", "p35", "Mulitpolygon where just one sub-polygon touches other polygon"),
+    (p35, mp3, "p35", "mp3", "Polygon that touches just one of the sub-polygons of multipolygon"),
+    (mp4, mp3, "mp4", "mp3", "Two multipolygons, which with two sub-polygons, where some of the sub-polygons intersect and some don't")
 ]
 
 const ϵ = 1e-10
diff --git a/test/runtests.jl b/test/runtests.jl
index b30ff4e51..eca987fd6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -39,5 +39,6 @@ const GO = GeometryOps
     @testset "Geometry correction" begin 
         include("transformations/correction/geometry_correction.jl")
         include("transformations/correction/closed_ring.jl") 
+        include("transformations/correction/intersecting_polygons.jl")
     end
 end
diff --git a/test/transformations/correction/intersecting_polygons.jl b/test/transformations/correction/intersecting_polygons.jl
new file mode 100644
index 000000000..24b3d97dc
--- /dev/null
+++ b/test/transformations/correction/intersecting_polygons.jl
@@ -0,0 +1,34 @@
+using Test
+
+import GeoInterface as GI
+import GeometryOps as GO
+
+p1 = GI.Polygon([[(0.0, 0.0), (3.0, 0.0), (3.0, 3.0), (0.0, 3.0), (0.0, 0.0)]])
+# (p1, p2) -> one polygon inside of the other, sharing an edge
+p2 = GI.Polygon([[(1.0, 0.0), (2.0, 0.0), (2.0, 1.0), (1.0, 1.0), (1.0, 0.0)]])
+# (p1, p3) -> polygons outside of one another, sharing an edge
+p3 = GI.Polygon([[(1.0, 0.0), (2.0, 0.0), (2.0, -1.0), (1.0, -1.0), (1.0, 0.0)]])
+# (p1, p4) -> polygons are completly disjoint (no holes)
+p4 = GI.Polygon([[(1.0, -1.0), (2.0, -1.0), (2.0, -2.0), (1.0, -2.0), (1.0, -1.0)]])
+
+mp1 = GI.MultiPolygon([p1])
+mp2 = GI.MultiPolygon([p1, p1])
+mp3 = GI.MultiPolygon([p1, p2, p3])
+mp4 = GI.MultiPolygon([p1, p4])
+
+fixed_mp1 = GO.UnionIntersectingPolygons()(mp1)
+@test GI.npolygon(fixed_mp1) == 1
+@test GO.equals(GI.getpolygon(fixed_mp1, 1), p1)
+
+fixed_mp2 = GO.UnionIntersectingPolygons()(mp2)
+@test GI.npolygon(fixed_mp2) == 1
+@test GO.equals(GI.getpolygon(fixed_mp2, 1), p1)
+
+fixed_mp3 = GO.UnionIntersectingPolygons()(mp3)
+@test GI.npolygon(fixed_mp3) == 1
+@test GO.coveredby(p1, fixed_mp3) && GO.coveredby(p2, fixed_mp3) && GO.coveredby(p3, fixed_mp3)
+
+fixed_mp4 = GO.UnionIntersectingPolygons()(mp4)
+@test GI.npolygon(fixed_mp4) == 2
+@test (GO.equals(GI.getpolygon(fixed_mp4, 1), p1) && GO.equals(GI.getpolygon(fixed_mp4, 2), p4)) ||
+    (GO.equals(GI.getpolygon(fixed_mp4, 2), p1) && GO.equals(GI.getpolygon(fixed_mp4, 1), p4))
\ No newline at end of file