diff --git a/Changelog.md b/Changelog.md
--- a/Changelog.md
+++ b/Changelog.md
@@ -1,5 +1,8 @@
 # Change log for esqueleto-postgis project
 
+## Version 2.1.0 
+* add st_unions
+
 ## Version 2.0.1 
 * drop haskell works hedghog dependency
 
diff --git a/esqueleto-postgis.cabal b/esqueleto-postgis.cabal
--- a/esqueleto-postgis.cabal
+++ b/esqueleto-postgis.cabal
@@ -5,7 +5,7 @@
 -- see: https://github.com/sol/hpack
 
 name:           esqueleto-postgis
-version:        2.0.1
+version:        2.1.0
 homepage:       https://github.com/jappeace/esqueleto-postgis#readme
 bug-reports:    https://github.com/jappeace/esqueleto-postgis/issues
 author:         Jappie Klooster
diff --git a/src/Database/Esqueleto/Postgis.hs b/src/Database/Esqueleto/Postgis.hs
--- a/src/Database/Esqueleto/Postgis.hs
+++ b/src/Database/Esqueleto/Postgis.hs
@@ -13,6 +13,7 @@
     st_contains,
     st_intersects,
     st_union,
+    st_unions,
 
     -- * points
     point,
@@ -46,7 +47,7 @@
 import Data.Text.Lazy.Builder (toLazyText)
 import Data.Text.Lazy.Builder qualified as Text
 import Database.Esqueleto.Experimental (SqlExpr, Value)
-import Database.Esqueleto.Internal.Internal (unsafeSqlFunction)
+import Database.Esqueleto.Internal.Internal (unsafeSqlFunction, unsafeSqlCastAs )
 import Database.Persist.Sql
 import Data.Base16.Types(assertBase16)
 import Data.ByteString(fromStrict)
@@ -71,9 +72,10 @@
 tshow :: (Show a) => a -> Text
 tshow = pack . show
 
--- | like 'GeospatialGeometry' but not partial, eg no empty geometries, also only works
---   in a single dimention, eg PostgisGeometry PointXY can't work with PostgisGeometry PointXYZ.
---   so PointXY indicates a 2 dimension space, and PointXYZ a three dimension space.
+-- | like 'GeospatialGeometry' but not partial, eg no empty geometries.
+--   Also can put an inveriant on dimensions if a function requires it.
+--   for example 'st_intersects' 'PostgisGeometry' 'PointXY' can't work with 'PostgisGeometry' 'PointXYZ'.
+--   PointXY indicates a 2 dimension space, and PointXYZ a three dimension space.
 data PostgisGeometry point
   = Point point
   | MultiPoint (NonEmpty point)
@@ -245,7 +247,10 @@
   SqlExpr (Value Bool)
 st_contains a b = unsafeSqlFunction "ST_CONTAINS" (a, b)
 
--- | allows union of geometries, eg group a bunch together, for example:
+
+-- | allows union of geometries, eg group a bunch together,
+--   https://postgis.net/docs/ST_Union.html
+--   for example:
 --
 -- @
 --  mCombined <- selectOne $ do
@@ -263,6 +268,16 @@
   SqlExpr (Value (PostgisGeometry a)) ->
   SqlExpr (Value (PostgisGeometry a))
 st_union a = unsafeSqlFunction "ST_union" a
+
+st_unions ::
+  SqlExpr (Value (PostgisGeometry a)) ->
+  SqlExpr (Value (PostgisGeometry a)) ->
+  SqlExpr (Value (PostgisGeometry a))
+st_unions a b =
+  -- casts to prevent
+  -- function st_union(unknown, unknown) is not unique", sqlErrorDetail = "", sqlErrorHint = "Could not choose a best candidate function. You might need to add explicit type casts.
+  -- TODO shouldn't we use sqlType here?
+  unsafeSqlFunction "ST_union" ((unsafeSqlCastAs "geometry" a), (unsafeSqlCastAs "geometry" b))
 
 -- | Returns true if two geometries intersect.
 --   Geometries intersect if they have any point in common.
diff --git a/test/Test.hs b/test/Test.hs
--- a/test/Test.hs
+++ b/test/Test.hs
@@ -218,6 +218,12 @@
                 grid <- from $ table @Grid
                 pure $ st_union $ grid ^. GridGeom
             unValue <$> result @?= (Just $ Polygon $ makeLinearRing (PointXY {_xyX = 0.0, _xyY = 2.0}) (PointXY {_xyX = 2.0, _xyY = 2.0}) (PointXY {_xyX = 4.0, _xyY = 2.0}) (Seq.fromList [PointXY {_xyX = 4.0, _xyY = 0.0}, PointXY {_xyX = 2.0, _xyY = 0.0}, PointXY {_xyX = 0.0, _xyY = 0.0}])),
+          testCase ("see if we can unions in PG and then get out some Haskell") $ do
+            result <- runDB $ do
+              selectOne $
+                pure $ st_unions (val (Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [(PointXY 2 0)])) $
+                        val $ Polygon $ makePolygon (PointXY 2 0) (PointXY 2 2) (PointXY 4 2) $ Seq.fromList [(PointXY 4 0)]
+            unValue <$> result @?= (Just $ Polygon $ makeLinearRing (PointXY {_xyX = 0.0, _xyY = 2.0}) (PointXY {_xyX = 2.0, _xyY = 2.0}) (PointXY {_xyX = 4.0, _xyY = 2.0}) (Seq.fromList [PointXY {_xyX = 4.0, _xyY = 0.0}, PointXY {_xyX = 2.0, _xyY = 0.0}, PointXY {_xyX = 0.0, _xyY = 0.0}])),
           testCase ("see if we can get just the units in the polygons") $ do
             result <- runDB $ do
               _ <-