diff --git a/geo/src/algorithm/bool_ops/mod.rs b/geo/src/algorithm/bool_ops/mod.rs
index 29d2ad1be..371a862c1 100644
--- a/geo/src/algorithm/bool_ops/mod.rs
+++ b/geo/src/algorithm/bool_ops/mod.rs
@@ -2,6 +2,24 @@ use geo_types::{MultiLineString, MultiPolygon};
 
 use crate::{CoordsIter, GeoFloat, GeoNum, Polygon};
 
+/// Enum to represent geometry types that can be used in boolean ops
+pub enum BopGeometry<'a, T: GeoNum> {
+    Polygon(&'a Polygon<T>),
+    MultiPolygon(&'a MultiPolygon<T>),
+}
+
+impl<'a, T: GeoNum> From<&'a Polygon<T>> for BopGeometry<'a, T> {
+    fn from(polygon: &'a Polygon<T>) -> Self {
+        BopGeometry::Polygon(polygon)
+    }
+}
+
+impl<'a, T: GeoNum> From<&'a MultiPolygon<T>> for BopGeometry<'a, T> {
+    fn from(multi_polygon: &'a MultiPolygon<T>) -> Self {
+        BopGeometry::MultiPolygon(multi_polygon)
+    }
+}
+
 /// Boolean Operations on geometry.
 ///
 /// Boolean operations are set operations on geometries considered as a subset
@@ -25,23 +43,46 @@ use crate::{CoordsIter, GeoFloat, GeoNum, Polygon};
 pub trait BooleanOps: Sized {
     type Scalar: GeoNum;
 
-    fn boolean_op(&self, other: &Self, op: OpType) -> MultiPolygon<Self::Scalar>;
-    fn intersection(&self, other: &Self) -> MultiPolygon<Self::Scalar> {
+    fn boolean_op<'a, G>(&self, other: &'a G, op: OpType) -> MultiPolygon<Self::Scalar>
+    where
+        BopGeometry<'a, Self::Scalar>: From<&'a G>,
+        <Self>::Scalar: 'a;
+
+    fn intersection<'a, G>(&self, other: &'a G) -> MultiPolygon<Self::Scalar>
+    where
+        BopGeometry<'a, Self::Scalar>: From<&'a G>,
+        <Self>::Scalar: 'a,
+    {
         self.boolean_op(other, OpType::Intersection)
     }
-    fn union(&self, other: &Self) -> MultiPolygon<Self::Scalar> {
+
+    fn union<'a, G>(&self, other: &'a G) -> MultiPolygon<Self::Scalar>
+    where
+        BopGeometry<'a, Self::Scalar>: From<&'a G>,
+        <Self>::Scalar: 'a,
+    {
         self.boolean_op(other, OpType::Union)
     }
-    fn xor(&self, other: &Self) -> MultiPolygon<Self::Scalar> {
+
+    fn xor<'a, G>(&self, other: &'a G) -> MultiPolygon<Self::Scalar>
+    where
+        BopGeometry<'a, Self::Scalar>: From<&'a G>,
+        <Self>::Scalar: 'a,
+    {
         self.boolean_op(other, OpType::Xor)
     }
-    fn difference(&self, other: &Self) -> MultiPolygon<Self::Scalar> {
+
+    fn difference<'a, G>(&self, other: &'a G) -> MultiPolygon<Self::Scalar>
+    where
+        BopGeometry<'a, Self::Scalar>: From<&'a G>,
+        <Self>::Scalar: 'a,
+    {
         self.boolean_op(other, OpType::Difference)
     }
 
     /// Clip a 1-D geometry with self.
     ///
