{-# LANGUAGE LambdaCase, ScopedTypeVariables #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
module Data.Geometry.Geos.Geometry
( Geometry(..)
, GeometryConstructionError
, Point
, point
, LinearRing
, linearRing
, LineString
, lineString
, Polygon
, polygon
, MultiPoint
, multiPoint
, MultiLineString
, multiLineString
, MultiPolygon
, multiPolygon
, GeometryCollection
, geometryCollection
, Some(..)
, Coordinate
, coordinate2
, coordinate3
, SRID
, binaryPredicate
, convertGeometryFromRaw
, convertGeometryToRaw
, convertMultiPolygonFromRaw
, ensurePoint
, ensureLineString
, ensureLinearRing
, ensurePolygon
, ensureMultiPoint
, ensureMultiPolygon
, ensureMultiLineString
, ensureGeometryCollection
, interpolate
, interpolateNormalized
, project
, projectNormalized
, equalsExact
, equals
, area
, geometryLength
, distance
, hausdorffDistance
, nearestPoints
, withSomeGeometry
, mapSomeGeometry
)
where
import Data.Data
import Data.List (intercalate)
import Data.Semigroup as Sem
import qualified Data.Vector as V
import qualified Data.Geometry.Geos.Raw.Geometry
as R
import qualified Data.Geometry.Geos.Raw.CoordSeq
as RC
import Data.Geometry.Geos.Raw.Base
import Control.Monad
import Control.Applicative ( (<*>) )
import GHC.Generics ( Generic )
type SRID = Maybe Int
data Some :: (* -> *) -> * where
Some ::f a -> Some f
instance Eq (Some Geometry) where
(Some (PointGeometry p s)) == (Some (PointGeometry p2 s2)) = p == p2 && s == s2
(Some (PointGeometry _ _)) == (Some _) = False
(Some (LineStringGeometry g s)) == (Some (LineStringGeometry g2 s2)) = g == g2 && s == s2
(Some (LineStringGeometry _ _)) == (Some _) = False
(Some (PolygonGeometry g s)) == (Some (PolygonGeometry g2 s2)) = g == g2 && s == s2
(Some (PolygonGeometry _ _)) == (Some _) = False
(Some (MultiPointGeometry g s)) == (Some (MultiPointGeometry g2 s2)) = g == g2 && s == s2
(Some (MultiPointGeometry _ _)) == (Some _) = False
(Some (MultiLineStringGeometry g s)) == (Some (MultiLineStringGeometry g2 s2)) = g == g2 && s == s2
(Some (MultiLineStringGeometry _ _)) == (Some _) = False
(Some (MultiPolygonGeometry g s)) == (Some (MultiPolygonGeometry g2 s2)) = g == g2 && s == s2
(Some (MultiPolygonGeometry _ _)) == (Some _) = False
(Some (CollectionGeometry g s)) == (Some (CollectionGeometry g2 s2)) = g == g2 && s == s2
(Some (CollectionGeometry _ _)) == (Some _) = False
(Some (LinearRingGeometry g s)) == (Some (LinearRingGeometry g2 s2)) = g == g2 && s == s2
(Some (LinearRingGeometry _ _)) == (Some _) = False
withSomeGeometry :: Some Geometry -> (forall a . Geometry a -> b) -> b
withSomeGeometry (Some p) f = f p
mapSomeGeometry :: (forall a . Geometry a -> b) -> Some Geometry -> b
mapSomeGeometry f (Some p) = f p
instance Show (Some Geometry) where
show (Some a) = "Some (" <> show a <> ")"
data GeometryConstructionError
= InvalidLinearRing
| InvalidLineString
| InvalidPolygon
data Geometry a where
PointGeometry :: Point -> SRID -> Geometry Point
LineStringGeometry :: LineString -> SRID -> Geometry LineString
LinearRingGeometry :: LinearRing -> SRID -> Geometry LinearRing
PolygonGeometry :: Polygon -> SRID -> Geometry Polygon
MultiPointGeometry :: MultiPoint -> SRID -> Geometry MultiPoint
MultiLineStringGeometry :: MultiLineString -> SRID -> Geometry MultiLineString
MultiPolygonGeometry :: MultiPolygon -> SRID -> Geometry MultiPolygon
CollectionGeometry :: GeometryCollection -> SRID -> Geometry GeometryCollection
deriving instance Eq (Geometry a)
deriving instance Show (Geometry a)
data Coordinate =
Coordinate2 {-# UNPACK #-} !Double {-# UNPACK #-} !Double
| Coordinate3 {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double
deriving (Read, Ord, Eq, Data, Typeable, Generic)
instance Show Coordinate where
show (Coordinate2 x y) = "(" ++ show x ++ ", " ++ show y ++ ")"
show (Coordinate3 x y z) =
"(" ++ show x ++ ", " ++ show y ++ ", " ++ show z ++ ")"
coordinate2 :: Double -> Double -> Coordinate
coordinate2 = Coordinate2
coordinate3 :: Double -> Double -> Double -> Coordinate
coordinate3 = Coordinate3
dimensionsCoordinate :: Coordinate -> Int
dimensionsCoordinate = length . gmapQ (const ())
type CoordinateSequence = V.Vector Coordinate
dimensionsCoordinateSequence :: CoordinateSequence -> Int
dimensionsCoordinateSequence = dimensionsCoordinate . V.head
newtype Point = Point Coordinate
deriving (Read, Ord, Show, Eq, Data, Typeable, Generic)
point :: Coordinate -> Point
point = Point
newtype LinearRing = LinearRing {
coordinateSequenceLinearRing :: CoordinateSequence } deriving (Read, Ord, Eq, Data, Typeable, Generic)
instance Show LinearRing where
show (LinearRing vec) =
let inner = intercalate ", " $ show <$> V.toList vec
in "[" ++ inner ++ "]"
linearRing :: CoordinateSequence -> Either GeometryConstructionError LinearRing
linearRing vec
| 1 <= V.length vec && V.length vec < 4 = Left InvalidLinearRing
| V.head vec /= V.last vec = Left InvalidLinearRing
| otherwise = Right $ LinearRing vec
instance Sem.Semigroup LinearRing where
(<>) (LinearRing a) (LinearRing b) = LinearRing (a <> b)
instance Monoid LinearRing where
mempty = LinearRing V.empty
newtype LineString = LineString {
coordinateSequenceLineString :: CoordinateSequence } deriving (Read, Ord, Eq, Data, Typeable, Generic)
instance Show LineString where
show (LineString vec) =
let inner = intercalate ", " $ show <$> V.toList vec
in "[" ++ inner ++ "]"
lineString :: CoordinateSequence -> Either GeometryConstructionError LineString
lineString vec
| V.length vec == 1 = Left InvalidLineString
| V.length vec == 2 && V.last vec == V.head vec = Left InvalidLineString
| otherwise = Right $ LineString vec
instance Sem.Semigroup LineString where
(<>) (LineString a) (LineString b) = LineString (a <> b)
instance Monoid LineString where
mempty = LineString V.empty
newtype Polygon = Polygon (V.Vector LinearRing)
deriving (Read, Ord, Show, Eq, Data, Typeable, Generic)
polygon :: V.Vector LinearRing -> Either GeometryConstructionError Polygon
polygon vec =
let shell = maybe V.empty coordinateSequenceLinearRing $ vec V.!? 0
holes = coordinateSequenceLinearRing <$> if V.null vec then V.empty else V.tail vec
in
do
_ <- if V.null shell && V.any (not . V.null) holes then Left InvalidPolygon else Right ()
return . Polygon $ V.cons (LinearRing shell) (fmap LinearRing holes)
newtype MultiPoint = MultiPoint (V.Vector Point)
deriving (Read, Ord, Show, Eq, Data, Typeable, Generic)
multiPoint :: V.Vector Point -> MultiPoint
multiPoint = MultiPoint
instance Sem.Semigroup MultiPoint where
(<>) (MultiPoint a) (MultiPoint b) = MultiPoint (a <> b)
instance Monoid MultiPoint where
mempty = MultiPoint V.empty
newtype MultiLineString = MultiLineString (V.Vector LineString)
deriving (Read, Ord, Show, Eq, Data, Typeable, Generic)
multiLineString :: V.Vector LineString -> MultiLineString
multiLineString = MultiLineString
newtype MultiPolygon = MultiPolygon (V.Vector Polygon)
deriving (Read, Ord, Show, Eq, Data, Typeable, Generic)
multiPolygon :: V.Vector Polygon -> MultiPolygon
multiPolygon = MultiPolygon
newtype GeometryCollection = GeometryCollection (V.Vector (Some Geometry))
deriving (Show, Eq, Typeable, Generic)
geometryCollection :: V.Vector (Some Geometry) -> GeometryCollection
geometryCollection = GeometryCollection
project :: Geometry LineString -> Geometry Point -> Double
project g1 g2 = runGeos $ do
g1' :: R.Geom <- convertGeometryToRaw g1
g2' :: R.Geom <- convertGeometryToRaw g2
R.project g1' g2'
projectNormalized :: Geometry LineString -> Geometry Point -> Double
projectNormalized g1 g2 = runGeos $ do
g1' :: R.Geom <- convertGeometryToRaw g1
g2' :: R.Geom <- convertGeometryToRaw g2
R.projectNormalized g1' g2'
interpolate :: Geometry LineString -> Double -> Geometry Point
interpolate g d = runGeos $ do
g' :: R.Geom <- convertGeometryToRaw g
s <- R.getSRID g'
p <- convertPointFromRaw =<< R.interpolate g' ( realToFrac d)
return $ PointGeometry p s
interpolateNormalized :: Geometry LineString -> Double -> Geometry Point
interpolateNormalized g d = runGeos $ do
g' :: R.Geom <- convertGeometryToRaw g
s <- R.getSRID g'
p <- convertPointFromRaw =<< R.interpolateNormalized g' ( realToFrac d)
return $ PointGeometry p s
binaryPredicate
:: (R.GeomConst -> R.GeomConst -> Geos Bool)
-> Geometry a
-> Geometry b
-> Bool
binaryPredicate f g1 g2 =
runGeos . join $ (f <$> convertGeometryToRaw g1 <*> convertGeometryToRaw g2)
equals :: Geometry a -> Geometry a -> Bool
equals = binaryPredicate R.equals
equalsExact :: Geometry a -> Geometry a -> Double -> Bool
equalsExact g1 g2 d =
binaryPredicate (\g1' g2' -> R.equalsExact g1' g2' d) g1 g2
convertGeometryToRaw
:: (R.Geometry a, R.CoordSeqInput a ~ cb, RC.CoordinateSequence cb)
=> Geometry b
-> Geos a
convertGeometryToRaw = \case
PointGeometry pg s -> convertPointToRaw pg s
LineStringGeometry lsg s -> convertLineStringToRaw lsg s
LinearRingGeometry lg s -> convertLinearRingToRaw lg s
PolygonGeometry pg s -> convertPolygonToRaw pg s
MultiPointGeometry mp s -> convertMultiPointToRaw mp s
MultiLineStringGeometry ml s -> convertMultiLineStringToRaw ml s
MultiPolygonGeometry mp s -> convertMultiPolygonToRaw mp s
CollectionGeometry gc s -> convertGeometryCollectionToRaw gc s
convertSomeGeometryToRaw
:: (R.Geometry a, R.CoordSeqInput a ~ cb, RC.CoordinateSequence cb)
=> Some Geometry
-> Geos a
convertSomeGeometryToRaw sg = withSomeGeometry sg convertGeometryToRaw
convertPointToRaw :: R.Geometry b => Point -> Maybe Int -> Geos b
convertPointToRaw (Point c) s = do
cs <- RC.createEmptyCoordinateSequence 1 (dimensionsCoordinate c)
setCoordinateSequence cs 0 c
R.createPoint cs >>= R.setSRID s
convertLinearRingToRaw :: R.Geometry b => LinearRing -> SRID -> Geos b
convertLinearRingToRaw =
convertLineal coordinateSequenceLinearRing R.createLinearRing
convertLineStringToRaw :: R.Geometry b => LineString -> SRID -> Geos b
convertLineStringToRaw =
convertLineal coordinateSequenceLineString R.createLineString
convertLineal
:: R.Geometry b
=> (a -> CoordinateSequence)
-> (RC.CoordSeqConst -> Geos b)
-> a
-> SRID
-> Geos b
convertLineal getCoordSeq convertGeom geom srid = do
let coordseq = getCoordSeq geom
len = V.length coordseq
csr <- RC.createEmptyCoordinateSequence
len
(dimensionsCoordinateSequence coordseq)
V.zipWithM_ (setCoordinateSequence csr) (V.enumFromN 0 len) coordseq
convertGeom csr >>= R.setSRID srid
convertPolygonToRaw :: R.Geometry a => Polygon -> SRID -> Geos a
convertPolygonToRaw (Polygon lrs) s = do
ext <- convertLinearRingToRaw (V.head lrs) s
inn <- flip convertLinearRingToRaw s `V.mapM` V.tail lrs
R.createPolygon ext (V.toList inn) >>= R.setSRID s
convertMultiPointToRaw :: R.Geometry a => MultiPoint -> SRID -> Geos a
convertMultiPointToRaw (MultiPoint vp) s = do
vr <- flip convertPointToRaw s `V.mapM` vp
R.createMultiPoint (V.toList vr) >>= R.setSRID s
convertMultiLineStringToRaw :: R.Geometry a => MultiLineString -> SRID -> Geos a
convertMultiLineStringToRaw (MultiLineString vl) s = do
vr <- flip convertLineStringToRaw s `V.mapM` vl
R.createMultiLineString (V.toList vr) >>= R.setSRID s
convertMultiPolygonToRaw :: R.Geometry a => MultiPolygon -> SRID -> Geos a
convertMultiPolygonToRaw (MultiPolygon vp) s = do
vr <- flip convertPolygonToRaw s `V.mapM` vp
R.createMultiPolygon (V.toList vr) >>= R.setSRID s
convertGeometryCollectionToRaw :: R.Geometry a => GeometryCollection -> SRID -> Geos a
convertGeometryCollectionToRaw (GeometryCollection vp) s = do
vr <- convertSomeGeometryToRaw `V.mapM` vp
R.createMultiPolygon (V.toList vr) >>= R.setSRID s
setCoordinateSequence
:: RC.CoordinateSequence a => a -> Int -> Coordinate -> Geos ()
setCoordinateSequence cs i (Coordinate2 x y) =
RC.setCoordinateSequenceX cs i x >> RC.setCoordinateSequenceY cs i y
setCoordinateSequence cs i (Coordinate3 x y z) =
RC.setCoordinateSequenceX cs i x
>> RC.setCoordinateSequenceY cs i y
>> RC.setCoordinateSequenceZ cs i z
convertGeometryFromRaw
:: (R.Geometry a, R.CoordSeqInput a ~ cb, RC.CoordinateSequence cb)
=> a
-> Geos (Some Geometry)
convertGeometryFromRaw rg = do
s <- R.getSRID rg
tid <- R.getTypeId rg
case tid of
R.PointTypeId -> do
p <- convertPointFromRaw rg
return $ Some (PointGeometry p s)
R.LineStringTypeId -> do
l <- convertLineStringFromRaw rg
return $ Some (LineStringGeometry l s)
R.LinearRingTypeId -> do
l <- convertLinearRingFromRaw rg
return $ Some (LinearRingGeometry l s)
R.PolygonTypeId -> do
p <- convertPolygonFromRaw rg
return $ Some (PolygonGeometry p s)
R.MultiPointTypeId -> do
mp <- convertMultiPointFromRaw rg
return $ Some (MultiPointGeometry mp s)
R.MultiLineStringTypeId -> do
ml <- convertMultiLineStringFromRaw rg
return $ Some (MultiLineStringGeometry ml s)
R.MultiPolygonTypeId -> do
mp <- convertMultiPolygonFromRaw rg
return $ Some (MultiPolygonGeometry mp s)
R.GeometryCollectionTypeId -> do
gc <- convertGeometryCollectionFromRaw rg
return $ Some (CollectionGeometry gc s)
ensurePoint :: Some Geometry -> Maybe (Geometry Point)
ensurePoint = mapSomeGeometry $ \case
g@(PointGeometry _ _) -> Just g
_ -> Nothing
ensureLineString :: Some Geometry -> Maybe (Geometry LineString)
ensureLineString = mapSomeGeometry $ \case
g@(LineStringGeometry _ _) -> Just g
_ -> Nothing
ensureLinearRing :: Some Geometry -> Maybe (Geometry LinearRing)
ensureLinearRing = mapSomeGeometry $ \case
g@(LinearRingGeometry _ _) -> Just g
_ -> Nothing
ensurePolygon :: Some Geometry -> Maybe (Geometry Polygon)
ensurePolygon = mapSomeGeometry $ \case
g@(PolygonGeometry _ _) -> Just g
_ -> Nothing
ensureMultiPoint :: Some Geometry -> Maybe (Geometry MultiPoint)
ensureMultiPoint = mapSomeGeometry $ \case
g@(MultiPointGeometry _ _) -> Just g
_ -> Nothing
ensureMultiLineString :: Some Geometry -> Maybe (Geometry MultiLineString)
ensureMultiLineString = mapSomeGeometry $ \case
g@(MultiLineStringGeometry _ _) -> Just g
_ -> Nothing
ensureMultiPolygon :: Some Geometry -> Maybe (Geometry MultiPolygon)
ensureMultiPolygon = mapSomeGeometry $ \case
g@(MultiPolygonGeometry _ _) -> Just g
_ -> Nothing
ensureGeometryCollection :: Some Geometry -> Maybe (Geometry GeometryCollection)
ensureGeometryCollection = mapSomeGeometry $ \case
g@(CollectionGeometry _ _) -> Just g
_ -> Nothing
getPosition :: RC.CoordinateSequence a => a -> Int -> Geos Coordinate
getPosition cs i = do
dim <- RC.getCoordinateSequenceDimensions cs
x <- RC.getCoordinateSequenceX cs i
y <- RC.getCoordinateSequenceY cs i
z <- if dim == 3
then Just <$> RC.getCoordinateSequenceZ cs i
else return Nothing
case z of
Nothing -> return $ Coordinate2 x y
Just z' -> return $ Coordinate3 x y z'
convertPointFromRaw
:: (R.Geometry a, R.CoordSeqInput a ~ ca, RC.CoordinateSequence ca)
=> a
-> Geos Point
convertPointFromRaw g = do
cs <- R.getCoordinateSequence g
Point <$> getPosition cs 0
convertSequenceFromRaw
:: (R.Geometry a, R.CoordSeqInput a ~ ca, RC.CoordinateSequence ca)
=> a
-> Geos CoordinateSequence
convertSequenceFromRaw g = do
cs <- R.getCoordinateSequence g
size <- R.getNumCoordinates g
V.generateM size (getPosition cs)
convertLineStringFromRaw
:: (R.Geometry a, R.CoordSeqInput a ~ ca, RC.CoordinateSequence ca)
=> a
-> Geos LineString
convertLineStringFromRaw g = LineString <$> convertSequenceFromRaw g
convertLinearRingFromRaw
:: (R.Geometry a, R.CoordSeqInput a ~ ca, RC.CoordinateSequence ca)
=> a
-> Geos LinearRing
convertLinearRingFromRaw g = LinearRing <$> convertSequenceFromRaw g
convertPolygonFromRaw :: R.Geometry a => a -> Geos Polygon
convertPolygonFromRaw g = do
is <- R.getNumInteriorRings g
ext :: R.GeomConst <- R.getExteriorRing g
ins :: V.Vector R.GeomConst <- V.generateM is (R.getInteriorRingN g)
Polygon <$> convertLinearRingFromRaw `V.mapM` (ext `V.cons` ins)
getGeometryN :: R.Geometry a => a -> Int -> Geos R.GeomConst
getGeometryN = R.getGeometryN
convertMultiPointFromRaw
:: (R.Geometry a, R.CoordSeqInput a ~ ca, RC.CoordinateSequence ca)
=> a
-> Geos MultiPoint
convertMultiPointFromRaw g = do
ng <- R.getNumGeometries g
MultiPoint <$> V.generateM ng (convertPointFromRaw <=< getGeometryN g)
convertMultiLineStringFromRaw
:: (R.Geometry a, R.CoordSeqInput a ~ ca, RC.CoordinateSequence ca)
=> a
-> Geos MultiLineString
convertMultiLineStringFromRaw g = do
ng <- R.getNumGeometries g
MultiLineString
<$> V.generateM ng (convertLineStringFromRaw <=< getGeometryN g)
convertMultiPolygonFromRaw
:: (R.Geometry a, R.CoordSeqInput a ~ ca, RC.CoordinateSequence ca)
=> a
-> Geos MultiPolygon
convertMultiPolygonFromRaw g = do
ng <- R.getNumGeometries g
MultiPolygon <$> V.generateM ng (convertPolygonFromRaw <=< getGeometryN g)
convertGeometryCollectionFromRaw
:: (R.Geometry a, R.CoordSeqInput a ~ ca, RC.CoordinateSequence ca)
=> a
-> Geos GeometryCollection
convertGeometryCollectionFromRaw g = do
ng <- R.getNumGeometries g
GeometryCollection <$> V.generateM ng (convertGeometryFromRaw <=< getGeometryN g)
area :: Geometry a -> Double
area g = runGeos $ do
r :: R.Geom <- convertGeometryToRaw g
R.area r
geometryLength :: Geometry a -> Double
geometryLength g = runGeos $ do
r :: R.Geom <- convertGeometryToRaw g
R.geometryLength r
distance :: Geometry a -> Geometry a -> Double
distance p g = runGeos $ do
p' :: R.Geom <- convertGeometryToRaw p
g' :: R.Geom <- convertGeometryToRaw g
R.distance p' g'
hausdorffDistance :: Geometry a -> Geometry a -> Double
hausdorffDistance p g = runGeos $ do
p' :: R.Geom <- convertGeometryToRaw p
g' :: R.Geom <- convertGeometryToRaw g
R.hausdorffDistance p' g'
nearestPoints :: Geometry a -> Geometry a -> (Coordinate, Coordinate)
nearestPoints g1 g2 = runGeos $ do
g1' :: R.Geom <- convertGeometryToRaw g1
g2' :: R.Geom <- convertGeometryToRaw g2
cs :: RC.CoordSeq <- R.nearestPoints g1' g2'
p1 <- getPosition cs 0
p2 <- getPosition cs 1
return (p1, p2)