{-# LANGUAGE LambdaCase, ScopedTypeVariables #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
module Data.Geometry.Geos.Geometry (
Geometry(..)
, Point(..)
, LinearRing(..)
, LineString(..)
, Polygon(..)
, MultiPoint(..)
, MultiLineString(..)
, MultiPolygon(..)
, Some(..)
, SRID
, binaryPredicate
, convertGeometryFromRaw
, convertGeometryToRaw
, convertMultiPolygonFromRaw
, ensurePoint
, ensureLineString
, ensureLinearRing
, ensurePolygon
, ensureMultiPoint
, ensureMultiPolygon
, ensureMultiLineString
, interpolate
, interpolateNormalized
, project
, projectNormalized
, equalsExact
, equals
, area
, geometryLength
, distance
, hausdorffDistance
, nearestPoints
, withSomeGeometry
) where
import Data.Data
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 ((<*>))
type SRID = Maybe Int
data Some :: (* -> *) -> * where
Some :: f a -> Some f
withSomeGeometry :: Some Geometry -> (forall a . Geometry a -> b) -> b
withSomeGeometry (Some p) f = f p
instance Show (Some Geometry) where
show (Some a) = "Some (" <> show a <> ")"
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
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, Show, Eq, Data, Typeable)
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)
newtype LinearRing = LinearRing CoordinateSequence
deriving (Read, Ord, Show, Eq, Data, Typeable)
instance Sem.Semigroup LinearRing where
(<>) (LinearRing a) (LinearRing b) = LinearRing (a <> b)
instance Monoid LinearRing where
mempty = LinearRing V.empty
newtype LineString = LineString CoordinateSequence
deriving (Read, Ord, Show, Eq, Data, Typeable)
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)
newtype MultiPoint = MultiPoint (V.Vector Point)
deriving (Read, Ord, Show, Eq, Data, Typeable)
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)
newtype MultiPolygon = MultiPolygon (V.Vector Polygon)
deriving (Read, Ord, Show, Eq, Data, Typeable)
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
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 (LinearRing cs) s = do
csr <- RC.createEmptyCoordinateSequence len (dimensionsCoordinateSequence cs)
V.zipWithM_ (setCoordinateSequence csr) (V.enumFromN 0 len) cs
R.createLinearRing csr >>= R.setSRID s
where
len = V.length cs
convertLineStringToRaw :: R.Geometry b => LineString -> SRID -> Geos b
convertLineStringToRaw (LineString cs) s = do
csr <- RC.createEmptyCoordinateSequence len (dimensionsCoordinateSequence cs)
V.zipWithM_ (setCoordinateSequence csr) ( V.enumFromN 0 len) cs
R.createLineString csr >>= R.setSRID s
where
len = V.length cs
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
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 -> error "GeometryCollection currently unsupported"
ensurePoint :: Some Geometry -> Geometry Point
ensurePoint g = withSomeGeometry g $ \g' -> case g' of
PointGeometry _ _ -> g'
_ -> error "This geometry was expected to be a Point"
ensureLineString :: Some Geometry -> Geometry LineString
ensureLineString g = withSomeGeometry g $ \g' -> case g' of
LineStringGeometry _ _ -> g'
_ -> error "This geometry was expected to be a LineString"
ensureLinearRing :: Some Geometry -> Geometry LinearRing
ensureLinearRing g = withSomeGeometry g $ \g' -> case g' of
LinearRingGeometry _ _ -> g'
_ -> error "This geometry was expected to be a LinearRing"
ensurePolygon :: Some Geometry -> Geometry Polygon
ensurePolygon g = withSomeGeometry g $ \g' -> case g' of
PolygonGeometry _ _ -> g'
_ -> error "This geometry was expected to be a Polygon"
ensureMultiPoint :: Some Geometry -> Geometry MultiPoint
ensureMultiPoint g = withSomeGeometry g $ \p' -> case p' of
MultiPointGeometry _ _ -> p'
_ -> error "This geometry was expected to be a MultiPoint"
ensureMultiLineString :: Some Geometry -> Geometry MultiLineString
ensureMultiLineString g = withSomeGeometry g $ \p' -> case p' of
MultiLineStringGeometry _ _ -> p'
_ -> error "This geometry was expected to be a MultiLineString"
ensureMultiPolygon :: Some Geometry -> Geometry MultiPolygon
ensureMultiPolygon g = withSomeGeometry g $ \p' -> case p' of
MultiPolygonGeometry _ _ -> p'
_ -> error "This geometry was expected to be a MultiPolygon"
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 (\i -> convertPointFromRaw =<< getGeometryN g i )
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 (\i -> convertLineStringFromRaw =<< getGeometryN g i)
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 (\i -> convertPolygonFromRaw =<< getGeometryN g i)
area :: Geometry a -> Double
area g = runGeos $ do
r :: R.Geom <- convertGeometryToRaw g
d <- R.area r
return d
geometryLength :: Geometry a -> Double
geometryLength g = runGeos $ do
r :: R.Geom <- convertGeometryToRaw g
l <- R.geometryLength r
return l
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)