-    /// Returns the portion of `ls` that lies within `self` (known as the set-theoeretic
+    /// Returns the portion of `ls` that lies within `self` (known as the set-theoretic
     /// intersection) if `invert` is false, and the difference (`ls - self`) otherwise.
     fn clip(
         &self,
@@ -61,11 +102,21 @@ pub enum OpType {
 impl<T: GeoFloat> BooleanOps for Polygon<T> {
     type Scalar = T;
 
-    fn boolean_op(&self, other: &Self, op: OpType) -> MultiPolygon<Self::Scalar> {
+    fn boolean_op<'a, G>(&self, other: &'a G, op: OpType) -> MultiPolygon<Self::Scalar>
+    where
+        BopGeometry<'a, Self::Scalar>: From<&'a G>,
+        <Self>::Scalar: 'a,
+    {
+        let other: BopGeometry<'a, Self::Scalar> = other.into();
         let spec = BoolOp::from(op);
         let mut bop = Proc::new(spec, self.coords_count() + other.coords_count());
         bop.add_polygon(self, 0);
-        bop.add_polygon(other, 1);
+        match other {
+            BopGeometry::Polygon(other_polygon) => bop.add_polygon(other_polygon, 1),
+            BopGeometry::MultiPolygon(other_multi_polygon) => {
+                bop.add_multi_polygon(other_multi_polygon, 1)
+            }
+        }
         bop.sweep()
     }
 
@@ -83,14 +134,25 @@ impl<T: GeoFloat> BooleanOps for Polygon<T> {
         bop.sweep()
     }
 }
+
 impl<T: GeoFloat> BooleanOps for MultiPolygon<T> {
     type Scalar = T;
 
-    fn boolean_op(&self, other: &Self, op: OpType) -> MultiPolygon<Self::Scalar> {
+    fn boolean_op<'a, G>(&self, other: &'a G, op: OpType) -> MultiPolygon<Self::Scalar>
+    where
+        BopGeometry<'a, Self::Scalar>: From<&'a G>,
+        <Self>::Scalar: 'a,
+    {
+        let other: BopGeometry<'a, Self::Scalar> = other.into();
         let spec = BoolOp::from(op);
         let mut bop = Proc::new(spec, self.coords_count() + other.coords_count());
         bop.add_multi_polygon(self, 0);
-        bop.add_multi_polygon(other, 1);
+        match other {
+            BopGeometry::Polygon(other_polygon) => bop.add_polygon(other_polygon, 1),
+            BopGeometry::MultiPolygon(other_multi_polygon) => {
+                bop.add_multi_polygon(other_multi_polygon, 1)
+            }
+        }
         bop.sweep()
     }
 
@@ -102,13 +164,20 @@ impl<T: GeoFloat> BooleanOps for MultiPolygon<T> {
         let spec = ClipOp::new(invert);
         let mut bop = Proc::new(spec, self.coords_count() + ls.coords_count());
         bop.add_multi_polygon(self, 0);
-        ls.0.iter().enumerate().for_each(|(idx, l)| {
-            bop.add_line_string(l, idx + 1);
-        });
+        ls.0.iter().for_each(|l| bop.add_line_string(l, 0));
         bop.sweep()
     }
 }
 
+impl<'a, T: GeoNum> BopGeometry<'a, T> {
+    fn coords_count(&self) -> usize {
+        match self {
+            BopGeometry::Polygon(polygon) => polygon.coords_count(),
+            BopGeometry::MultiPolygon(multi_polygon) => multi_polygon.coords_count(),
+        }
+    }
+}
+
 mod op;
 use op::*;
 mod assembly;
diff --git a/geo/src/algorithm/bool_ops/tests.rs b/geo/src/algorithm/bool_ops/tests.rs
index 11373dcbd..33b04ffd9 100644
--- a/geo/src/algorithm/bool_ops/tests.rs
+++ b/geo/src/algorithm/bool_ops/tests.rs
@@ -1,4 +1,4 @@
-use crate::{LineString, MultiPolygon, Polygon};
+use crate::{LineString, MultiPolygon, Polygon, Relate};
 use log::{error, info};
 
 use std::{
@@ -116,6 +116,39 @@ fn test_complex_rects() -> Result<()> {
     }
     Ok(())
 }
+
+#[test]
+fn single_and_multi() {
+    let wkt1 = "POLYGON ((110 310, 220 310, 220 210, 110 210, 110 310))";
+    // multipolygon containing a single polygon
+    let wkt2 = "MULTIPOLYGON (((80 260, 90 260, 90 250, 80 250, 80 260)))";
+    // From JTS union op
+    let res = "MULTIPOLYGON (((110 310, 220 310, 220 210, 110 210, 110 310)), ((80 260, 90 260, 90 250, 80 250, 80 260)))";
+    let poly = crate::wkt!(POLYGON ((110 310, 220 310, 220 210, 110 210, 110 310)));
+    let mpoly = MultiPolygon::<f64>::try_from_wkt_str(wkt2).unwrap();
+    let respoly = MultiPolygon::<f64>::try_from_wkt_str(res).unwrap();
+    let union = mpoly.union(&poly);
+    let intersection_matrix = respoly.relate(&union);
+    // coords will be arranged differently, but we only care about topology
+    assert!(intersection_matrix.is_equal_topo());
+}
+
+#[test]
+fn multi_and_single() {
+    let wkt1 = "POLYGON ((110 310, 220 310, 220 210, 110 210, 110 310))";
+    // multipolygon containing a single polygon
+    let wkt2 = "MULTIPOLYGON (((80 260, 90 260, 90 250, 80 250, 80 260)))";
+    // From JTS union op
+    let res = "MULTIPOLYGON (((110 310, 220 310, 220 210, 110 210, 110 310)), ((80 260, 90 260, 90 250, 80 250, 80 260)))";
+    let poly = Polygon::<f64>::try_from_wkt_str(wkt1).unwrap();
+    let mpoly = MultiPolygon::<f64>::try_from_wkt_str(wkt2).unwrap();
+    let respoly = MultiPolygon::<f64>::try_from_wkt_str(res).unwrap();
+    let union = poly.union(&mpoly);
+    let intersection_matrix = respoly.relate(&union);
+    // coords will be arranged differently, but we only care about topology
+    assert!(intersection_matrix.is_equal_topo());
+}
+
 #[test]
 fn test_complex_rects1() -> Result<()> {
     let wkt1 = "MULTIPOLYGON(((-1 -2,-1.0000000000000002 2,-0.8823529411764707 2,-0.8823529411764706 -2,-1 -2)))";