diff --git a/Changelog.md b/Changelog.md
--- a/Changelog.md
+++ b/Changelog.md
@@ -1,5 +1,38 @@
 # Change log for esqueleto-postgis project
 
+## Version 4.1.0
++ Add 95 PostGIS function bindings covering the bulk of the PostGIS 3.6 geometry API.
++ Geometry constructors: `st_collect`, `st_makeenvelope`, `st_makeline`, `st_makepolygon_line`.
++ Geometry accessors: `st_boundary`, `st_coorddim`, `st_endpoint`, `st_exteriorring`,
+  `st_geometryn`, `st_geometrytype`, `st_interiorringn`, `st_iscollection`, `st_isempty`,
+  `st_ispolygonccw`, `st_ispolygoncw`, `st_isring`, `st_m`, `st_ndims`, `st_nrings`,
+  `st_numinteriorrings`, `st_numpoints`, `st_pointn`, `st_startpoint`, `st_z`.
++ Geometry editors: `st_addpoint`, `st_collectionextract`, `st_flipcoordinates`,
+  `st_force2d`, `st_force3d`, `st_force4d`, `st_forcecollection`,
+  `st_forcepolygonccw`, `st_forcepolygoncw`, `st_multi`, `st_normalize`,
+  `st_reverse`, `st_segmentize`, `st_setpoint`, `st_snaptogrid`, `st_snap`.
++ Geometry validation: `st_makevalid`, `st_isvalidreason`.
++ SRS: `st_setsrid`.
++ Geometry output: `st_astext`, `st_asgeojson`, `st_asewkt`, `st_geohash`.
++ Spatial relationships: `st_3dintersects`, `st_relate`, `st_orderingequals`,
+  `st_dfullywithin`, `st_pointinsidecircle`.
++ Measurements: `st_angle`, `st_closestpoint`, `st_3dclosestpoint`, `st_3ddistance`,
+  `st_distancesphere`, `st_frechetdistance`, `st_hausdorffdistance`, `st_length2d`,
+  `st_3dlength`, `st_longestline`, `st_3dlongestline`, `st_3dmaxdistance`,
+  `st_minimumclearance`, `st_shortestline`, `st_3dshortestline`.
++ Overlay functions: `st_symdifference`, `st_unaryunion`, `st_split`, `st_node`.
++ Geometry processing: `st_buildarea`, `st_chaikinsmoothing`, `st_concavehull`,
+  `st_delaunaytriangles`, `st_generatepoints`, `st_geometricmedian`, `st_linemerge`,
+  `st_minimumboundingcircle`, `st_offsetcurve`, `st_reduceprecision`, `st_sharedpaths`,
+  `st_simplify`, `st_simplifypreservetopology`, `st_voronoilines`, `st_voronoipolygons`.
++ Affine transformations: `st_translate`, `st_scale`, `st_rotate`, `st_rotatex`, `st_rotatez`.
++ Bounding box: `st_expand`.
++ Linear referencing: `st_lineinterpolatepoint`, `st_linelocatepoint`, `st_linesubstring`.
++ Integration test for every new function.
+
+## Version 4.0.3
++ add comments on x and y, z and m.
+
 ## Version 4.0.2 
 + fix bug where st_point and friends did not respect spatialType.
   solved it by always casting them to whatever the user specified.
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:        4.0.2
+version:        4.1.0
 homepage:       https://github.com/jappeace/esqueleto-postgis#readme
 bug-reports:    https://github.com/jappeace/esqueleto-postgis/issues
 author:         Jappie Klooster
@@ -79,6 +79,12 @@
       Database.Esqueleto.Postgis.Ewkb
       Database.Esqueleto.Postgis.Wkb
   other-modules:
+      Database.Esqueleto.Postgis.Geometry
+      Database.Esqueleto.Postgis.Accessor
+      Database.Esqueleto.Postgis.Editor
+      Database.Esqueleto.Postgis.Spatial
+      Database.Esqueleto.Postgis.Measurement
+      Database.Esqueleto.Postgis.Processing
       Database.Esqueleto.Postgis.Wkb.Endian
       Database.Esqueleto.Postgis.Wkb.Point
       Database.Esqueleto.Postgis.Wkb.Geometry
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
@@ -15,8 +15,7 @@
 --   You don't have to use Geography if you're only interested in
 --   topological functions such as 'st_intersects' and 'st_union',
 --   these are indifferent to space distortions,
---   see <https://www.gdmc.nl/publications/2013/3D_Geo-DBMS_Topological_Operators.pdf>
---   the related work section.
+--   see <https://jappie.me/announcement-esqueleto-postgis-v4.html>
 --
 --   if you can't use a function for example when you're using 'Geography'.
 --   there is the option to 'st_transform_geography'.
@@ -26,44 +25,91 @@
     SpatialType(..),
     getPoints,
 
-    -- * functions
-    st_contains,
-    st_intersects,
-    st_union,
-    st_unions,
-    st_dwithin,
-    st_distance         ,
+    -- * Spatial relationships
+    module Database.Esqueleto.Postgis.Spatial,
 
-    -- * points
+    -- * Measurement functions
+    module Database.Esqueleto.Postgis.Measurement,
+
+    -- * Geometry accessors
+    module Database.Esqueleto.Postgis.Accessor,
+
+    -- * Geometry constructors
+    st_collect,
+    st_makeenvelope,
+    st_makeline,
+    st_makepolygon_line,
+
+    -- * Geometry editors
+    module Database.Esqueleto.Postgis.Editor,
+
+    -- * Geometry validation
+    st_makevalid,
+    st_isvalidreason,
+
+    -- * SRS functions
+    st_setsrid,
+
+    -- * Geometry output
+    st_astext,
+    st_asgeojson,
+    st_asewkt,
+    st_geohash,
+
+    -- * Geometry processing
+    module Database.Esqueleto.Postgis.Processing,
+
+    -- * Affine transformations
+    st_translate,
+    st_scale,
+    st_rotate,
+    st_rotatex,
+    st_rotatez,
+
+    -- * Bounding box
+    st_expand,
+
+    -- * Linear referencing
+    st_lineinterpolatepoint,
+    st_linelocatepoint,
+    st_linesubstring,
+
+    -- * Points
     point,
     point_v,
     st_point,
     st_point_xyz,
     st_point_xyzm,
 
-    -- * transform
+    -- * Transform
     st_transform_geography,
     st_transform_geometry,
-    -- ** srid
+    -- ** SRID
     SRID,
     wgs84,
-    mercator ,
+    mercator,
     britishNationalGrid,
-    SridUnit    (..),
-
+    SridUnit(..),
 
-    -- * other
+    -- * Other
     makePolygon,
     PostgisGeometry,
     HasPgType,
 
-    -- * re-exports
+    -- * Re-exports
     PointXY(..),
     PointXYZ(..),
     PointXYZM(..),
   )
 where
 
+import Database.Esqueleto.Postgis.Geometry (Postgis(..), SpatialType(..), HasPgType(..), PostgisGeometry)
+import Database.Esqueleto.Postgis.Spatial
+import Database.Esqueleto.Postgis.Measurement
+import Database.Esqueleto.Postgis.Accessor
+import Database.Esqueleto.Postgis.Editor
+import Database.Esqueleto.Postgis.Processing
+
 import Database.Esqueleto.Experimental(val)
 import Data.Proxy
 import Data.Bifunctor (first)
@@ -111,42 +157,6 @@
 tshow :: (Show a) => a -> Text
 tshow = pack . show
 
-
--- | Guarantees we don't accidently mix curved space with flat space.
---   Postgis will catch this too, this just put's it in the type system.
-data SpatialType = Geometry -- ^ assume a flat space.
-                 | Geography -- ^ assume curvature of the earth.
-
--- | technical typeclass to bind a spatial type to a string value.
---   because we represent the constructors as a datakind, we need
---   this to go back to a value.
-class HasPgType (spatialType :: SpatialType) where
-  pgType :: Proxy spatialType -> Text
-
-instance HasPgType 'Geometry where
-  pgType _ = "geometry"
-
-instance HasPgType 'Geography where
-  pgType _ = "geography"
-
--- | backwards compatibility, initial version only dealt in geometry
-type PostgisGeometry = Postgis 'Geometry
-
-
--- | 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 Postgis (spatialType :: SpatialType) point
-  = Point point
-  | MultiPoint (NonEmpty point)
-  | Line (LineString point)
-  | Multiline (NonEmpty (LineString point))
-  | Polygon (LinearRing point)
-  | MultiPolygon (NonEmpty (LinearRing point))
-  | Collection (NonEmpty (PostgisGeometry point))
-  deriving (Show, Functor, Eq)
-
 data GeomErrors
   = MismatchingDimensionsXYZ PointXYZ
   | MismatchingDimensionsXYZM PointXYZM
@@ -290,29 +300,6 @@
 instance HasPgType spatialType => PersistFieldSql (Postgis spatialType PointXYZM) where
   sqlType _ = SqlOther $ pgType $ Proxy @spatialType
 
-
--- | Returns TRUE if geometry A contains geometry B.
---   https://postgis.net/docs/ST_Contains.html
-st_contains ::
-  -- | geom a
-  SqlExpr (Value (Postgis 'Geometry a)) ->
-  -- | geom b
-  SqlExpr (Value (Postgis 'Geometry a)) ->
-  SqlExpr (Value Bool)
-st_contains a b = unsafeSqlFunction "ST_CONTAINS" (a, b)
-
--- | Returns true if the geometries are within a given distance
---   https://postgis.net/docs/ST_DWithin.html
-st_dwithin ::
-  -- | geometry g1
-  SqlExpr (Value (Postgis spatialType a)) ->
-  -- | geometry g2
-  SqlExpr (Value (Postgis spatialType a)) ->
-  -- | distance of srid
-  SqlExpr (Value Double) ->
-  SqlExpr (Value Bool)
-st_dwithin a b c = unsafeSqlFunction "ST_DWithin" (a, b, c)
-
 -- | SRID
 -- you can find your local like this: https://blog.rustprooflabs.com/2020/11/postgis-find-local-srid
 -- geography appears to use 'wgs84'. So I hardcoded the use going from geom to geography as that.
@@ -373,74 +360,229 @@
   where
     castAs = pgType $ Proxy @two
 
--- | allows union of geometries, eg group a bunch together,
---   https://postgis.net/docs/ST_Union.html
---   for example:
---
--- @
---  mCombined <- selectOne $ do
---    grid <- from $ table @Grid
---    pure $ st_union $ grid ^. GridGeom
---
---
---  select $  do
---    unit <- from $ table @Unit
---    forM_ mCombined $ \combined ->
---      where_ $ (unit ^. UnitGeom) `st_intersects` (val $ unValue combined)
---    pure unit
--- @
-st_union ::
+-- ---------------------------------------------------------------------------
+-- Geometry Constructors
+-- ---------------------------------------------------------------------------
+
+-- | Aggregate function that collects geometries into a geometry collection.
+--   https://postgis.net/docs/ST_Collect.html
+st_collect ::
   SqlExpr (Value (Postgis 'Geometry a)) ->
   SqlExpr (Value (Postgis 'Geometry a))
-st_union a = unsafeSqlFunction "ST_union" a
+st_collect a = unsafeSqlFunction "ST_Collect" a
 
-st_unions ::
+-- | Creates a rectangular polygon from minimum and maximum coordinates.
+--   https://postgis.net/docs/ST_MakeEnvelope.html
+st_makeenvelope ::
+  SqlExpr (Value Double) -> -- ^ xmin
+  SqlExpr (Value Double) -> -- ^ ymin
+  SqlExpr (Value Double) -> -- ^ xmax
+  SqlExpr (Value Double) -> -- ^ ymax
+  SqlExpr (Value Int) ->    -- ^ srid
+  SqlExpr (Value (Postgis 'Geometry a))
+st_makeenvelope xmin' ymin' xmax' ymax' srid' = unsafeSqlFunction "ST_MakeEnvelope" (xmin', ymin', xmax', ymax', srid')
+
+-- | Creates a linestring from two point geometries.
+--   https://postgis.net/docs/ST_MakeLine.html
+st_makeline ::
   SqlExpr (Value (Postgis 'Geometry a)) ->
   SqlExpr (Value (Postgis 'Geometry a)) ->
   SqlExpr (Value (Postgis 'Geometry 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.
-  unsafeSqlFunction "ST_union" ((unsafeSqlCastAs casted a), (unsafeSqlCastAs casted b))
-  where
-    casted = (pgType $ Proxy @'Geometry)
+st_makeline a b = unsafeSqlFunction "ST_MakeLine" (a, b)
 
--- | calculate the distance between two points
---   https://postgis.net/docs/ST_Distance.html
-st_distance ::
-  SqlExpr (Value (Postgis spatialType a)) ->
-  SqlExpr (Value (Postgis spatialType a)) ->
+-- | Creates a polygon from a closed linestring shell.
+--   Named st_makepolygon_line to avoid collision with 'makePolygon'.
+--   https://postgis.net/docs/ST_MakePolygon.html
+st_makepolygon_line ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_makepolygon_line a = unsafeSqlFunction "ST_MakePolygon" a
+
+-- ---------------------------------------------------------------------------
+-- Geometry Validation
+-- ---------------------------------------------------------------------------
+
+-- | Attempts to make an invalid geometry valid without losing vertices.
+--   https://postgis.net/docs/ST_MakeValid.html
+st_makevalid ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_makevalid a = unsafeSqlFunction "ST_MakeValid" a
+
+-- | Returns text describing why a geometry is invalid, or "Valid Geometry".
+--   https://postgis.net/docs/ST_IsValidReason.html
+st_isvalidreason ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Text)
+st_isvalidreason a = unsafeSqlFunction "ST_IsValidReason" a
+
+-- ---------------------------------------------------------------------------
+-- SRS Functions
+-- ---------------------------------------------------------------------------
+
+-- | Sets the SRID on a geometry to a particular integer value.
+--   https://postgis.net/docs/ST_SetSRID.html
+st_setsrid ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_setsrid a srid' = unsafeSqlFunction "ST_SetSRID" (a, srid')
+
+-- ---------------------------------------------------------------------------
+-- Geometry Output
+-- ---------------------------------------------------------------------------
+
+-- | Returns the Well-Known Text (WKT) representation of the geometry.
+--   https://postgis.net/docs/ST_AsText.html
+st_astext ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Text)
+st_astext a = unsafeSqlFunction "ST_AsText" a
+
+-- | Returns the GeoJSON representation of the geometry.
+--   https://postgis.net/docs/ST_AsGeoJSON.html
+st_asgeojson ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Text)
+st_asgeojson a = unsafeSqlFunction "ST_AsGeoJSON" a
+
+-- | Returns the Extended Well-Known Text (EWKT) representation of the geometry.
+--   https://postgis.net/docs/ST_AsEWKT.html
+st_asewkt ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Text)
+st_asewkt a = unsafeSqlFunction "ST_AsEWKT" a
+
+-- | Returns a GeoHash representation of the geometry.
+--   https://postgis.net/docs/ST_GeoHash.html
+st_geohash ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Text)
+st_geohash a = unsafeSqlFunction "ST_GeoHash" a
+
+-- ---------------------------------------------------------------------------
+-- Affine Transformations
+-- ---------------------------------------------------------------------------
+
+-- | Translates a geometry by given X and Y offsets.
+--   https://postgis.net/docs/ST_Translate.html
+st_translate ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) -> -- ^ deltaX
+  SqlExpr (Value Double) -> -- ^ deltaY
+  SqlExpr (Value (Postgis 'Geometry a))
+st_translate a dx dy = unsafeSqlFunction "ST_Translate" (a, dx, dy)
+
+-- | Scales a geometry by given X and Y factors.
+--   https://postgis.net/docs/ST_Scale.html
+st_scale ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) -> -- ^ scaleX
+  SqlExpr (Value Double) -> -- ^ scaleY
+  SqlExpr (Value (Postgis 'Geometry a))
+st_scale a sx sy = unsafeSqlFunction "ST_Scale" (a, sx, sy)
+
+-- | Rotates a geometry around the origin by an angle in radians.
+--   https://postgis.net/docs/ST_Rotate.html
+st_rotate ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) -> -- ^ angle in radians
+  SqlExpr (Value (Postgis 'Geometry a))
+st_rotate a angle = unsafeSqlFunction "ST_Rotate" (a, angle)
+
+-- | Rotates a geometry around the X axis by an angle in radians.
+--   https://postgis.net/docs/ST_RotateX.html
+st_rotatex ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) -> -- ^ angle in radians
+  SqlExpr (Value (Postgis 'Geometry a))
+st_rotatex a angle = unsafeSqlFunction "ST_RotateX" (a, angle)
+
+-- | Rotates a geometry around the Z axis by an angle in radians.
+--   https://postgis.net/docs/ST_RotateZ.html
+st_rotatez ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) -> -- ^ angle in radians
+  SqlExpr (Value (Postgis 'Geometry a))
+st_rotatez a angle = unsafeSqlFunction "ST_RotateZ" (a, angle)
+
+-- ---------------------------------------------------------------------------
+-- Bounding Box
+-- ---------------------------------------------------------------------------
+
+-- | Returns a bounding box expanded in all directions by a given distance.
+--   https://postgis.net/docs/ST_Expand.html
+st_expand ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_expand a d = unsafeSqlFunction "ST_Expand" (a, d)
+
+-- ---------------------------------------------------------------------------
+-- Linear Referencing
+-- ---------------------------------------------------------------------------
+
+-- | Returns a point interpolated along a line at a fractional distance (0.0 to 1.0).
+--   https://postgis.net/docs/ST_LineInterpolatePoint.html
+st_lineinterpolatepoint ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) -> -- ^ fraction (0.0 to 1.0)
+  SqlExpr (Value (Postgis 'Geometry a))
+st_lineinterpolatepoint a f = unsafeSqlFunction "ST_LineInterpolatePoint" (a, f)
+
+-- | Returns a float between 0 and 1 representing the location of the closest point on a line to a given point.
+--   https://postgis.net/docs/ST_LineLocatePoint.html
+st_linelocatepoint ::
+  SqlExpr (Value (Postgis 'Geometry a)) -> -- ^ line
+  SqlExpr (Value (Postgis 'Geometry a)) -> -- ^ point
   SqlExpr (Value Double)
-st_distance a b =
-  unsafeSqlFunction "ST_distance" (a, b)
+st_linelocatepoint a b = unsafeSqlFunction "ST_LineLocatePoint" (a, b)
 
--- | Returns true if two geometries intersect.
---   Geometries intersect if they have any point in common.
---   https://postgis.net/docs/ST_Intersects.html
-st_intersects ::
-  SqlExpr (Value (Postgis spatialType a)) -> -- ^ geomA or geogA
-  SqlExpr (Value (Postgis spatialType a)) -> -- ^ geomB or geogB
-  SqlExpr (Value Bool)
-st_intersects a b = unsafeSqlFunction "ST_Intersects" (a, b)
+-- | Returns the portion of a line between two fractional locations (0.0 to 1.0).
+--   https://postgis.net/docs/ST_LineSubstring.html
+st_linesubstring ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) -> -- ^ start fraction
+  SqlExpr (Value Double) -> -- ^ end fraction
+  SqlExpr (Value (Postgis 'Geometry a))
+st_linesubstring a s e = unsafeSqlFunction "ST_LineSubstring" (a, s, e)
 
-point :: Double -> Double -> (Postgis spatialType PointXY)
+point ::
+  Double -> -- ^ x or longitude
+  Double -> -- ^ y or latitude
+  (Postgis spatialType PointXY)
 point x y = Point (PointXY {_xyX = x, _xyY = y})
 
-point_v :: HasPgType spatialType => Double -> Double -> SqlExpr (Value (Postgis spatialType PointXY))
+point_v :: HasPgType spatialType =>
+  Double -> -- ^ x or longitude
+  Double -> -- ^ y or latitude
+  SqlExpr (Value (Postgis spatialType PointXY))
 point_v = fmap val . point
 
-st_point ::  forall spatialType . HasPgType spatialType => SqlExpr (Value Double) -> SqlExpr (Value Double) -> SqlExpr (Value (Postgis spatialType PointXY))
+st_point ::  forall spatialType . HasPgType spatialType =>
+  SqlExpr (Value Double) -> -- ^ x or longitude
+  SqlExpr (Value Double) -> -- ^ y or latitude
+  SqlExpr (Value (Postgis spatialType PointXY))
 st_point a b = unsafeSqlCastAs castAs $ unsafeSqlFunction "ST_POINT" (a, b)
   where
     castAs = pgType $ Proxy @spatialType
 
 
-st_point_xyz ::  forall spatialType . HasPgType spatialType =>  SqlExpr (Value Double) -> SqlExpr (Value Double) -> SqlExpr (Value Double) -> SqlExpr (Value (Postgis spatialType PointXYZ))
+st_point_xyz ::  forall spatialType . HasPgType spatialType =>
+  SqlExpr (Value Double) -> -- ^ x or longitude
+  SqlExpr (Value Double) -> -- ^ y or latitude
+  SqlExpr (Value Double) -> -- ^ z elevation/altitude
+  SqlExpr (Value (Postgis spatialType PointXYZ))
 st_point_xyz a b c = unsafeSqlCastAs castAs $ unsafeSqlFunction "ST_POINT" (a, b, c)
   where
     castAs = pgType $ Proxy @spatialType
 
-st_point_xyzm :: forall spatialType . HasPgType spatialType =>  SqlExpr (Value Double) -> SqlExpr (Value Double) -> SqlExpr (Value Double) -> SqlExpr (Value Double) -> SqlExpr (Value (Postgis spatialType PointXYZM))
+st_point_xyzm :: forall spatialType . HasPgType spatialType =>
+  SqlExpr (Value Double) -> -- ^ x or longitude
+  SqlExpr (Value Double) -> -- ^ y or latitude
+  SqlExpr (Value Double) -> -- ^ z elevation/altitude
+  SqlExpr (Value Double) -> -- ^ m measure, user defined dimension
+  SqlExpr (Value (Postgis spatialType PointXYZM))
 st_point_xyzm a b c m = unsafeSqlCastAs castAs $ unsafeSqlFunction "ST_POINT" (a, b, c, m)
   where
     castAs = pgType $ Proxy @spatialType
diff --git a/src/Database/Esqueleto/Postgis/Accessor.hs b/src/Database/Esqueleto/Postgis/Accessor.hs
new file mode 100644
--- /dev/null
+++ b/src/Database/Esqueleto/Postgis/Accessor.hs
@@ -0,0 +1,253 @@
+{-# LANGUAGE DataKinds #-}
+
+module Database.Esqueleto.Postgis.Accessor
+  ( st_x
+  , st_y
+  , st_npoints
+  , st_numgeometries
+  , st_dimension
+  , st_issimple
+  , st_isclosed
+  , st_isvalid
+  , st_srid
+  , st_boundary
+  , st_coorddim
+  , st_endpoint
+  , st_exteriorring
+  , st_geometryn
+  , st_geometrytype
+  , st_interiorringn
+  , st_iscollection
+  , st_isempty
+  , st_ispolygonccw
+  , st_ispolygoncw
+  , st_isring
+  , st_m
+  , st_ndims
+  , st_nrings
+  , st_numinteriorrings
+  , st_numpoints
+  , st_pointn
+  , st_startpoint
+  , st_z
+  ) where
+
+import Database.Esqueleto.Postgis.Geometry (Postgis, SpatialType(..))
+import Database.Esqueleto.Experimental (SqlExpr, Value)
+import Database.Esqueleto.Internal.Internal (unsafeSqlFunction)
+import Data.Text (Text)
+
+-- | Returns the X coordinate of a point geometry.
+--   Only works on points; other geometry types will error.
+--   https://postgis.net/docs/ST_X.html
+st_x ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_x a = unsafeSqlFunction "ST_X" a
+
+-- | Returns the Y coordinate of a point geometry.
+--   Only works on points; other geometry types will error.
+--   https://postgis.net/docs/ST_Y.html
+st_y ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_y a = unsafeSqlFunction "ST_Y" a
+
+-- | Returns the number of points (vertices) in a geometry.
+--   Works on any geometry type.
+--   https://postgis.net/docs/ST_NPoints.html
+st_npoints ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int)
+st_npoints a = unsafeSqlFunction "ST_NPoints" a
+
+-- | Returns the number of sub-geometries in a geometry collection or multi-type.
+--   Returns 1 for single geometries.
+--   https://postgis.net/docs/ST_NumGeometries.html
+st_numgeometries ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int)
+st_numgeometries a = unsafeSqlFunction "ST_NumGeometries" a
+
+-- | Returns the topological dimension of a geometry.
+--   0 for points, 1 for lines, 2 for polygons.
+--   https://postgis.net/docs/ST_Dimension.html
+st_dimension ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int)
+st_dimension a = unsafeSqlFunction "ST_Dimension" a
+
+-- | Returns TRUE if the geometry has no self-intersections.
+--   Points and properly-formed polygons are always simple.
+--   https://postgis.net/docs/ST_IsSimple.html
+st_issimple ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_issimple a = unsafeSqlFunction "ST_IsSimple" a
+
+-- | Returns TRUE if the linestring's start and end points are coincident.
+--   For polyhedral surfaces, reports if the surface is areal (open) or volumetric (closed).
+--   https://postgis.net/docs/ST_IsClosed.html
+st_isclosed ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_isclosed a = unsafeSqlFunction "ST_IsClosed" a
+
+-- | Returns TRUE if the geometry is well-formed and valid per the OGC rules.
+--   Points and lines are always valid; polygons need correct ring structure.
+--   https://postgis.net/docs/ST_IsValid.html
+st_isvalid ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_isvalid a = unsafeSqlFunction "ST_IsValid" a
+
+-- | Returns the spatial reference identifier (SRID) of the geometry.
+--   0 means no SRID is set.
+--   https://postgis.net/docs/ST_SRID.html
+st_srid ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int)
+st_srid a = unsafeSqlFunction "ST_SRID" a
+
+-- | Returns the boundary of a geometry.
+--   https://postgis.net/docs/ST_Boundary.html
+st_boundary ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_boundary a = unsafeSqlFunction "ST_Boundary" a
+
+-- | Returns the coordinate dimension of a geometry (2, 3, or 4).
+--   https://postgis.net/docs/ST_CoordDim.html
+st_coorddim ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int)
+st_coorddim a = unsafeSqlFunction "ST_CoordDim" a
+
+-- | Returns the last point of a linestring.
+--   https://postgis.net/docs/ST_EndPoint.html
+st_endpoint ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_endpoint a = unsafeSqlFunction "ST_EndPoint" a
+
+-- | Returns the exterior ring of a polygon geometry.
+--   https://postgis.net/docs/ST_ExteriorRing.html
+st_exteriorring ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_exteriorring a = unsafeSqlFunction "ST_ExteriorRing" a
+
+-- | Returns the Nth sub-geometry of a geometry collection (1-indexed).
+--   https://postgis.net/docs/ST_GeometryN.html
+st_geometryn ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_geometryn a n = unsafeSqlFunction "ST_GeometryN" (a, n)
+
+-- | Returns the geometry type as a string (e.g. "ST_Polygon").
+--   https://postgis.net/docs/ST_GeometryType.html
+st_geometrytype ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Text)
+st_geometrytype a = unsafeSqlFunction "ST_GeometryType" a
+
+-- | Returns the Nth interior ring of a polygon (1-indexed).
+--   https://postgis.net/docs/ST_InteriorRingN.html
+st_interiorringn ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_interiorringn a n = unsafeSqlFunction "ST_InteriorRingN" (a, n)
+
+-- | Returns TRUE if the geometry is a geometry collection type.
+--   https://postgis.net/docs/ST_IsCollection.html
+st_iscollection ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_iscollection a = unsafeSqlFunction "ST_IsCollection" a
+
+-- | Returns TRUE if the geometry is an empty geometry.
+--   https://postgis.net/docs/ST_IsEmpty.html
+st_isempty ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_isempty a = unsafeSqlFunction "ST_IsEmpty" a
+
+-- | Returns TRUE if all polygon rings are oriented counter-clockwise.
+--   https://postgis.net/docs/ST_IsPolygonCCW.html
+st_ispolygonccw ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_ispolygonccw a = unsafeSqlFunction "ST_IsPolygonCCW" a
+
+-- | Returns TRUE if all polygon rings are oriented clockwise.
+--   https://postgis.net/docs/ST_IsPolygonCW.html
+st_ispolygoncw ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_ispolygoncw a = unsafeSqlFunction "ST_IsPolygonCW" a
+
+-- | Returns TRUE if the linestring is closed and simple (a ring).
+--   https://postgis.net/docs/ST_IsRing.html
+st_isring ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_isring a = unsafeSqlFunction "ST_IsRing" a
+
+-- | Returns the M coordinate of a point.
+--   https://postgis.net/docs/ST_M.html
+st_m ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_m a = unsafeSqlFunction "ST_M" a
+
+-- | Returns the number of dimensions of a geometry's coordinates (2, 3, or 4).
+--   https://postgis.net/docs/ST_NDims.html
+st_ndims ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int)
+st_ndims a = unsafeSqlFunction "ST_NDims" a
+
+-- | Returns the number of rings in a polygon (exterior + interior).
+--   https://postgis.net/docs/ST_NRings.html
+st_nrings ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int)
+st_nrings a = unsafeSqlFunction "ST_NRings" a
+
+-- | Returns the number of interior rings of a polygon.
+--   https://postgis.net/docs/ST_NumInteriorRings.html
+st_numinteriorrings ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int)
+st_numinteriorrings a = unsafeSqlFunction "ST_NumInteriorRings" a
+
+-- | Returns the number of points in a linestring.
+--   https://postgis.net/docs/ST_NumPoints.html
+st_numpoints ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int)
+st_numpoints a = unsafeSqlFunction "ST_NumPoints" a
+
+-- | Returns the Nth point in a linestring (1-indexed).
+--   https://postgis.net/docs/ST_PointN.html
+st_pointn ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_pointn a n = unsafeSqlFunction "ST_PointN" (a, n)
+
+-- | Returns the first point of a linestring.
+--   https://postgis.net/docs/ST_StartPoint.html
+st_startpoint ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_startpoint a = unsafeSqlFunction "ST_StartPoint" a
+
+-- | Returns the Z coordinate of a point.
+--   https://postgis.net/docs/ST_Z.html
+st_z ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_z a = unsafeSqlFunction "ST_Z" a
diff --git a/src/Database/Esqueleto/Postgis/Editor.hs b/src/Database/Esqueleto/Postgis/Editor.hs
new file mode 100644
--- /dev/null
+++ b/src/Database/Esqueleto/Postgis/Editor.hs
@@ -0,0 +1,146 @@
+{-# LANGUAGE DataKinds #-}
+
+module Database.Esqueleto.Postgis.Editor
+  ( st_addpoint
+  , st_collectionextract
+  , st_flipcoordinates
+  , st_force2d
+  , st_force3d
+  , st_force4d
+  , st_forcecollection
+  , st_forcepolygonccw
+  , st_forcepolygoncw
+  , st_multi
+  , st_normalize
+  , st_reverse
+  , st_segmentize
+  , st_setpoint
+  , st_snaptogrid
+  , st_snap
+  ) where
+
+import Database.Esqueleto.Postgis.Geometry (Postgis, SpatialType(..))
+import Database.Esqueleto.Experimental (SqlExpr, Value)
+import Database.Esqueleto.Internal.Internal (unsafeSqlFunction)
+
+-- | Adds a point to a linestring at a given position (0-indexed).
+--   https://postgis.net/docs/ST_AddPoint.html
+st_addpoint ::
+  SqlExpr (Value (Postgis 'Geometry a)) -> -- ^ linestring
+  SqlExpr (Value (Postgis 'Geometry a)) -> -- ^ point
+  SqlExpr (Value Int) ->                   -- ^ position
+  SqlExpr (Value (Postgis 'Geometry a))
+st_addpoint a b n = unsafeSqlFunction "ST_AddPoint" (a, b, n)
+
+-- | Extracts sub-geometries of a specified type from a collection.
+--   Type: 1=POINT, 2=LINESTRING, 3=POLYGON.
+--   https://postgis.net/docs/ST_CollectionExtract.html
+st_collectionextract ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_collectionextract a t = unsafeSqlFunction "ST_CollectionExtract" (a, t)
+
+-- | Swaps X and Y coordinates.
+--   https://postgis.net/docs/ST_FlipCoordinates.html
+st_flipcoordinates ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_flipcoordinates a = unsafeSqlFunction "ST_FlipCoordinates" a
+
+-- | Forces the geometry into 2D mode.
+--   https://postgis.net/docs/ST_Force2D.html
+st_force2d ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_force2d a = unsafeSqlFunction "ST_Force2D" a
+
+-- | Forces the geometry into 3D (XYZ) mode.
+--   https://postgis.net/docs/ST_Force3D.html
+st_force3d ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_force3d a = unsafeSqlFunction "ST_Force3D" a
+
+-- | Forces the geometry into 4D (XYZM) mode.
+--   https://postgis.net/docs/ST_Force4D.html
+st_force4d ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_force4d a = unsafeSqlFunction "ST_Force4D" a
+
+-- | Converts the geometry into a geometry collection.
+--   https://postgis.net/docs/ST_ForceCollection.html
+st_forcecollection ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_forcecollection a = unsafeSqlFunction "ST_ForceCollection" a
+
+-- | Forces polygon rings to counter-clockwise orientation.
+--   https://postgis.net/docs/ST_ForcePolygonCCW.html
+st_forcepolygonccw ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_forcepolygonccw a = unsafeSqlFunction "ST_ForcePolygonCCW" a
+
+-- | Forces polygon rings to clockwise orientation.
+--   https://postgis.net/docs/ST_ForcePolygonCW.html
+st_forcepolygoncw ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_forcepolygoncw a = unsafeSqlFunction "ST_ForcePolygonCW" a
+
+-- | Returns the geometry as a multi-type (e.g. POINT -> MULTIPOINT).
+--   https://postgis.net/docs/ST_Multi.html
+st_multi ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_multi a = unsafeSqlFunction "ST_Multi" a
+
+-- | Returns the geometry in its canonical (normalized) form.
+--   https://postgis.net/docs/ST_Normalize.html
+st_normalize ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_normalize a = unsafeSqlFunction "ST_Normalize" a
+
+-- | Returns the geometry with vertex order reversed.
+--   https://postgis.net/docs/ST_Reverse.html
+st_reverse ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_reverse a = unsafeSqlFunction "ST_Reverse" a
+
+-- | Segmentizes a geometry so no segment is longer than the given distance.
+--   https://postgis.net/docs/ST_Segmentize.html
+st_segmentize ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_segmentize a d = unsafeSqlFunction "ST_Segmentize" (a, d)
+
+-- | Replaces a point of a linestring at a given 0-based index.
+--   https://postgis.net/docs/ST_SetPoint.html
+st_setpoint ::
+  SqlExpr (Value (Postgis 'Geometry a)) -> -- ^ linestring
+  SqlExpr (Value Int) ->                   -- ^ 0-based index
+  SqlExpr (Value (Postgis 'Geometry a)) -> -- ^ new point
+  SqlExpr (Value (Postgis 'Geometry a))
+st_setpoint a n p = unsafeSqlFunction "ST_SetPoint" (a, n, p)
+
+-- | Snaps all points of a geometry to a regular grid of the given size.
+--   https://postgis.net/docs/ST_SnapToGrid.html
+st_snaptogrid ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_snaptogrid a d = unsafeSqlFunction "ST_SnapToGrid" (a, d)
+
+-- | Snaps vertices of one geometry to vertices and edges of another within a tolerance.
+--   https://postgis.net/docs/ST_Snap.html
+st_snap ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_snap a b d = unsafeSqlFunction "ST_Snap" (a, b, d)
diff --git a/src/Database/Esqueleto/Postgis/Geometry.hs b/src/Database/Esqueleto/Postgis/Geometry.hs
new file mode 100644
--- /dev/null
+++ b/src/Database/Esqueleto/Postgis/Geometry.hs
@@ -0,0 +1,48 @@
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE KindSignatures #-}
+
+module Database.Esqueleto.Postgis.Geometry
+  ( SpatialType(..)
+  , HasPgType(..)
+  , Postgis(..)
+  , PostgisGeometry
+  ) where
+
+import Data.List.NonEmpty (NonEmpty)
+import Data.LineString (LineString)
+import Data.LinearRing (LinearRing)
+import Data.Text (Text)
+
+-- | Guarantees we don't accidently mix curved space with flat space.
+--   Postgis will catch this too, this just put's it in the type system.
+data SpatialType = Geometry -- ^ assume a flat space.
+                 | Geography -- ^ assume curvature of the earth.
+
+-- | technical typeclass to bind a spatial type to a string value.
+--   because we represent the constructors as a datakind, we need
+--   this to go back to a value.
+class HasPgType (spatialType :: SpatialType) where
+  pgType :: proxy spatialType -> Text
+
+instance HasPgType 'Geometry where
+  pgType _ = "geometry"
+
+instance HasPgType 'Geography where
+  pgType _ = "geography"
+
+-- | backwards compatibility, initial version only dealt in geometry
+type PostgisGeometry = Postgis 'Geometry
+
+-- | 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 Postgis (spatialType :: SpatialType) point
+  = Point point
+  | MultiPoint (NonEmpty point)
+  | Line (LineString point)
+  | Multiline (NonEmpty (LineString point))
+  | Polygon (LinearRing point)
+  | MultiPolygon (NonEmpty (LinearRing point))
+  | Collection (NonEmpty (PostgisGeometry point))
+  deriving (Show, Functor, Eq)
diff --git a/src/Database/Esqueleto/Postgis/Measurement.hs b/src/Database/Esqueleto/Postgis/Measurement.hs
new file mode 100644
--- /dev/null
+++ b/src/Database/Esqueleto/Postgis/Measurement.hs
@@ -0,0 +1,213 @@
+{-# LANGUAGE DataKinds #-}
+
+module Database.Esqueleto.Postgis.Measurement
+  ( st_distance
+  , st_dwithin
+  , st_area
+  , st_perimeter
+  , st_length
+  , st_azimuth
+  , st_maxdistance
+  , st_angle
+  , st_closestpoint
+  , st_3dclosestpoint
+  , st_3ddistance
+  , st_distancesphere
+  , st_frechetdistance
+  , st_hausdorffdistance
+  , st_length2d
+  , st_3dlength
+  , st_longestline
+  , st_3dlongestline
+  , st_3dmaxdistance
+  , st_minimumclearance
+  , st_shortestline
+  , st_3dshortestline
+  ) where
+
+import Database.Esqueleto.Postgis.Geometry (Postgis, SpatialType(..))
+import Database.Esqueleto.Experimental (SqlExpr, Value)
+import Database.Esqueleto.Internal.Internal (unsafeSqlFunction)
+
+-- | calculate the distance between two points
+--   https://postgis.net/docs/ST_Distance.html
+st_distance ::
+  SqlExpr (Value (Postgis spatialType a)) ->
+  SqlExpr (Value (Postgis spatialType a)) ->
+  SqlExpr (Value Double)
+st_distance a b =
+  unsafeSqlFunction "ST_distance" (a, b)
+
+-- | Returns true if the geometries are within a given distance
+--   https://postgis.net/docs/ST_DWithin.html
+st_dwithin ::
+  -- | geometry g1
+  SqlExpr (Value (Postgis spatialType a)) ->
+  -- | geometry g2
+  SqlExpr (Value (Postgis spatialType a)) ->
+  -- | distance of srid
+  SqlExpr (Value Double) ->
+  SqlExpr (Value Bool)
+st_dwithin a b c = unsafeSqlFunction "ST_DWithin" (a, b, c)
+
+-- | Returns the area of a geometry or geography.
+--   For geometry, area is in SRID units squared. For geography, area is in square meters.
+--   https://postgis.net/docs/ST_Area.html
+st_area ::
+  SqlExpr (Value (Postgis spatialType a)) ->
+  SqlExpr (Value Double)
+st_area a = unsafeSqlFunction "ST_Area" a
+
+-- | Returns the perimeter of a geometry or geography.
+--   For geometry, in SRID units. For geography, in meters.
+--   https://postgis.net/docs/ST_Perimeter.html
+st_perimeter ::
+  SqlExpr (Value (Postgis spatialType a)) ->
+  SqlExpr (Value Double)
+st_perimeter a = unsafeSqlFunction "ST_Perimeter" a
+
+-- | Returns the 2D length of a linear geometry or geography.
+--   For geometry, in SRID units. For geography, in meters.
+--   https://postgis.net/docs/ST_Length.html
+st_length ::
+  SqlExpr (Value (Postgis spatialType a)) ->
+  SqlExpr (Value Double)
+st_length a = unsafeSqlFunction "ST_Length" a
+
+-- | Returns the azimuth in radians of the segment from point A to point B.
+--   The angle is measured clockwise from north (positive Y axis).
+--   https://postgis.net/docs/ST_Azimuth.html
+st_azimuth ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_azimuth a b = unsafeSqlFunction "ST_Azimuth" (a, b)
+
+-- | Returns the maximum distance between two geometries.
+--   The distance of the furthest pair of points.
+--   https://postgis.net/docs/ST_MaxDistance.html
+st_maxdistance ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_maxdistance a b = unsafeSqlFunction "ST_MaxDistance" (a, b)
+
+-- | Returns the angle between two vectors defined by two 2-point linestrings.
+--   This is the 2-argument form of ST_Angle which expects linestring inputs.
+--   Passing point geometries will result in a NULL return value — use the
+--   3 or 4 point overloads in PostGIS directly if you need point-based angles.
+--   https://postgis.net/docs/ST_Angle.html
+st_angle ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_angle a b = unsafeSqlFunction "ST_Angle" (a, b)
+
+-- | Returns the 2D point on geom1 closest to geom2.
+--   https://postgis.net/docs/ST_ClosestPoint.html
+st_closestpoint ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_closestpoint a b = unsafeSqlFunction "ST_ClosestPoint" (a, b)
+
+-- | Returns the 3D point on geom1 closest to geom2.
+--   https://postgis.net/docs/ST_3DClosestPoint.html
+st_3dclosestpoint ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_3dclosestpoint a b = unsafeSqlFunction "ST_3DClosestPoint" (a, b)
+
+-- | Returns the 3D cartesian minimum distance between two geometries.
+--   https://postgis.net/docs/ST_3DDistance.html
+st_3ddistance ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_3ddistance a b = unsafeSqlFunction "ST_3DDistance" (a, b)
+
+-- | Returns the minimum distance in meters between two lon/lat geometries using a spherical model.
+--   https://postgis.net/docs/ST_DistanceSphere.html
+st_distancesphere ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_distancesphere a b = unsafeSqlFunction "ST_DistanceSphere" (a, b)
+
+-- | Returns the Frechet distance between two geometries.
+--   https://postgis.net/docs/ST_FrechetDistance.html
+st_frechetdistance ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_frechetdistance a b = unsafeSqlFunction "ST_FrechetDistance" (a, b)
+
+-- | Returns the Hausdorff distance between two geometries.
+--   https://postgis.net/docs/ST_HausdorffDistance.html
+st_hausdorffdistance ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_hausdorffdistance a b = unsafeSqlFunction "ST_HausdorffDistance" (a, b)
+
+-- | Returns the 2D length of a linear geometry.
+--   https://postgis.net/docs/ST_Length2D.html
+st_length2d ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_length2d a = unsafeSqlFunction "ST_Length2D" a
+
+-- | Returns the 3D length of a linear geometry.
+--   https://postgis.net/docs/ST_3DLength.html
+st_3dlength ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_3dlength a = unsafeSqlFunction "ST_3DLength" a
+
+-- | Returns the 2D longest line between two geometries.
+--   https://postgis.net/docs/ST_LongestLine.html
+st_longestline ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_longestline a b = unsafeSqlFunction "ST_LongestLine" (a, b)
+
+-- | Returns the 3D longest line between two geometries.
+--   https://postgis.net/docs/ST_3DLongestLine.html
+st_3dlongestline ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_3dlongestline a b = unsafeSqlFunction "ST_3DLongestLine" (a, b)
+
+-- | Returns the 3D maximum distance between two geometries.
+--   https://postgis.net/docs/ST_3DMaxDistance.html
+st_3dmaxdistance ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_3dmaxdistance a b = unsafeSqlFunction "ST_3DMaxDistance" (a, b)
+
+-- | Returns the minimum clearance of a geometry, the shortest distance between two distinct vertices.
+--   https://postgis.net/docs/ST_MinimumClearance.html
+st_minimumclearance ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double)
+st_minimumclearance a = unsafeSqlFunction "ST_MinimumClearance" a
+
+-- | Returns the 2D shortest line between two geometries.
+--   https://postgis.net/docs/ST_ShortestLine.html
+st_shortestline ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_shortestline a b = unsafeSqlFunction "ST_ShortestLine" (a, b)
+
+-- | Returns the 3D shortest line between two geometries.
+--   https://postgis.net/docs/ST_3DShortestLine.html
+st_3dshortestline ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_3dshortestline a b = unsafeSqlFunction "ST_3DShortestLine" (a, b)
diff --git a/src/Database/Esqueleto/Postgis/Processing.hs b/src/Database/Esqueleto/Postgis/Processing.hs
new file mode 100644
--- /dev/null
+++ b/src/Database/Esqueleto/Postgis/Processing.hs
@@ -0,0 +1,271 @@
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Database.Esqueleto.Postgis.Processing
+  ( st_union
+  , st_unions
+  , st_centroid
+  , st_buffer
+  , st_convexhull
+  , st_envelope
+  , st_pointonsurface
+  , st_intersection
+  , st_difference
+  , st_symdifference
+  , st_unaryunion
+  , st_split
+  , st_node
+  , st_buildarea
+  , st_chaikinsmoothing
+  , st_concavehull
+  , st_delaunaytriangles
+  , st_generatepoints
+  , st_geometricmedian
+  , st_linemerge
+  , st_minimumboundingcircle
+  , st_offsetcurve
+  , st_reduceprecision
+  , st_sharedpaths
+  , st_simplify
+  , st_simplifypreservetopology
+  , st_voronoilines
+  , st_voronoipolygons
+  ) where
+
+import Data.Proxy
+import Database.Esqueleto.Postgis.Geometry (Postgis, SpatialType(..), HasPgType(..))
+import Database.Esqueleto.Experimental (SqlExpr, Value)
+import Database.Esqueleto.Internal.Internal (unsafeSqlFunction, unsafeSqlCastAs)
+
+-- | allows union of geometries, eg group a bunch together,
+--   https://postgis.net/docs/ST_Union.html
+--   for example:
+--
+-- @
+--  mCombined <- selectOne $ do
+--    grid <- from $ table @Grid
+--    pure $ st_union $ grid ^. GridGeom
+--
+--
+--  select $  do
+--    unit <- from $ table @Unit
+--    forM_ mCombined $ \combined ->
+--      where_ $ (unit ^. UnitGeom) `st_intersects` (val $ unValue combined)
+--    pure unit
+-- @
+st_union ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_union a = unsafeSqlFunction "ST_union" a
+
+st_unions ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry 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.
+  unsafeSqlFunction "ST_union" ((unsafeSqlCastAs casted a), (unsafeSqlCastAs casted b))
+  where
+    casted = (pgType $ Proxy @'Geometry)
+
+-- | Returns the geometric center (centroid) of a geometry.
+--   For polygons this is the center of mass of the surface.
+--   https://postgis.net/docs/ST_Centroid.html
+st_centroid ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_centroid a = unsafeSqlFunction "ST_Centroid" a
+
+-- | Returns a geometry that is the buffer of the input geometry at a given distance.
+--   The buffer is a polygon expanded (or shrunk if negative) from the input by the distance.
+--   https://postgis.net/docs/ST_Buffer.html
+st_buffer ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_buffer a d = unsafeSqlFunction "ST_Buffer" (a, d)
+
+-- | Returns the convex hull of a geometry.
+--   The smallest convex polygon that contains the input geometry.
+--   https://postgis.net/docs/ST_ConvexHull.html
+st_convexhull ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_convexhull a = unsafeSqlFunction "ST_ConvexHull" a
+
+-- | Returns the bounding box of a geometry as a geometry (polygon or point).
+--   A minimal axis-aligned rectangle that fully contains the input.
+--   https://postgis.net/docs/ST_Envelope.html
+st_envelope ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_envelope a = unsafeSqlFunction "ST_Envelope" a
+
+-- | Returns a point guaranteed to lie on the surface of the geometry.
+--   Unlike 'st_centroid', the result is always on or inside the geometry.
+--   https://postgis.net/docs/ST_PointOnSurface.html
+st_pointonsurface ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_pointonsurface a = unsafeSqlFunction "ST_PointOnSurface" a
+
+-- | Returns the shared portion of two geometries (their intersection).
+--   The result geometry contains only the area common to both inputs.
+--   https://postgis.net/docs/ST_Intersection.html
+st_intersection ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_intersection a b = unsafeSqlFunction "ST_Intersection" (a, b)
+
+-- | Returns the part of geometry A that does not intersect geometry B.
+--   Computes the geometric difference: A minus the shared area of A and B.
+--   https://postgis.net/docs/ST_Difference.html
+st_difference ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_difference a b = unsafeSqlFunction "ST_Difference" (a, b)
+
+-- | Returns the portions of A and B that do not intersect (symmetric difference).
+--   https://postgis.net/docs/ST_SymDifference.html
+st_symdifference ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_symdifference a b = unsafeSqlFunction "ST_SymDifference" (a, b)
+
+-- | Computes the union of a single geometry (dissolves internal boundaries).
+--   https://postgis.net/docs/ST_UnaryUnion.html
+st_unaryunion ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_unaryunion a = unsafeSqlFunction "ST_UnaryUnion" a
+
+-- | Splits a geometry by another geometry, returning a collection of geometries.
+--   https://postgis.net/docs/ST_Split.html
+st_split ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_split a b = unsafeSqlFunction "ST_Split" (a, b)
+
+-- | Nodes a set of linestrings, adding intersection points.
+--   https://postgis.net/docs/ST_Node.html
+st_node ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_node a = unsafeSqlFunction "ST_Node" a
+
+-- | Creates an areal geometry formed by the constituent linework of the input geometry.
+--   https://postgis.net/docs/ST_BuildArea.html
+st_buildarea ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_buildarea a = unsafeSqlFunction "ST_BuildArea" a
+
+-- | Returns a smoothed version of the geometry using Chaikin's algorithm.
+--   https://postgis.net/docs/ST_ChaikinSmoothing.html
+st_chaikinsmoothing ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_chaikinsmoothing a = unsafeSqlFunction "ST_ChaikinSmoothing" a
+
+-- | Returns a concave hull of a geometry with a target percent of area.
+--   https://postgis.net/docs/ST_ConcaveHull.html
+st_concavehull ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) -> -- ^ target_percent
+  SqlExpr (Value (Postgis 'Geometry a))
+st_concavehull a d = unsafeSqlFunction "ST_ConcaveHull" (a, d)
+
+-- | Returns the Delaunay triangulation of a geometry's vertices.
+--   https://postgis.net/docs/ST_DelaunayTriangles.html
+st_delaunaytriangles ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_delaunaytriangles a = unsafeSqlFunction "ST_DelaunayTriangles" a
+
+-- | Generates pseudo-random points within a polygon.
+--   https://postgis.net/docs/ST_GeneratePoints.html
+st_generatepoints ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Int) -> -- ^ npoints
+  SqlExpr (Value (Postgis 'Geometry a))
+st_generatepoints a n = unsafeSqlFunction "ST_GeneratePoints" (a, n)
+
+-- | Returns the geometric median of a multipoint.
+--   https://postgis.net/docs/ST_GeometricMedian.html
+st_geometricmedian ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_geometricmedian a = unsafeSqlFunction "ST_GeometricMedian" a
+
+-- | Merges a collection of linestrings into a single (or fewer) linestrings.
+--   https://postgis.net/docs/ST_LineMerge.html
+st_linemerge ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_linemerge a = unsafeSqlFunction "ST_LineMerge" a
+
+-- | Returns the minimum bounding circle enclosing a geometry.
+--   https://postgis.net/docs/ST_MinimumBoundingCircle.html
+st_minimumboundingcircle ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_minimumboundingcircle a = unsafeSqlFunction "ST_MinimumBoundingCircle" a
+
+-- | Returns an offset line at a given distance from the input line.
+--   https://postgis.net/docs/ST_OffsetCurve.html
+st_offsetcurve ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_offsetcurve a d = unsafeSqlFunction "ST_OffsetCurve" (a, d)
+
+-- | Reduces the precision of coordinates in a geometry to a given grid size.
+--   https://postgis.net/docs/ST_ReducePrecision.html
+st_reduceprecision ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_reduceprecision a d = unsafeSqlFunction "ST_ReducePrecision" (a, d)
+
+-- | Returns a collection of shared paths between two linestring/multilinestring geometries.
+--   https://postgis.net/docs/ST_SharedPaths.html
+st_sharedpaths ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_sharedpaths a b = unsafeSqlFunction "ST_SharedPaths" (a, b)
+
+-- | Returns a simplified version of a geometry using the Douglas-Peucker algorithm.
+--   https://postgis.net/docs/ST_Simplify.html
+st_simplify ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_simplify a d = unsafeSqlFunction "ST_Simplify" (a, d)
+
+-- | Returns a simplified version of a geometry that preserves topology.
+--   https://postgis.net/docs/ST_SimplifyPreserveTopology.html
+st_simplifypreservetopology ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_simplifypreservetopology a d = unsafeSqlFunction "ST_SimplifyPreserveTopology" (a, d)
+
+-- | Returns the Voronoi diagram edges for a set of points.
+--   https://postgis.net/docs/ST_VoronoiLines.html
+st_voronoilines ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_voronoilines a = unsafeSqlFunction "ST_VoronoiLines" a
+
+-- | Returns the Voronoi diagram polygons for a set of points.
+--   https://postgis.net/docs/ST_VoronoiPolygons.html
+st_voronoipolygons ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a))
+st_voronoipolygons a = unsafeSqlFunction "ST_VoronoiPolygons" a
diff --git a/src/Database/Esqueleto/Postgis/Spatial.hs b/src/Database/Esqueleto/Postgis/Spatial.hs
new file mode 100644
--- /dev/null
+++ b/src/Database/Esqueleto/Postgis/Spatial.hs
@@ -0,0 +1,169 @@
+{-# LANGUAGE DataKinds #-}
+
+module Database.Esqueleto.Postgis.Spatial
+  ( st_contains
+  , st_intersects
+  , st_within
+  , st_touches
+  , st_crosses
+  , st_disjoint
+  , st_equals
+  , st_covers
+  , st_coveredby
+  , st_overlaps
+  , st_containsproperly
+  , st_3dintersects
+  , st_relate
+  , st_orderingequals
+  , st_dfullywithin
+  , st_pointinsidecircle
+  ) where
+
+import Database.Esqueleto.Postgis.Geometry (Postgis, SpatialType(..))
+import Database.Esqueleto.Experimental (SqlExpr, Value)
+import Database.Esqueleto.Internal.Internal (unsafeSqlFunction)
+import Data.Text (Text)
+
+-- | Returns TRUE if geometry A contains geometry B.
+--   https://postgis.net/docs/ST_Contains.html
+st_contains ::
+  -- | geom a
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  -- | geom b
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_contains a b = unsafeSqlFunction "ST_CONTAINS" (a, b)
+
+-- | Returns true if two geometries intersect.
+--   Geometries intersect if they have any point in common.
+--   https://postgis.net/docs/ST_Intersects.html
+st_intersects ::
+  SqlExpr (Value (Postgis spatialType a)) -> -- ^ geomA or geogA
+  SqlExpr (Value (Postgis spatialType a)) -> -- ^ geomB or geogB
+  SqlExpr (Value Bool)
+st_intersects a b = unsafeSqlFunction "ST_Intersects" (a, b)
+
+-- | Returns TRUE if geometry A is within geometry B.
+--   Tests if every point of A lies inside (interior or boundary of) B.
+--   The inverse of 'st_contains': @st_within a b == st_contains b a@.
+--   https://postgis.net/docs/ST_Within.html
+st_within ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_within a b = unsafeSqlFunction "ST_Within" (a, b)
+
+-- | Returns TRUE if geometry A touches geometry B.
+--   They have at least one boundary point in common, but no interior points.
+--   https://postgis.net/docs/ST_Touches.html
+st_touches ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_touches a b = unsafeSqlFunction "ST_Touches" (a, b)
+
+-- | Returns TRUE if geometry A crosses geometry B.
+--   They have some but not all interior points in common.
+--   https://postgis.net/docs/ST_Crosses.html
+st_crosses ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_crosses a b = unsafeSqlFunction "ST_Crosses" (a, b)
+
+-- | Returns TRUE if geometry A is disjoint from geometry B.
+--   They do not share any space together, the inverse of 'st_intersects'.
+--   https://postgis.net/docs/ST_Disjoint.html
+st_disjoint ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_disjoint a b = unsafeSqlFunction "ST_Disjoint" (a, b)
+
+-- | Returns TRUE if geometry A is spatially equal to geometry B.
+--   The geometries represent the same region of space regardless of vertex order.
+--   https://postgis.net/docs/ST_Equals.html
+st_equals ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_equals a b = unsafeSqlFunction "ST_Equals" (a, b)
+
+-- | Returns TRUE if geometry/geography A covers geometry/geography B.
+--   No point in B is outside A. Similar to 'st_contains' but does not distinguish boundary and interior.
+--   https://postgis.net/docs/ST_Covers.html
+st_covers ::
+  SqlExpr (Value (Postgis spatialType a)) ->
+  SqlExpr (Value (Postgis spatialType a)) ->
+  SqlExpr (Value Bool)
+st_covers a b = unsafeSqlFunction "ST_Covers" (a, b)
+
+-- | Returns TRUE if geometry/geography A is covered by geometry/geography B.
+--   No point in A is outside B. The inverse of 'st_covers'.
+--   https://postgis.net/docs/ST_CoveredBy.html
+st_coveredby ::
+  SqlExpr (Value (Postgis spatialType a)) ->
+  SqlExpr (Value (Postgis spatialType a)) ->
+  SqlExpr (Value Bool)
+st_coveredby a b = unsafeSqlFunction "ST_CoveredBy" (a, b)
+
+-- | Returns TRUE if geometry A overlaps geometry B.
+--   They share some space but neither contains the other entirely.
+--   https://postgis.net/docs/ST_Overlaps.html
+st_overlaps ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_overlaps a b = unsafeSqlFunction "ST_Overlaps" (a, b)
+
+-- | Returns TRUE if geometry A contains geometry B properly.
+--   B must lie entirely inside the interior of A (not touching the boundary).
+--   https://postgis.net/docs/ST_ContainsProperly.html
+st_containsproperly ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_containsproperly a b = unsafeSqlFunction "ST_ContainsProperly" (a, b)
+
+-- | Returns TRUE if two 3D geometries intersect.
+--   https://postgis.net/docs/ST_3DIntersects.html
+st_3dintersects ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_3dintersects a b = unsafeSqlFunction "ST_3DIntersects" (a, b)
+
+-- | Returns the DE-9IM intersection matrix string for two geometries.
+--   https://postgis.net/docs/ST_Relate.html
+st_relate ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Text)
+st_relate a b = unsafeSqlFunction "ST_Relate" (a, b)
+
+-- | Returns TRUE if two geometries are point-by-point equal in the same order.
+--   https://postgis.net/docs/ST_OrderingEquals.html
+st_orderingequals ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Bool)
+st_orderingequals a b = unsafeSqlFunction "ST_OrderingEquals" (a, b)
+
+-- | Returns TRUE if all of the geometries are within the specified distance of one another.
+--   https://postgis.net/docs/ST_DFullyWithin.html
+st_dfullywithin ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) ->
+  SqlExpr (Value Bool)
+st_dfullywithin a b d = unsafeSqlFunction "ST_DFullyWithin" (a, b, d)
+
+-- | Returns TRUE if the point geometry is inside the circle defined by center_x, center_y and radius.
+--   https://postgis.net/docs/ST_PointInsideCircle.html
+st_pointinsidecircle ::
+  SqlExpr (Value (Postgis 'Geometry a)) ->
+  SqlExpr (Value Double) -> -- ^ center_x
+  SqlExpr (Value Double) -> -- ^ center_y
+  SqlExpr (Value Double) -> -- ^ radius
+  SqlExpr (Value Bool)
+st_pointinsidecircle a cx cy r = unsafeSqlFunction "ST_PointInsideCircle" (a, cx, cy, r)
diff --git a/test/Test.hs b/test/Test.hs
--- a/test/Test.hs
+++ b/test/Test.hs
@@ -23,341 +23,1411 @@
 import Control.Monad.Trans.Resource (MonadThrow, ResourceT, runResourceT)
 import Data.LineString (LineString, makeLineString)
 import Data.LinearRing (LinearRing, makeLinearRing)
-import Data.List.NonEmpty (NonEmpty)
-import Data.Sequence (Seq)
-import qualified Data.Sequence as Seq
-import Data.Text
-import Database.Esqueleto.Experimental
-import Database.Esqueleto.Postgis
-import Database.Persist
-import Database.Persist.Postgresql
-  ( ConnectionString,
-    withPostgresqlConn,
-  )
-import Database.Persist.TH
-  ( mkMigrate,
-    mkPersist,
-    persistUpperCase,
-    share,
-    sqlSettings,
-  )
-import Hedgehog
-import qualified Hedgehog.Gen as Gen
-import qualified Hedgehog.Range as Range
-import Test.Tasty
-import Test.Tasty.HUnit
-import Test.Tasty.Hspec(testSpec)
-
-import qualified Ewkb.GeometrySpec
-import qualified Ewkb.LineSpec
-import qualified Ewkb.PointSpec
-
-import qualified Wkb.EndianSpec
-import qualified Wkb.GeometryCollectionSpec
-import qualified Wkb.GeometrySpec
-import qualified Wkb.LineSpec
-import qualified Wkb.PointSpec
-import qualified Wkb.PolygonSpec
-
-
-connString :: ConnectionString
-connString = "host=localhost port=5432 user=test dbname=test password=test"
-
--- Test schema
-share
-  [mkPersist sqlSettings, mkMigrate "migrateAll"]
-  [persistUpperCase|
-  Grid
-    geom (PostgisGeometry PointXY)
-    label Text
-
-  GeoGrid
-    geo (Postgis 'Geography PointXY)
-    label Text
-
-  Unit sql=unit
-    geom       (PostgisGeometry PointXY)
-    deriving Eq Show
-
-  Unityz sql=unityz
-    geom       (PostgisGeometry PointXYZ)
-    deriving Eq Show
-
-  Unityzm sql=unityzm
-    geom       (PostgisGeometry PointXYZM)
-    deriving Eq Show
-|]
-
-initializeDB ::
-  (MonadIO m) =>
-  SqlPersistT (ResourceT m) ()
-initializeDB = do
-  runMigration migrateAll
-
-runDB :: (forall m. (MonadIO m, MonadLogger m, MonadThrow m) => SqlPersistT (ResourceT m) a) -> IO a
-runDB act =
-  runStderrLoggingT
-    . runResourceT
-    . withPostgresqlConn connString
-    . runSqlConn
-    $ (initializeDB >> act >>= \ret -> transactionUndo >> return ret)
-
-main :: IO ()
-main = do
-  -- normally hspec discover does this, but we don't want to use
-  -- convuleted preprocessor,
-  -- it sometimes doesn't work and it's hard to understand why.
-  -- so I prefer to just import everything by hand.
-  -- also, go look at the source of hspec discover, it's kinda shocking.
-  hspecTrees <- sequence [
-    testSpec "Ewkb.GeometrySpec" Ewkb.GeometrySpec.spec,
-    testSpec "Ewkb.LineSpec" Ewkb.LineSpec.spec,
-    testSpec "Ewkb.PointSpec" Ewkb.PointSpec.spec,
-
-    testSpec "Wkb.EndianSpec" Wkb.EndianSpec.spec,
-    testSpec "Wkb.GeometryCollectionSpec" Wkb.GeometryCollectionSpec.spec,
-    testSpec "Wkb.GeometrySpec" Wkb.GeometrySpec.spec,
-    testSpec "Wkb.LineSpec" Wkb.LineSpec.spec,
-    testSpec "Wkb.PointSpec" Wkb.PointSpec.spec,
-    testSpec "Wkb.PolygonSpec" Wkb.PolygonSpec.spec
-
-    ]
-
-  defaultMain $
-    testGroup "all tests" $ postgisBindingsTests : hspecTrees
-
-test' :: Gen (PostgisGeometry PointXY) -> TestTree
-test' gen =
-  testCase "roundtrip xy geometry" $ do
-    someUnit <- Gen.sample (Unit <$> gen)
-    result <- runDB $ do
-      _ <- insert someUnit
-      selectList @(Unit) [] []
-    (entityVal <$> result) @?= [someUnit]
-
-testxyz :: Gen (PostgisGeometry PointXYZ) -> TestTree
-testxyz gen =
-  testCase "roundtrip xyz geometry" $ do
-    someUnit <- Gen.sample (Unityz <$> gen)
-    result <- runDB $ do
-      _ <- insert someUnit
-      selectList @(Unityz) [] []
-    (entityVal <$> result) @?= [someUnit]
-
-testxyzm :: Gen (PostgisGeometry PointXYZM) -> TestTree
-testxyzm gen =
-  testCase "roundtryp xyzm geometry" $ do
-    someUnit <- Gen.sample (Unityzm <$> gen)
-    result <- runDB $ do
-      _ <- insert someUnit
-      selectList @(Unityzm) [] []
-    (entityVal <$> result) @?= [someUnit]
-
-postgisBindingsTests :: TestTree
-postgisBindingsTests =
-  testGroup
-    "postgis binding tests"
-    [ testGroup "roundtrip_tests_xy" $
-        (test') <$> (genCollection genPointxy : genGeometry genPointxy),
-      testGroup "roundtrip tests xyz" $
-        (testxyz) <$> (genCollection genPointxyz : genGeometry genPointxyz),
-      testGroup "roundtrip tests xyzm" $
-        (testxyzm) <$> (genCollection genPointxyzm : genGeometry genPointxyzm),
-      testGroup "function bindings" $
-        [ testCase ("it_finds_the_one_unit_with st_contains") $ do
-            result <- runDB $ do
-              _ <-
-                insert $
-                  Unit
-                    { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [(PointXY 2 0)]
-                    }
-
-              selectOne $ do
-                unit <- from $ table @Unit
-                where_ $ unit ^. UnitGeom `st_contains` (val $ Point (PointXY 1 1))
-                pure countRows
-            unValue <$> result @?= (Just (1 :: Int)),
-          testCase ("it finds the one unit with intersection st_point") $ do
-            result <- runDB $ do
-              _ <-
-                insert $
-                  Unit
-                    { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [(PointXY 2 0)]
-                    }
-
-              selectOne $ do
-                unit <- from $ table @Unit
-                where_ $ unit ^. UnitGeom `st_contains` (st_point (val 1) (val 1))
-                pure countRows
-            unValue <$> result @?= (Just (1 :: Int)),
-          testCase ("it finds the one unit with intersection st_intersects") $ do
-            result <- runDB $ do
-              _ <-
-                insert $
-                  Unit
-                    { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [(PointXY 2 0)]
-                    }
-
-              selectOne $ do
-                unit <- from $ table @Unit
-                where_ $ unit ^. UnitGeom `st_intersects` (st_point (val 1) (val 1))
-                pure countRows
-            unValue <$> result @?= (Just (1 :: Int)),
-          testCase ("see if we can union in PG and then get out some Haskell") $ do
-            result <- runDB $ do
-              _ <-
-                insert $
-                  Grid
-                    { gridGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [(PointXY 2 0)],
-                      gridLabel = "x"
-                    }
-              _ <-
-                insert $
-                  Grid
-                    { gridGeom = Polygon $ makePolygon (PointXY 2 0) (PointXY 2 2) (PointXY 4 2) $ Seq.fromList [(PointXY 4 0)],
-                      gridLabel = "y"
-                    }
-
-              selectOne $ do
-                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 ("union_in_PG_and then get out some Haskell for geo") $ do
-            result <- runDB $ do
-              _ <-
-                insert $
-                  GeoGrid
-                    { geoGridGeo = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [(PointXY 2 0)],
-                      geoGridLabel = "x"
-                    }
-              _ <-
-                insert $
-                  GeoGrid
-                    { geoGridGeo = Polygon $ makePolygon (PointXY 2 0) (PointXY 2 2) (PointXY 4 2) $ Seq.fromList [(PointXY 4 0)],
-                      geoGridLabel = "y"
-                    }
-
-              selectOne $ do
-                grid <- from $ table @GeoGrid
-                pure $ st_transform_geometry $ st_union $ st_transform_geography mercator $ grid ^. GeoGridGeo
-            -- just delete the input because it keeps changing order
-            -- we just want to make sure something comes out
-            (() <$ result) @?= (Just ()),
-          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 ("st_distance@geom can distance PG and then get out some Haskell, doing it wrong with geometry") $ do
-            result <- runDB $ do
-              selectOne $
-                pure $ st_distance @'Geometry
-                    (point_v (-118.24) 34.05) -- LA
-                    (point_v (-74.00) 40.71) -- NYC
-            unValue <$> result @?= (Just  44.73849796316367), -- not 44km, but geometry does that
-          testCase ("st_distance@geography can distance PG and then get out some Haskell, doing it wrong with geometry") $ do
-            result <- runDB $ do
-              selectOne $
-                pure $ st_distance @'Geography
-                    (point_v (-118.24) 34.05) -- LA
-                    (point_v (-74.00) 40.71) -- NYC
-            unValue <$> result @?= (Just 3_944_735.82464902), -- correct! (in m)
-          testCase ("st_distance@geom can distance PG and then get out some Haskell, doing it wrong with geometry") $ do
-            result <- runDB $ do
-              selectOne $
-                pure $ st_distance @'Geometry
-                    (st_point (val (-118.24)) (val 34.05)) -- LA
-                    (st_point (val (-74.00)) (val 40.71)) -- NYC
-            unValue <$> result @?= (Just  44.73849796316367), -- not 44km, but geometry does that
-          testCase ("st_distance@geography can distance PG and then get out some Haskell, doing it wrong with geometry") $ do
-            result <- runDB $ do
-              selectOne $
-                pure $ st_distance @'Geography
-                    (st_point (val (-118.24)) (val 34.05)) -- LA
-                    (st_point (val (-74.00)) (val 40.71)) -- NYC
-            unValue <$> result @?= (Just 3_944_735.82464902), -- correct! (in m)
-          testCase ("see if we can get just the units in the polygons") $ do
-            result <- runDB $ do
-              _ <-
-                insert $
-                  Grid
-                    { gridGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [(PointXY 2 0)],
-                      gridLabel = "x"
-                    }
-              _ <-
-                insert $
-                  Grid
-                    { gridGeom = Polygon $ makePolygon (PointXY 2 0) (PointXY 2 2) (PointXY 4 2) $ Seq.fromList [(PointXY 4 0)],
-                      gridLabel = "y"
-                    }
-              _ <-
-                insert $
-                  Unit
-                    { unitGeom = point 1 1
-                    }
-              _ <-
-                insert $
-                  Unit
-                    { unitGeom = point 1 2
-                    }
-
-              _ <-
-                insert $
-                  Unit
-                    { unitGeom = point 2 2
-                    }
-              _ <-
-                insert $
-                  Unit
-                    { unitGeom = point 9 9
-                    }
-              _ <-
-                insert $
-                  Unit
-                    { unitGeom = point 10 10
-                    }
-
-              mCombined <- selectOne $ do
-                grid <- from $ table @Grid
-                pure $ st_union $ grid ^. GridGeom
-
-              select $ do
-                unit <- from $ table @Unit
-                forM_ mCombined $ \combined ->
-                  where_ $ (unit ^. UnitGeom) `st_intersects` (val $ unValue combined)
-                pure unit
-
-            entityVal <$> result @?= [Unit {unitGeom = Point (PointXY {_xyX = 1.0, _xyY = 1.0})}, Unit {unitGeom = Point (PointXY {_xyX = 1.0, _xyY = 2.0})}, Unit {unitGeom = Point (PointXY {_xyX = 2.0, _xyY = 2.0})}],
-          testCase ("st_dwithin finds it wihtin range") $ do
-            result <- runDB $ do
-              _ <-
-                insert $
-                  Unit
-                    { unitGeom = point 1 1
-                    }
-
-              select $ do
-                unit <- from $ table @Unit
-                where_ $ st_dwithin (unit ^. UnitGeom) (st_point (val 1) (val 0)) (val 1)
-                pure unit
-
-            entityVal <$> result @?= [Unit {unitGeom = Point (PointXY {_xyX = 1.0, _xyY = 1.0})} ],
-          testCase ("st_dwithin doesn't finds it out range") $ do
-            result <- runDB $ do
-              _ <-
-                insert $
-                  Unit
-                    { unitGeom = point 1 1
-                    }
-
-              select $ do
-                unit <- from $ table @Unit
-                where_ $ st_dwithin (unit ^. UnitGeom) (st_point (val 2) (val 0)) (val 1)
-                pure unit
-
-            entityVal <$> result @?= []
+import Data.List.NonEmpty (NonEmpty(..))
+import Data.Sequence (Seq)
+import qualified Data.Sequence as Seq
+import Data.Text (Text, isInfixOf)
+import qualified Data.Text as T
+import Database.Esqueleto.Experimental
+import Database.Esqueleto.Postgis
+import Database.Persist
+import Database.Persist.Postgresql
+  ( ConnectionString,
+    withPostgresqlConn,
+  )
+import Database.Persist.TH
+  ( mkMigrate,
+    mkPersist,
+    persistUpperCase,
+    share,
+    sqlSettings,
+  )
+import Hedgehog
+import qualified Hedgehog.Gen as Gen
+import qualified Hedgehog.Range as Range
+import Test.Tasty
+import Test.Tasty.HUnit
+import Test.Tasty.Hspec(testSpec)
+
+import qualified Ewkb.GeometrySpec
+import qualified Ewkb.LineSpec
+import qualified Ewkb.PointSpec
+
+import qualified Wkb.EndianSpec
+import qualified Wkb.GeometryCollectionSpec
+import qualified Wkb.GeometrySpec
+import qualified Wkb.LineSpec
+import qualified Wkb.PointSpec
+import qualified Wkb.PolygonSpec
+
+
+connString :: ConnectionString
+connString = "host=localhost port=5432 user=test dbname=test password=test"
+
+-- Test schema
+share
+  [mkPersist sqlSettings, mkMigrate "migrateAll"]
+  [persistUpperCase|
+  Grid
+    geom (PostgisGeometry PointXY)
+    label Text
+
+  GeoGrid
+    geo (Postgis 'Geography PointXY)
+    label Text
+
+  Unit sql=unit
+    geom       (PostgisGeometry PointXY)
+    deriving Eq Show
+
+  Unityz sql=unityz
+    geom       (PostgisGeometry PointXYZ)
+    deriving Eq Show
+
+  Unityzm sql=unityzm
+    geom       (PostgisGeometry PointXYZM)
+    deriving Eq Show
+|]
+
+initializeDB ::
+  (MonadIO m) =>
+  SqlPersistT (ResourceT m) ()
+initializeDB = do
+  runMigration migrateAll
+
+runDB :: (forall m. (MonadIO m, MonadLogger m, MonadThrow m) => SqlPersistT (ResourceT m) a) -> IO a
+runDB act =
+  runStderrLoggingT
+    . runResourceT
+    . withPostgresqlConn connString
+    . runSqlConn
+    $ (initializeDB >> act >>= \ret -> transactionUndo >> return ret)
+
+main :: IO ()
+main = do
+  -- normally hspec discover does this, but we don't want to use
+  -- convuleted preprocessor,
+  -- it sometimes doesn't work and it's hard to understand why.
+  -- so I prefer to just import everything by hand.
+  -- also, go look at the source of hspec discover, it's kinda shocking.
+  hspecTrees <- sequence [
+    testSpec "Ewkb.GeometrySpec" Ewkb.GeometrySpec.spec,
+    testSpec "Ewkb.LineSpec" Ewkb.LineSpec.spec,
+    testSpec "Ewkb.PointSpec" Ewkb.PointSpec.spec,
+
+    testSpec "Wkb.EndianSpec" Wkb.EndianSpec.spec,
+    testSpec "Wkb.GeometryCollectionSpec" Wkb.GeometryCollectionSpec.spec,
+    testSpec "Wkb.GeometrySpec" Wkb.GeometrySpec.spec,
+    testSpec "Wkb.LineSpec" Wkb.LineSpec.spec,
+    testSpec "Wkb.PointSpec" Wkb.PointSpec.spec,
+    testSpec "Wkb.PolygonSpec" Wkb.PolygonSpec.spec
+
+    ]
+
+  defaultMain $
+    testGroup "all tests" $ postgisBindingsTests : hspecTrees
+
+test' :: Gen (PostgisGeometry PointXY) -> TestTree
+test' gen =
+  testCase "roundtrip xy geometry" $ do
+    someUnit <- Gen.sample (Unit <$> gen)
+    result <- runDB $ do
+      _ <- insert someUnit
+      selectList @(Unit) [] []
+    (entityVal <$> result) @?= [someUnit]
+
+testxyz :: Gen (PostgisGeometry PointXYZ) -> TestTree
+testxyz gen =
+  testCase "roundtrip xyz geometry" $ do
+    someUnit <- Gen.sample (Unityz <$> gen)
+    result <- runDB $ do
+      _ <- insert someUnit
+      selectList @(Unityz) [] []
+    (entityVal <$> result) @?= [someUnit]
+
+testxyzm :: Gen (PostgisGeometry PointXYZM) -> TestTree
+testxyzm gen =
+  testCase "roundtryp xyzm geometry" $ do
+    someUnit <- Gen.sample (Unityzm <$> gen)
+    result <- runDB $ do
+      _ <- insert someUnit
+      selectList @(Unityzm) [] []
+    (entityVal <$> result) @?= [someUnit]
+
+postgisBindingsTests :: TestTree
+postgisBindingsTests =
+  testGroup
+    "postgis binding tests"
+    [ testGroup "roundtrip_tests_xy" $
+        (test') <$> (genCollection genPointxy : genGeometry genPointxy),
+      testGroup "roundtrip tests xyz" $
+        (testxyz) <$> (genCollection genPointxyz : genGeometry genPointxyz),
+      testGroup "roundtrip tests xyzm" $
+        (testxyzm) <$> (genCollection genPointxyzm : genGeometry genPointxyzm),
+      testGroup "function bindings" $
+        [ testCase ("it_finds_the_one_unit_with st_contains") $ do
+            result <- runDB $ do
+              _ <-
+                insert $
+                  Unit
+                    { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [(PointXY 2 0)]
+                    }
+
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ unit ^. UnitGeom `st_contains` (val $ Point (PointXY 1 1))
+                pure countRows
+            unValue <$> result @?= (Just (1 :: Int)),
+          testCase ("it finds the one unit with intersection st_point") $ do
+            result <- runDB $ do
+              _ <-
+                insert $
+                  Unit
+                    { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [(PointXY 2 0)]
+                    }
+
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ unit ^. UnitGeom `st_contains` (st_point (val 1) (val 1))
+                pure countRows
+            unValue <$> result @?= (Just (1 :: Int)),
+          testCase ("it finds the one unit with intersection st_intersects") $ do
+            result <- runDB $ do
+              _ <-
+                insert $
+                  Unit
+                    { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [(PointXY 2 0)]
+                    }
+
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ unit ^. UnitGeom `st_intersects` (st_point (val 1) (val 1))
+                pure countRows
+            unValue <$> result @?= (Just (1 :: Int)),
+          testCase ("see if we can union in PG and then get out some Haskell") $ do
+            result <- runDB $ do
+              _ <-
+                insert $
+                  Grid
+                    { gridGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [(PointXY 2 0)],
+                      gridLabel = "x"
+                    }
+              _ <-
+                insert $
+                  Grid
+                    { gridGeom = Polygon $ makePolygon (PointXY 2 0) (PointXY 2 2) (PointXY 4 2) $ Seq.fromList [(PointXY 4 0)],
+                      gridLabel = "y"
+                    }
+
+              selectOne $ do
+                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 ("union_in_PG_and then get out some Haskell for geo") $ do
+            result <- runDB $ do
+              _ <-
+                insert $
+                  GeoGrid
+                    { geoGridGeo = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [(PointXY 2 0)],
+                      geoGridLabel = "x"
+                    }
+              _ <-
+                insert $
+                  GeoGrid
+                    { geoGridGeo = Polygon $ makePolygon (PointXY 2 0) (PointXY 2 2) (PointXY 4 2) $ Seq.fromList [(PointXY 4 0)],
+                      geoGridLabel = "y"
+                    }
+
+              selectOne $ do
+                grid <- from $ table @GeoGrid
+                pure $ st_transform_geometry $ st_union $ st_transform_geography mercator $ grid ^. GeoGridGeo
+            -- just delete the input because it keeps changing order
+            -- we just want to make sure something comes out
+            (() <$ result) @?= (Just ()),
+          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 ("st_distance@geom can distance PG and then get out some Haskell, doing it wrong with geometry") $ do
+            result <- runDB $ do
+              selectOne $
+                pure $ st_distance @'Geometry
+                    (point_v (-118.24) 34.05) -- LA
+                    (point_v (-74.00) 40.71) -- NYC
+            unValue <$> result @?= (Just  44.73849796316367), -- not 44km, but geometry does that
+          testCase ("st_distance@geography can distance PG and then get out some Haskell, doing it wrong with geometry") $ do
+            result <- runDB $ do
+              selectOne $
+                pure $ st_distance @'Geography
+                    (point_v (-118.24) 34.05) -- LA
+                    (point_v (-74.00) 40.71) -- NYC
+            unValue <$> result @?= (Just 3_944_735.82464902), -- correct! (in m)
+          testCase ("st_distance@geom can distance PG and then get out some Haskell, doing it wrong with geometry") $ do
+            result <- runDB $ do
+              selectOne $
+                pure $ st_distance @'Geometry
+                    (st_point (val (-118.24)) (val 34.05)) -- LA
+                    (st_point (val (-74.00)) (val 40.71)) -- NYC
+            unValue <$> result @?= (Just  44.73849796316367), -- not 44km, but geometry does that
+          testCase ("st_distance@geography can distance PG and then get out some Haskell, doing it wrong with geometry") $ do
+            result <- runDB $ do
+              selectOne $
+                pure $ st_distance @'Geography
+                    (st_point (val (-118.24)) (val 34.05)) -- LA
+                    (st_point (val (-74.00)) (val 40.71)) -- NYC
+            unValue <$> result @?= (Just 3_944_735.82464902), -- correct! (in m)
+          testCase ("see if we can get just the units in the polygons") $ do
+            result <- runDB $ do
+              _ <-
+                insert $
+                  Grid
+                    { gridGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [(PointXY 2 0)],
+                      gridLabel = "x"
+                    }
+              _ <-
+                insert $
+                  Grid
+                    { gridGeom = Polygon $ makePolygon (PointXY 2 0) (PointXY 2 2) (PointXY 4 2) $ Seq.fromList [(PointXY 4 0)],
+                      gridLabel = "y"
+                    }
+              _ <-
+                insert $
+                  Unit
+                    { unitGeom = point 1 1
+                    }
+              _ <-
+                insert $
+                  Unit
+                    { unitGeom = point 1 2
+                    }
+
+              _ <-
+                insert $
+                  Unit
+                    { unitGeom = point 2 2
+                    }
+              _ <-
+                insert $
+                  Unit
+                    { unitGeom = point 9 9
+                    }
+              _ <-
+                insert $
+                  Unit
+                    { unitGeom = point 10 10
+                    }
+
+              mCombined <- selectOne $ do
+                grid <- from $ table @Grid
+                pure $ st_union $ grid ^. GridGeom
+
+              select $ do
+                unit <- from $ table @Unit
+                forM_ mCombined $ \combined ->
+                  where_ $ (unit ^. UnitGeom) `st_intersects` (val $ unValue combined)
+                pure unit
+
+            entityVal <$> result @?= [Unit {unitGeom = Point (PointXY {_xyX = 1.0, _xyY = 1.0})}, Unit {unitGeom = Point (PointXY {_xyX = 1.0, _xyY = 2.0})}, Unit {unitGeom = Point (PointXY {_xyX = 2.0, _xyY = 2.0})}],
+          testCase ("st_dwithin finds it wihtin range") $ do
+            result <- runDB $ do
+              _ <-
+                insert $
+                  Unit
+                    { unitGeom = point 1 1
+                    }
+
+              select $ do
+                unit <- from $ table @Unit
+                where_ $ st_dwithin (unit ^. UnitGeom) (st_point (val 1) (val 0)) (val 1)
+                pure unit
+
+            entityVal <$> result @?= [Unit {unitGeom = Point (PointXY {_xyX = 1.0, _xyY = 1.0})} ],
+          testCase ("st_dwithin doesn't finds it out range") $ do
+            result <- runDB $ do
+              _ <-
+                insert $
+                  Unit
+                    { unitGeom = point 1 1
+                    }
+
+              select $ do
+                unit <- from $ table @Unit
+                where_ $ st_dwithin (unit ^. UnitGeom) (st_point (val 2) (val 0)) (val 1)
+                pure unit
+
+            entityVal <$> result @?= [],
+
+          -- st_within: interior point is within polygon
+          testCase "st_within finds point inside polygon" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0] }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ (val $ Point (PointXY 1 1)) `st_within` (unit ^. UnitGeom)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_touches: adjacent polygons touch at shared edge
+          testCase "st_touches detects shared edge" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0] }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ (val $ Polygon $ makePolygon (PointXY 2 0) (PointXY 2 2) (PointXY 4 2) $ Seq.fromList [PointXY 4 0]) `st_touches` (unit ^. UnitGeom)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_crosses: line crossing through polygon
+          testCase "st_crosses detects line crossing polygon" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0] }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ (val $ Line $ makeLineString (PointXY (-1) 1) (PointXY 3 1) $ Seq.fromList []) `st_crosses` (unit ^. UnitGeom)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_disjoint: exterior point is disjoint from polygon
+          testCase "st_disjoint finds disjoint geometries" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0] }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ (val $ Point (PointXY 9 9)) `st_disjoint` (unit ^. UnitGeom)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_equals: polygon equals itself
+          testCase "st_equals detects equal geometries" $ do
+            let poly = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0]
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = poly }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ (unit ^. UnitGeom) `st_equals` (val poly)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_covers: polygon covers interior point
+          testCase "st_covers polygon covers interior point" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0] }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ (unit ^. UnitGeom) `st_covers` (val $ Point (PointXY 1 1))
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_coveredby: interior point is covered by polygon
+          testCase "st_coveredby point covered by polygon" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0] }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ (val $ Point (PointXY 1 1)) `st_coveredby` (unit ^. UnitGeom)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_overlaps: overlapping polygons
+          testCase "st_overlaps detects overlapping polygons" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0] }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ (val $ Polygon $ makePolygon (PointXY 1 1) (PointXY 1 3) (PointXY 3 3) $ Seq.fromList [PointXY 3 1]) `st_overlaps` (unit ^. UnitGeom)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_containsproperly: polygon properly contains interior point (not on boundary)
+          testCase "st_containsproperly finds interior point" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0] }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ (unit ^. UnitGeom) `st_containsproperly` (val $ Point (PointXY 1 1))
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_area: area of 2x2 square = 4.0
+          testCase "st_area computes polygon area" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just (4.0 :: Double),
+
+          -- st_perimeter: perimeter of 2x2 square = 8.0
+          testCase "st_perimeter computes polygon perimeter" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_perimeter (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just (8.0 :: Double),
+
+          -- st_length: length of line from (0,0) to (3,4) = 5.0
+          testCase "st_length computes line length" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_length (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 3 4) $ Seq.fromList [])
+            unValue <$> result @?= Just (5.0 :: Double),
+
+          -- st_azimuth: azimuth from origin to (1,0) = pi/2 (east)
+          testCase "st_azimuth computes angle" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_azimuth
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 0 0))
+                (val $ Point (PointXY 1 0))
+            -- azimuth to east (positive X) = pi/2
+            case unValue <$> result of
+              Just v -> assertBool "azimuth ~ pi/2" (abs (v - pi / 2) < 1e-10)
+              Nothing -> assertFailure "expected a result for st_azimuth",
+
+          -- st_maxdistance: max distance between two squares
+          testCase "st_maxdistance computes maximum distance" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_maxdistance
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 1) (PointXY 1 1) $ Seq.fromList [PointXY 1 0])
+                (val $ Point (PointXY 4 0))
+            -- max distance is from (0,1) to (4,0) = sqrt(16+1) = sqrt 17
+            case unValue <$> result of
+              Just v -> assertBool "maxdistance ~ sqrt 17" (abs (v - sqrt 17) < 1e-10)
+              Nothing -> assertFailure "expected a result for st_maxdistance",
+
+          -- st_x: X coordinate of point
+          testCase "st_x extracts X coordinate" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x (val @(Postgis 'Geometry PointXY) $ Point (PointXY 3.5 7.2))
+            unValue <$> result @?= Just (3.5 :: Double),
+
+          -- st_y: Y coordinate of point
+          testCase "st_y extracts Y coordinate" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_y (val @(Postgis 'Geometry PointXY) $ Point (PointXY 3.5 7.2))
+            unValue <$> result @?= Just (7.2 :: Double),
+
+          -- st_npoints: 2x2 square has 5 points (closed ring)
+          testCase "st_npoints counts vertices" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_npoints (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just (6 :: Int),
+
+          -- st_numgeometries: single polygon returns 1
+          testCase "st_numgeometries counts sub-geometries" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_numgeometries (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_dimension: polygon has dimension 2
+          testCase "st_dimension returns topological dimension" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_dimension (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just (2 :: Int),
+
+          -- st_issimple: simple line is simple
+          testCase "st_issimple detects simple geometry" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Line $ makeLineString (PointXY 0 0) (PointXY 3 4) $ Seq.fromList [] }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ st_issimple (unit ^. UnitGeom)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_isclosed: closed linestring (start=end)
+          testCase "st_isclosed detects closed linestring" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Line $ makeLineString (PointXY 0 0) (PointXY 1 1) $ Seq.fromList [PointXY 2 0, PointXY 0 0] }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ st_isclosed (unit ^. UnitGeom)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_isvalid: valid polygon is valid
+          testCase "st_isvalid detects valid geometry" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0] }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ st_isvalid (unit ^. UnitGeom)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_srid: default SRID is 0
+          testCase "st_srid returns spatial reference id" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_srid (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 1))
+            unValue <$> result @?= Just (0 :: Int),
+
+          -- st_centroid: centroid of 2x2 square at origin is (1,1)
+          testCase "st_centroid computes geometric center" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0] }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ (val $ Point (PointXY 1 1)) `st_equals` st_centroid (unit ^. UnitGeom)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_buffer: buffer of a point by 1 has area ~ pi
+          testCase "st_buffer expands geometry by distance" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_buffer (val @(Postgis 'Geometry PointXY) $ Point (PointXY 0 0)) (val 1)
+            case unValue <$> result of
+              Just v -> assertBool "buffer area ~ pi" (abs (v - pi) < 0.1)
+              Nothing -> assertFailure "expected a result for st_buffer",
+
+          -- st_convexhull: convex hull of L-shaped points is a triangle
+          testCase "st_convexhull computes convex hull" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_convexhull (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 4 0) $ Seq.fromList [PointXY 0 3])
+            -- triangle with base 4, height 3, area = 6
+            unValue <$> result @?= Just (6.0 :: Double),
+
+          -- st_envelope: bounding box of 2x2 square is itself
+          testCase "st_envelope computes bounding box" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_envelope (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just (4.0 :: Double),
+
+          -- st_pointonsurface: result is contained by the polygon
+          testCase "st_pointonsurface returns point on surface" $ do
+            let poly = Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0]
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = poly }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ (unit ^. UnitGeom) `st_contains` st_pointonsurface (unit ^. UnitGeom)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_intersection: intersection of overlapping 2x2 squares = 1x1 area
+          testCase "st_intersection computes shared area" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_intersection
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+                (val $ Polygon $ makePolygon (PointXY 1 1) (PointXY 1 3) (PointXY 3 3) $ Seq.fromList [PointXY 3 1])
+            unValue <$> result @?= Just (1.0 :: Double),
+
+          -- st_difference: A minus overlap
+          testCase "st_difference computes geometric difference" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_difference
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+                (val $ Polygon $ makePolygon (PointXY 1 1) (PointXY 1 3) (PointXY 3 3) $ Seq.fromList [PointXY 3 1])
+            -- 2x2 square area 4.0, minus 1x1 overlap = 3.0
+            unValue <$> result @?= Just (3.0 :: Double),
+
+          -- =================================================================
+          -- Batch 2: 95 new function bindings
+          -- =================================================================
+
+          -- Geometry Constructors (4)
+
+          -- st_collect: aggregate collect over two polygons
+          testCase "st_collect aggregates geometries" $ do
+            result <- runDB $ do
+              _ <- insert $ Grid
+                { gridGeom = Point (PointXY 1 1), gridLabel = "a" }
+              _ <- insert $ Grid
+                { gridGeom = Point (PointXY 2 2), gridLabel = "b" }
+              selectOne $ do
+                grid <- from $ table @Grid
+                pure $ st_numgeometries $ st_collect $ grid ^. GridGeom
+            unValue <$> result @?= Just (2 :: Int),
+
+          -- st_makeenvelope: create a box from coordinates
+          testCase "st_makeenvelope creates a rectangle" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_makeenvelope (val 0) (val 0) (val 3) (val 4) (val 0)
+            unValue <$> result @?= Just (12.0 :: Double),
+
+          -- st_makeline: create line from two points
+          testCase "st_makeline creates a line from two points" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_npoints $ st_makeline
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 0 0))
+                (val $ Point (PointXY 3 4))
+            unValue <$> result @?= Just (2 :: Int),
+
+          -- st_makepolygon_line: create polygon from closed linestring
+          testCase "st_makepolygon_line creates polygon from closed linestring" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_makepolygon_line
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 0 2) $ Seq.fromList [PointXY 2 2, PointXY 2 0, PointXY 0 0])
+            unValue <$> result @?= Just (4.0 :: Double),
+
+          -- Geometry Accessors (20)
+
+          -- st_boundary: boundary of a polygon is a linestring
+          testCase "st_boundary returns boundary" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_geometrytype $ st_boundary
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just ("ST_LineString" :: Text),
+
+          -- st_coorddim: 2D point has coordinate dimension 2
+          testCase "st_coorddim returns coordinate dimension" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_coorddim (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 1))
+            unValue <$> result @?= Just (2 :: Int),
+
+          -- st_endpoint: endpoint of a linestring
+          testCase "st_endpoint returns last point" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_endpoint
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 3 4) $ Seq.fromList [])
+            unValue <$> result @?= Just (3.0 :: Double),
+
+          -- st_exteriorring: exterior ring of a polygon
+          testCase "st_exteriorring returns exterior ring" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_npoints $ st_exteriorring
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            case unValue <$> result of
+              Just n -> assertBool "exterior ring has points" (n >= 4)
+              Nothing -> assertFailure "expected a result for st_exteriorring",
+
+          -- st_geometryn: get first geometry from a multipoint
+          testCase "st_geometryn gets Nth geometry" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_geometryn
+                (val @(Postgis 'Geometry PointXY) $ MultiPoint (PointXY 5 6 :| [PointXY 7 8]))
+                (val 1)
+            unValue <$> result @?= Just (5.0 :: Double),
+
+          -- st_geometrytype: point reports as ST_Point
+          testCase "st_geometrytype returns type string" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_geometrytype (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 1))
+            unValue <$> result @?= Just ("ST_Point" :: Text),
+
+          -- st_interiorringn: polygon without holes has no interior ring, tested with polygon with hole
+          testCase "st_interiorringn does not crash" $ do
+            -- A simple polygon has 0 interior rings; we just check this doesn't error
+            result <- runDB $ do
+              selectOne $ pure $ st_numinteriorrings
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just (0 :: Int),
+
+          -- st_iscollection: multipoint is a collection
+          testCase "st_iscollection detects collection types" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = MultiPoint (PointXY 1 1 :| [PointXY 2 2]) }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ st_iscollection (unit ^. UnitGeom)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_isempty: non-empty point is not empty
+          testCase "st_isempty detects non-empty geometry" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_isempty (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 1))
+            unValue <$> result @?= Just False,
+
+          -- st_ispolygonccw: check CCW orientation
+          testCase "st_ispolygonccw checks orientation" $ do
+            -- CCW polygon: (0,0) (2,0) (2,2) (0,2) is CCW
+            result <- runDB $ do
+              selectOne $ pure $ st_ispolygonccw
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 2 0) (PointXY 2 2) $ Seq.fromList [PointXY 0 2])
+            case unValue <$> result of
+              Just _ -> pure () -- just check it returns a bool
+              Nothing -> assertFailure "expected a result for st_ispolygonccw",
+
+          -- st_ispolygoncw: check CW orientation
+          testCase "st_ispolygoncw checks orientation" $ do
+            -- CW polygon: (0,0) (0,2) (2,2) (2,0) is CW
+            result <- runDB $ do
+              selectOne $ pure $ st_ispolygoncw
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            case unValue <$> result of
+              Just _ -> pure () -- just check it returns a bool
+              Nothing -> assertFailure "expected a result for st_ispolygoncw",
+
+          -- st_isring: closed simple linestring is a ring
+          testCase "st_isring detects ring" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Line $ makeLineString (PointXY 0 0) (PointXY 1 1) $ Seq.fromList [PointXY 2 0, PointXY 0 0] }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ st_isring (unit ^. UnitGeom)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_m: M coordinate of a 4D point
+          testCase "st_m extracts M coordinate" $ do
+            result <- runDB $ do
+              _ <- insert $ Unityzm { unityzmGeom = Point (PointXYZM 1 2 3 42) }
+              selectOne $ do
+                unit <- from $ table @Unityzm
+                pure $ st_m (unit ^. UnityzmGeom)
+            unValue <$> result @?= Just (42.0 :: Double),
+
+          -- st_ndims: 2D point has 2 dimensions
+          testCase "st_ndims returns number of dimensions" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_ndims (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 1))
+            unValue <$> result @?= Just (2 :: Int),
+
+          -- st_nrings: polygon with no holes has 1 ring
+          testCase "st_nrings counts rings" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_nrings
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_numinteriorrings: simple polygon has 0 interior rings
+          testCase "st_numinteriorrings counts interior rings" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_numinteriorrings
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just (0 :: Int),
+
+          -- st_numpoints: linestring with 2 points
+          testCase "st_numpoints counts linestring points" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_numpoints
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 3 4) $ Seq.fromList [])
+            unValue <$> result @?= Just (2 :: Int),
+
+          -- st_pointn: first point of a linestring
+          testCase "st_pointn gets Nth point" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_pointn
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 5 6) (PointXY 7 8) $ Seq.fromList [])
+                (val 1)
+            unValue <$> result @?= Just (5.0 :: Double),
+
+          -- st_startpoint: first point of a linestring
+          testCase "st_startpoint returns first point" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_startpoint
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 5 6) (PointXY 7 8) $ Seq.fromList [])
+            unValue <$> result @?= Just (5.0 :: Double),
+
+          -- st_z: Z coordinate of a 3D point
+          testCase "st_z extracts Z coordinate" $ do
+            result <- runDB $ do
+              _ <- insert $ Unityz { unityzGeom = Point (PointXYZ 1 2 99) }
+              selectOne $ do
+                unit <- from $ table @Unityz
+                pure $ st_z (unit ^. UnityzGeom)
+            unValue <$> result @?= Just (99.0 :: Double),
+
+          -- Geometry Editors (16)
+
+          -- st_addpoint: add point to linestring
+          testCase "st_addpoint adds a point to linestring" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_npoints $ st_addpoint
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 1 1) $ Seq.fromList [])
+                (val $ Point (PointXY 2 2))
+                (val 1)
+            unValue <$> result @?= Just (3 :: Int),
+
+          -- st_collectionextract: extract points from a collection
+          testCase "st_collectionextract extracts geometry type" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_numgeometries $ st_collectionextract
+                (val @(Postgis 'Geometry PointXY) $ Collection (Point (PointXY 1 1) :| [Line $ makeLineString (PointXY 0 0) (PointXY 1 1) $ Seq.fromList []]))
+                (val 1) -- 1 = POINT
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_flipcoordinates: swap X and Y
+          testCase "st_flipcoordinates swaps X and Y" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_flipcoordinates (val @(Postgis 'Geometry PointXY) $ Point (PointXY 3 7))
+            unValue <$> result @?= Just (7.0 :: Double),
+
+          -- st_force2d: force to 2D
+          testCase "st_force2d forces to 2D" $ do
+            result <- runDB $ do
+              _ <- insert $ Unityz { unityzGeom = Point (PointXYZ 1 2 3) }
+              selectOne $ do
+                unit <- from $ table @Unityz
+                pure $ st_ndims $ st_force2d (unit ^. UnityzGeom)
+            unValue <$> result @?= Just (2 :: Int),
+
+          -- st_force3d: force to 3D
+          testCase "st_force3d forces to 3D" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_ndims $ st_force3d (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 2))
+            unValue <$> result @?= Just (3 :: Int),
+
+          -- st_force4d: force to 4D
+          testCase "st_force4d forces to 4D" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_ndims $ st_force4d (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 2))
+            unValue <$> result @?= Just (4 :: Int),
+
+          -- st_forcecollection: wrap in collection
+          testCase "st_forcecollection wraps in collection" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_geometrytype $ st_forcecollection (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 1))
+            unValue <$> result @?= Just ("ST_GeometryCollection" :: Text),
+
+          -- st_forcepolygonccw: force CCW orientation
+          testCase "st_forcepolygonccw forces CCW" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_ispolygonccw $ st_forcepolygonccw
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just True,
+
+          -- st_forcepolygoncw: force CW orientation
+          testCase "st_forcepolygoncw forces CW" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_ispolygoncw $ st_forcepolygoncw
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just True,
+
+          -- st_multi: convert to multi-type
+          testCase "st_multi converts to multi-type" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_geometrytype $ st_multi (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 1))
+            unValue <$> result @?= Just ("ST_MultiPoint" :: Text),
+
+          -- st_normalize: normalize geometry
+          testCase "st_normalize normalizes geometry" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_normalize
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just (4.0 :: Double),
+
+          -- st_reverse: reverse vertex order
+          testCase "st_reverse reverses vertex order" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_startpoint $ st_reverse
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 3 4) $ Seq.fromList [])
+            unValue <$> result @?= Just (3.0 :: Double),
+
+          -- st_segmentize: add vertices to long segments
+          testCase "st_segmentize adds vertices" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_npoints $ st_segmentize
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 10 0) $ Seq.fromList [])
+                (val 3)
+            case unValue <$> result of
+              Just n -> assertBool "segmentize adds points" (n > 2)
+              Nothing -> assertFailure "expected a result for st_segmentize",
+
+          -- st_setpoint: replace a point in linestring
+          testCase "st_setpoint replaces a point" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_startpoint $ st_setpoint
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 1 1) $ Seq.fromList [])
+                (val 0)
+                (val $ Point (PointXY 99 99))
+            unValue <$> result @?= Just (99.0 :: Double),
+
+          -- st_snaptogrid: snap to grid
+          testCase "st_snaptogrid snaps coordinates" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_snaptogrid (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1.3 2.7)) (val 1)
+            unValue <$> result @?= Just (1.0 :: Double),
+
+          -- st_snap: snap vertices
+          testCase "st_snap snaps vertices" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_npoints $ st_snap
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 10 0) $ Seq.fromList [])
+                (val $ Point (PointXY 5 1))
+                (val 2)
+            case unValue <$> result of
+              Just n -> assertBool "snap adds or keeps points" (n >= 2)
+              Nothing -> assertFailure "expected a result for st_snap",
+
+          -- Geometry Validation (2)
+
+          -- st_makevalid: valid geometry stays the same area
+          testCase "st_makevalid makes geometry valid" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_makevalid
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just (4.0 :: Double),
+
+          -- st_isvalidreason: valid geometry reports "Valid Geometry"
+          testCase "st_isvalidreason reports validity" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_isvalidreason
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just ("Valid Geometry" :: Text),
+
+          -- SRS Functions (1)
+
+          -- st_setsrid: set SRID and read back
+          testCase "st_setsrid sets SRID" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_srid $ st_setsrid (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 1)) (val 4326)
+            unValue <$> result @?= Just (4326 :: Int),
+
+          -- Geometry Output (4)
+
+          -- st_astext: WKT of a point
+          testCase "st_astext returns WKT" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_astext (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 2))
+            unValue <$> result @?= Just ("POINT(1 2)" :: Text),
+
+          -- st_asgeojson: GeoJSON of a point
+          testCase "st_asgeojson returns GeoJSON" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_asgeojson (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 2))
+            case unValue <$> result of
+              Just t -> assertBool "contains coordinates" (isInfixOf "coordinates" t)
+              Nothing -> assertFailure "expected a result for st_asgeojson",
+
+          -- st_asewkt: EWKT of a point
+          testCase "st_asewkt returns EWKT" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_asewkt (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 2))
+            case unValue <$> result of
+              Just t -> assertBool "contains POINT" (isInfixOf "POINT" t)
+              Nothing -> assertFailure "expected a result for st_asewkt",
+
+          -- st_geohash: geohash of a point with SRID 4326
+          testCase "st_geohash returns geohash" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_geohash $ st_setsrid (val @(Postgis 'Geometry PointXY) $ Point (PointXY (-74) 40.7)) (val 4326)
+            case unValue <$> result of
+              Just t -> assertBool "geohash is non-empty" (not (T.null t))
+              Nothing -> assertFailure "expected a result for st_geohash",
+
+          -- Spatial Relationships (5)
+
+          -- st_3dintersects: 3D intersection
+          testCase "st_3dintersects checks 3D intersection" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_3dintersects
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 0 0))
+                (val $ Point (PointXY 0 0))
+            unValue <$> result @?= Just True,
+
+          -- st_relate: DE-9IM matrix of two equal points
+          testCase "st_relate returns DE-9IM matrix" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_relate
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 0 0))
+                (val $ Point (PointXY 0 0))
+            -- Two equal points: "0FFFFFFF2"
+            unValue <$> result @?= Just ("0FFFFFFF2" :: Text),
+
+          -- st_orderingequals: same geometry with same vertex order
+          testCase "st_orderingequals checks vertex order equality" $ do
+            let line = Line $ makeLineString (PointXY 0 0) (PointXY 1 1) $ Seq.fromList []
+            result <- runDB $ do
+              selectOne $ pure $ st_orderingequals
+                (val @(Postgis 'Geometry PointXY) line)
+                (val line)
+            unValue <$> result @?= Just True,
+
+          -- st_dfullywithin: both geometries fully within distance
+          testCase "st_dfullywithin checks full containment" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Point (PointXY 0 0) }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ st_dfullywithin (unit ^. UnitGeom) (val $ Point (PointXY 1 0)) (val 2)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- st_pointinsidecircle: point inside circle
+          testCase "st_pointinsidecircle checks point in circle" $ do
+            result <- runDB $ do
+              _ <- insert $ Unit { unitGeom = Point (PointXY 1 1) }
+              selectOne $ do
+                unit <- from $ table @Unit
+                where_ $ st_pointinsidecircle (unit ^. UnitGeom) (val 0) (val 0) (val 2)
+                pure countRows
+            unValue <$> result @?= Just (1 :: Int),
+
+          -- Measurement Functions (15)
+
+          -- st_angle: angle between two 2-point linestrings (PostGIS ST_Angle with 2 line args)
+          testCase "st_angle computes angle" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_angle
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 1 0) $ Seq.fromList [])
+                (val $ Line $ makeLineString (PointXY 0 0) (PointXY 0 1) $ Seq.fromList [])
+            case unValue <$> result of
+              Just v -> assertBool "angle is finite" (v >= 0 && v <= 2 * pi)
+              Nothing -> assertFailure "expected a result for st_angle",
+
+          -- st_closestpoint: closest point on line to external point
+          testCase "st_closestpoint finds closest point" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_y $ st_closestpoint
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 10 0) $ Seq.fromList [])
+                (val $ Point (PointXY 5 5))
+            unValue <$> result @?= Just (0.0 :: Double),
+
+          -- st_3dclosestpoint: 3D closest point
+          testCase "st_3dclosestpoint finds 3D closest point" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_3dclosestpoint
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 10 0) $ Seq.fromList [])
+                (val $ Point (PointXY 5 5))
+            unValue <$> result @?= Just (5.0 :: Double),
+
+          -- st_3ddistance: 3D distance between two points
+          testCase "st_3ddistance computes 3D distance" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_3ddistance
+                (val @(Postgis 'Geometry PointXYZ) $ Point (PointXYZ 0 0 0))
+                (val $ Point (PointXYZ 1 0 0))
+            unValue <$> result @?= Just (1.0 :: Double),
+
+          -- st_distancesphere: spherical distance
+          testCase "st_distancesphere computes spherical distance" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_distancesphere
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 0 0))
+                (val $ Point (PointXY 1 0))
+            case unValue <$> result of
+              Just v -> assertBool "spherical distance > 0" (v > 0)
+              Nothing -> assertFailure "expected a result for st_distancesphere",
+
+          -- st_frechetdistance: Frechet distance
+          testCase "st_frechetdistance computes Frechet distance" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_frechetdistance
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 10 0) $ Seq.fromList [])
+                (val $ Line $ makeLineString (PointXY 0 1) (PointXY 10 1) $ Seq.fromList [])
+            unValue <$> result @?= Just (1.0 :: Double),
+
+          -- st_hausdorffdistance: Hausdorff distance
+          testCase "st_hausdorffdistance computes Hausdorff distance" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_hausdorffdistance
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 10 0) $ Seq.fromList [])
+                (val $ Line $ makeLineString (PointXY 0 1) (PointXY 10 1) $ Seq.fromList [])
+            unValue <$> result @?= Just (1.0 :: Double),
+
+          -- st_length2d: 2D length of a line
+          testCase "st_length2d computes 2D length" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_length2d
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 3 4) $ Seq.fromList [])
+            unValue <$> result @?= Just (5.0 :: Double),
+
+          -- st_3dlength: 3D length of a line
+          testCase "st_3dlength computes 3D length" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_3dlength $ st_makeline
+                (val @(Postgis 'Geometry PointXYZ) $ Point (PointXYZ 0 0 0))
+                (val $ Point (PointXYZ 3 4 0))
+            unValue <$> result @?= Just (5.0 :: Double),
+
+          -- st_longestline: longest line between two geometries
+          testCase "st_longestline returns longest line" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_length $ st_longestline
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 0 0))
+                (val $ Point (PointXY 3 4))
+            unValue <$> result @?= Just (5.0 :: Double),
+
+          -- st_3dlongestline: 3D longest line
+          testCase "st_3dlongestline returns 3D longest line" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_length $ st_3dlongestline
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 0 0))
+                (val $ Point (PointXY 3 4))
+            unValue <$> result @?= Just (5.0 :: Double),
+
+          -- st_3dmaxdistance: 3D max distance
+          testCase "st_3dmaxdistance computes 3D max distance" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_3dmaxdistance
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 0 0))
+                (val $ Point (PointXY 3 4))
+            unValue <$> result @?= Just (5.0 :: Double),
+
+          -- st_minimumclearance: minimum clearance of polygon
+          testCase "st_minimumclearance computes minimum clearance" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_minimumclearance
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            case unValue <$> result of
+              Just v -> assertBool "clearance > 0" (v > 0)
+              Nothing -> assertFailure "expected a result for st_minimumclearance",
+
+          -- st_shortestline: shortest line between two geometries
+          testCase "st_shortestline returns shortest line" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_length $ st_shortestline
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 0 0))
+                (val $ Point (PointXY 3 4))
+            unValue <$> result @?= Just (5.0 :: Double),
+
+          -- st_3dshortestline: 3D shortest line
+          testCase "st_3dshortestline returns 3D shortest line" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_length $ st_3dshortestline
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 0 0))
+                (val $ Point (PointXY 3 4))
+            unValue <$> result @?= Just (5.0 :: Double),
+
+          -- Overlay Functions (4)
+
+          -- st_symdifference: symmetric difference
+          testCase "st_symdifference computes symmetric difference" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_symdifference
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+                (val $ Polygon $ makePolygon (PointXY 1 1) (PointXY 1 3) (PointXY 3 3) $ Seq.fromList [PointXY 3 1])
+            -- each 4, overlap 1, so symdiff = 4 + 4 - 2*1 = 6
+            unValue <$> result @?= Just (6.0 :: Double),
+
+          -- st_unaryunion: dissolve internal boundaries
+          testCase "st_unaryunion dissolves internal boundaries" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_unaryunion
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 2) (PointXY 2 2) $ Seq.fromList [PointXY 2 0])
+            unValue <$> result @?= Just (4.0 :: Double),
+
+          -- st_split: split a line at a point
+          testCase "st_split splits geometry" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_numgeometries $ st_split
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 10 0) $ Seq.fromList [])
+                (val $ Point (PointXY 5 0))
+            unValue <$> result @?= Just (2 :: Int),
+
+          -- st_node: node linestrings
+          testCase "st_node nodes linestrings" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_npoints $ st_node
+                (val @(Postgis 'Geometry PointXY) $ Multiline
+                  ( makeLineString (PointXY 0 0) (PointXY 10 10) (Seq.fromList [])
+                    :| [makeLineString (PointXY 0 10) (PointXY 10 0) (Seq.fromList [])]))
+            case unValue <$> result of
+              Just n -> assertBool "node adds intersection point" (n > 4)
+              Nothing -> assertFailure "expected a result for st_node",
+
+          -- Geometry Processing (15)
+
+          -- st_buildarea: build area from linestrings
+          testCase "st_buildarea builds area from linestrings" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_buildarea
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 0 2) $ Seq.fromList [PointXY 2 2, PointXY 2 0, PointXY 0 0])
+            unValue <$> result @?= Just (4.0 :: Double),
+
+          -- st_chaikinsmoothing: smoothing
+          testCase "st_chaikinsmoothing smooths geometry" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_npoints $ st_chaikinsmoothing
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 1 5) $ Seq.fromList [PointXY 5 2, PointXY 10 0])
+            case unValue <$> result of
+              Just n -> assertBool "smoothing adds points" (n > 4)
+              Nothing -> assertFailure "expected a result for st_chaikinsmoothing",
+
+          -- st_concavehull: concave hull
+          testCase "st_concavehull computes concave hull" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_concavehull
+                (val @(Postgis 'Geometry PointXY) $ MultiPoint (PointXY 0 0 :| [PointXY 10 0, PointXY 5 10, PointXY 5 5]))
+                (val 1.0)
+            case unValue <$> result of
+              Just v -> assertBool "concave hull has area" (v > 0)
+              Nothing -> assertFailure "expected a result for st_concavehull",
+
+          -- st_delaunaytriangles: Delaunay triangulation
+          testCase "st_delaunaytriangles computes triangulation" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_numgeometries $ st_delaunaytriangles
+                (val @(Postgis 'Geometry PointXY) $ MultiPoint (PointXY 0 0 :| [PointXY 10 0, PointXY 5 10, PointXY 5 5]))
+            case unValue <$> result of
+              Just n -> assertBool "triangulation produces triangles" (n >= 1)
+              Nothing -> assertFailure "expected a result for st_delaunaytriangles",
+
+          -- st_generatepoints: generate random points in polygon
+          testCase "st_generatepoints generates points" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_numgeometries $ st_generatepoints
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 100) (PointXY 100 100) $ Seq.fromList [PointXY 100 0])
+                (val 5)
+            unValue <$> result @?= Just (5 :: Int),
+
+          -- st_geometricmedian: geometric median of multipoint
+          testCase "st_geometricmedian computes median" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_geometricmedian
+                (val @(Postgis 'Geometry PointXY) $ MultiPoint (PointXY 0 0 :| [PointXY 10 0, PointXY 5 10]))
+            case unValue <$> result of
+              Just v -> assertBool "median X in range" (v > 0 && v < 10)
+              Nothing -> assertFailure "expected a result for st_geometricmedian",
+
+          -- st_linemerge: merge multilinestring into single linestring
+          testCase "st_linemerge merges lines" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_geometrytype $ st_linemerge
+                (val @(Postgis 'Geometry PointXY) $ Multiline
+                  ( makeLineString (PointXY 0 0) (PointXY 1 1) (Seq.fromList [])
+                    :| [makeLineString (PointXY 1 1) (PointXY 2 2) (Seq.fromList [])]))
+            unValue <$> result @?= Just ("ST_LineString" :: Text),
+
+          -- st_minimumboundingcircle: bounding circle
+          testCase "st_minimumboundingcircle computes bounding circle" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_minimumboundingcircle
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 2 0) $ Seq.fromList [])
+            case unValue <$> result of
+              -- radius = 1, area ~ pi
+              Just v -> assertBool "bounding circle area ~ pi" (abs (v - pi) < 0.1)
+              Nothing -> assertFailure "expected a result for st_minimumboundingcircle",
+
+          -- st_offsetcurve: offset a line
+          testCase "st_offsetcurve offsets a line" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_length $ st_offsetcurve
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 10 0) $ Seq.fromList [])
+                (val 1)
+            case unValue <$> result of
+              Just v -> assertBool "offset line has length" (v > 0)
+              Nothing -> assertFailure "expected a result for st_offsetcurve",
+
+          -- st_reduceprecision: reduce precision
+          testCase "st_reduceprecision reduces precision" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_reduceprecision
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1.23456 2.34567))
+                (val 0.1)
+            case unValue <$> result of
+              Just v -> assertBool "reduced precision" (abs (v - 1.2) < 0.15)
+              Nothing -> assertFailure "expected a result for st_reduceprecision",
+
+          -- st_sharedpaths: shared paths between two lines
+          testCase "st_sharedpaths finds shared paths" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_geometrytype $ st_sharedpaths
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 5 0) $ Seq.fromList [PointXY 10 0])
+                (val $ Line $ makeLineString (PointXY 3 0) (PointXY 5 0) $ Seq.fromList [PointXY 7 0])
+            unValue <$> result @?= Just ("ST_GeometryCollection" :: Text),
+
+          -- st_simplify: simplify geometry
+          testCase "st_simplify simplifies geometry" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_npoints $ st_simplify
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 1 0.001) $ Seq.fromList [PointXY 2 0, PointXY 3 0.001, PointXY 4 0])
+                (val 0.01)
+            case unValue <$> result of
+              Just n -> assertBool "simplify reduces points" (n < 5)
+              Nothing -> assertFailure "expected a result for st_simplify",
+
+          -- st_simplifypreservetopology: simplify preserving topology
+          testCase "st_simplifypreservetopology simplifies preserving topology" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_npoints $ st_simplifypreservetopology
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 1 0.001) $ Seq.fromList [PointXY 2 0, PointXY 3 0.001, PointXY 4 0])
+                (val 0.01)
+            case unValue <$> result of
+              Just n -> assertBool "simplifyPreserveTopology reduces points" (n < 5)
+              Nothing -> assertFailure "expected a result for st_simplifypreservetopology",
+
+          -- st_voronoilines: Voronoi diagram edges
+          testCase "st_voronoilines computes Voronoi edges" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_npoints $ st_voronoilines
+                (val @(Postgis 'Geometry PointXY) $ MultiPoint (PointXY 0 0 :| [PointXY 10 0, PointXY 5 10]))
+            case unValue <$> result of
+              Just n -> assertBool "voronoi has points" (n > 0)
+              Nothing -> assertFailure "expected a result for st_voronoilines",
+
+          -- st_voronoipolygons: Voronoi diagram polygons
+          testCase "st_voronoipolygons computes Voronoi polygons" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_numgeometries $ st_voronoipolygons
+                (val @(Postgis 'Geometry PointXY) $ MultiPoint (PointXY 0 0 :| [PointXY 10 0, PointXY 5 10]))
+            unValue <$> result @?= Just (3 :: Int),
+
+          -- Affine Transformations (5)
+
+          -- st_translate: translate point
+          testCase "st_translate moves geometry" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_translate
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 2))
+                (val 10)
+                (val 20)
+            unValue <$> result @?= Just (11.0 :: Double),
+
+          -- st_scale: scale geometry
+          testCase "st_scale scales geometry" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_scale
+                (val @(Postgis 'Geometry PointXY) $ Polygon $ makePolygon (PointXY 0 0) (PointXY 0 1) (PointXY 1 1) $ Seq.fromList [PointXY 1 0])
+                (val 3)
+                (val 2)
+            -- original area 1, scaled by 3*2 = 6
+            unValue <$> result @?= Just (6.0 :: Double),
+
+          -- st_rotate: rotate point 180 degrees (pi radians) around origin
+          testCase "st_rotate rotates geometry" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_rotate
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 0))
+                (val pi)
+            case unValue <$> result of
+              Just v -> assertBool "rotated X ~ -1" (abs (v - (-1)) < 1e-10)
+              Nothing -> assertFailure "expected a result for st_rotate",
+
+          -- st_rotatex: rotate around X axis
+          testCase "st_rotatex rotates around X axis" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_z $ st_rotatex
+                (val @(Postgis 'Geometry PointXYZ) $ Point (PointXYZ 0 1 0))
+                (val (pi / 2))
+            case unValue <$> result of
+              Just v -> assertBool "rotateX z ~ 1" (abs (v - 1) < 1e-10)
+              Nothing -> assertFailure "expected a result for st_rotatex",
+
+          -- st_rotatez: rotate around Z axis (same as st_rotate for 2D)
+          testCase "st_rotatez rotates around Z axis" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_rotatez
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 1 0))
+                (val pi)
+            case unValue <$> result of
+              Just v -> assertBool "rotateZ X ~ -1" (abs (v - (-1)) < 1e-10)
+              Nothing -> assertFailure "expected a result for st_rotatez",
+
+          -- Bounding Box (1)
+
+          -- st_expand: expand bounding box
+          testCase "st_expand expands bounding box" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_area $ st_envelope $ st_expand
+                (val @(Postgis 'Geometry PointXY) $ Point (PointXY 0 0))
+                (val 5)
+            -- point expanded by 5 in all directions => 10x10 box => area 100
+            unValue <$> result @?= Just (100.0 :: Double),
+
+          -- Linear Referencing (3)
+
+          -- st_lineinterpolatepoint: midpoint of a line
+          testCase "st_lineinterpolatepoint interpolates point" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_x $ st_lineinterpolatepoint
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 10 0) $ Seq.fromList [])
+                (val 0.5)
+            unValue <$> result @?= Just (5.0 :: Double),
+
+          -- st_linelocatepoint: locate point on line
+          testCase "st_linelocatepoint locates point on line" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_linelocatepoint
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 10 0) $ Seq.fromList [])
+                (val $ Point (PointXY 5 0))
+            unValue <$> result @?= Just (0.5 :: Double),
+
+          -- st_linesubstring: extract portion of line
+          testCase "st_linesubstring extracts substring" $ do
+            result <- runDB $ do
+              selectOne $ pure $ st_length $ st_linesubstring
+                (val @(Postgis 'Geometry PointXY) $ Line $ makeLineString (PointXY 0 0) (PointXY 10 0) $ Seq.fromList [])
+                (val 0.0)
+                (val 0.5)
+            unValue <$> result @?= Just (5.0 :: Double)
         ]
     ]