diff --git a/hgeos.cabal b/hgeos.cabal --- a/hgeos.cabal +++ b/hgeos.cabal @@ -1,5 +1,5 @@ name: hgeos -version: 0.1.2.0 +version: 0.1.3.0 synopsis: Simple Haskell bindings to GEOS C API description: Simple Haskell bindings to the @@ -37,13 +37,13 @@ test-suite hgeos-app type: exitcode-stdio-1.0 - hs-source-dirs: src/app + hs-source-dirs: src/test main-is: Main.hs default-language: Haskell2010 build-depends: MissingH , base >= 4.7 && < 5 , hgeos - other-modules: HighLevelAPI - , LowLevelAPI - , Main - , Sample + other-modules: GEOSTest.Arith + , GEOSTest.HighLevelAPI + , GEOSTest.LowLevelAPI + , GEOSTest.Sample diff --git a/src/app/HighLevelAPI.hs b/src/app/HighLevelAPI.hs deleted file mode 100644 --- a/src/app/HighLevelAPI.hs +++ /dev/null @@ -1,18 +0,0 @@ -module HighLevelAPI (demo) where - -import Data.Geolocation.GEOS - --- Demonstrates use of high-level API --- Lifetimes of GEOS objects are automatically managed by the context objects --- which guarantees that they are released when the context goes out of scope -demo :: IO () -demo = do - withGEOS $ \ctx -> do - reader <- mkReader ctx - g0 <- readGeometry reader "POLYGON (( 10 10, 10 20, 20 20, 20 10, 10 10 ))" - g1 <- readGeometry reader "POLYGON (( 11 11, 11 12, 12 12, 12 11, 11 11 ))" - g2 <- intersection g0 g1 - writer <- mkWriter ctx - str <- writeGeometry writer g2 - putStrLn str - putStrLn "HighLevelAPI.demo done" diff --git a/src/app/LowLevelAPI.hs b/src/app/LowLevelAPI.hs deleted file mode 100644 --- a/src/app/LowLevelAPI.hs +++ /dev/null @@ -1,39 +0,0 @@ -module LowLevelAPI (demo) where - -import Control.Exception -import Data.Geolocation.GEOS.Imports -import Foreign.C - --- Demonstrates direct use of imports --- Lifetimes of various GEOS objects, including readers, writers and --- geometries, must be managed explicitly by clients using explicit calls to --- various destroy/free functions -demo :: IO () -demo = do - wkt0 <- newCString "POLYGON (( 10 10, 10 20, 20 20, 20 10, 10 10 ))" - wkt1 <- newCString "POLYGON (( 11 11, 11 12, 12 12, 12 11, 11 11 ))" - - withGEOS $ \ctx -> do - withWKTReader ctx $ \reader -> do - withGeometry ctx reader wkt0 $ \g0 -> do - withGeometry ctx reader wkt1 $ \g1 -> do - bracket (c_GEOSIntersection_r ctx g0 g1) (c_GEOSGeom_destroy_r ctx) $ \g2 -> do - withWKTWriter ctx $ \writer -> do - str <- bracket - (c_GEOSWKTWriter_write_r ctx writer g2) - (c_GEOSFree_r_CString ctx) - peekCString - putStrLn str - putStrLn "LowLevelAPI.demo done" - where - withGEOS :: (GEOSContextHandle_t -> IO a) -> IO a - withGEOS = bracket c_initializeGEOSWithHandlers c_finishGEOS_r - withWKTReader :: GEOSContextHandle_t -> (GEOSWKTReaderPtr -> IO a) -> IO a - withWKTReader ctx = bracket (c_GEOSWKTReader_create_r ctx) (c_GEOSWKTReader_destroy_r ctx) - withWKTWriter :: GEOSContextHandle_t -> (GEOSWKTWriterPtr -> IO a) -> IO a - withWKTWriter ctx = bracket (c_GEOSWKTWriter_create_r ctx) (c_GEOSWKTWriter_destroy_r ctx) - withGeometry :: GEOSContextHandle_t -> GEOSWKTReaderPtr -> CString -> (GEOSGeometryPtr -> IO a) -> IO a - withGeometry ctx reader wkt = - bracket - (c_GEOSWKTReader_read_r ctx reader wkt) - (c_GEOSGeom_destroy_r ctx) diff --git a/src/app/Main.hs b/src/app/Main.hs deleted file mode 100644 --- a/src/app/Main.hs +++ /dev/null @@ -1,17 +0,0 @@ -{-# LANGUAGE RecordWildCards #-} - -module Main (main) where - -import Data.Geolocation.GEOS -import Foreign.C -import qualified HighLevelAPI -import qualified LowLevelAPI -import qualified Sample - -main :: IO () -main = do - v <- version - putStrLn v - LowLevelAPI.demo - HighLevelAPI.demo - Sample.demo diff --git a/src/app/Sample.hs b/src/app/Sample.hs deleted file mode 100644 --- a/src/app/Sample.hs +++ /dev/null @@ -1,150 +0,0 @@ -{-# LANGUAGE RecordWildCards #-} - -module Sample (demo) where - -import Control.Monad -import Data.Geolocation.GEOS -import Data.List -import Data.String.Utils -import Paths_hgeos -import Text.Printf - -printGeometry :: Writer -> Geometry -> IO () -printGeometry r g = writeGeometry r g >>= putStrLn - -type Coordinate = (Double, Double, Double) - -getXYZs :: CoordinateSequence -> IO (Maybe [Coordinate]) -getXYZs coordSeq = do - maybeSize <- getSize coordSeq - case maybeSize of - Nothing -> return Nothing - Just size -> do - xyzs <- forM [0..(size - 1)] $ \i -> do - (Just x) <- getX coordSeq i - (Just y) <- getY coordSeq i - (Just z) <- getZ coordSeq i - return (x, y, z) - return $ Just xyzs - -data Extent = Extent - { minX :: Double - , maxX :: Double - , minY :: Double - , maxY :: Double - , minZ :: Double - , maxZ :: Double - } deriving Show - -extent :: [Coordinate] -> Extent -extent ((x, y, z) : cs) = - let (minX, maxX, minY, maxY, minZ, maxZ) = - foldr (\(x', y', z') (minX, maxX, minY, maxY, minZ, maxZ) -> - (if x' < minX then x' else minX, - if x' > maxX then x' else maxX, - if y' < minY then y' else minY, - if y' > maxY then y' else maxY, - if z' < minZ then z' else minZ, - if z' > maxX then z' else maxZ)) (x, x, y, y, z, z) cs - in Extent minX maxX minY maxY minZ maxZ - -mfloor :: Double -> Double -> Double -mfloor m x = (fromInteger $ floor (x / m)) * m - -frange :: Double -> Double -> Double -> [Double] -frange lower upper step = - let count = (upper - lower) / step - in [lower + (fromIntegral i) * step | i <- [0..(round count - 1)]] - -type Longitude = Double -type Latitude = Double -type Resolution = Double - -mkSquare :: Reader -> Longitude -> Latitude -> Resolution -> IO Geometry -mkSquare reader longitude latitude resolution = - let points = [ - (longitude, latitude), - (longitude + resolution, latitude), - (longitude + resolution, latitude + resolution), - (longitude, latitude + resolution), - (longitude, latitude) - ] - wkt = printf "POLYGON ((%s))" (intercalate "," (map (\(a, b) -> printf "%f %f" a b) points)) - in readGeometry reader wkt - -getGeometries :: Geometry -> IO [Geometry] -getGeometries geometry = do - count <- getNumGeometries geometry - forM [0..(count - 1)] (getGeometry geometry) - -findBiggestPolygon :: Geometry -> IO Geometry -findBiggestPolygon geometry = do - geometries@(g : gs) <- getGeometries geometry - a : as <- sequence $ map area geometries - let (g', _) = foldr (\p@(_, a') biggest@(_, aBiggest) -> if a' > aBiggest then p else biggest) (g, a) (zip gs as) - return g' - -resolution :: Double -resolution = 1.0 - -processPolygon :: String -> Resolution -> Longitude -> Latitude -> Geometry -> IO () -processPolygon tableName resolution longitude latitude p = do - let pId = polygonId resolution longitude latitude - shell <- exteriorRing p - coordSeq <- coordinateSequence shell - (Just xyzs) <- getXYZs coordSeq - forM_ (zip [0..] xyzs) $ \(pointId, (x, y, _)) -> do - let s = printf - "INSERT INTO %s (polygon_id, point_id, longitude, latitude) VALUES ('%s', %d, %f, %f);\n" - tableName - pId - (pointId :: Int) - x - y - putStr s - -polygonId :: Resolution -> Longitude -> Latitude -> String -polygonId resolution longitude latitude = "id_" ++ formatValue longitude ++ "_" ++ formatValue latitude - where - formatValue :: Double -> String - formatValue = replace "-" "m" . replace "." "_" . (printf "%.2f") . mfloor resolution - --- A more involved example of use of API -demo :: IO () -demo = do - fileName <- getDataFileName "data/namibia.wkt" - wkt <- readFile fileName - withGEOS $ \ctx -> do - reader <- mkReader ctx - writer <- mkWriter ctx - let p = printGeometry writer - - country <- readGeometry reader wkt - env <- envelope country - shell <- exteriorRing env - coordSeq <- coordinateSequence shell - (Just xyzs) <- getXYZs coordSeq - forM_ xyzs $ \(x, y, z) -> print (x, y, z) - let Extent{..} = extent xyzs - mfloorRes = mfloor resolution - longitudeBegin = mfloorRes minX - longitudeEnd = mfloorRes maxX + resolution - latitudeBegin = mfloorRes minY - latitudeEnd = mfloorRes maxY + resolution - print longitudeBegin - print longitudeEnd - print latitudeBegin - print latitudeEnd - let longitudes = frange longitudeBegin longitudeEnd 1.0 - latitudes = frange latitudeBegin latitudeEnd 1.0 - forM_ [(i, j) | i <- longitudes, j <- latitudes] $ \(longitude, latitude) -> do - square <- mkSquare reader longitude latitude resolution - overlap <- intersection square country - x <- isEmpty overlap - unless x $ do - id <- geometryTypeId overlap - polygon <- case id of - MultiPolygon -> findBiggestPolygon overlap - Polygon -> return overlap - processPolygon "foo" resolution longitude latitude polygon - putStrLn "Sample.demo done" diff --git a/src/lib/Data/Geolocation/GEOS.hs b/src/lib/Data/Geolocation/GEOS.hs --- a/src/lib/Data/Geolocation/GEOS.hs +++ b/src/lib/Data/Geolocation/GEOS.hs @@ -12,9 +12,9 @@ For the low-level FFI bindings, see "Data.Geolocation.GEOS.Imports". - + - + -} {-# LANGUAGE RecordWildCards #-} @@ -23,14 +23,14 @@ ( Context () , CoordinateSequence () , Geometry () - , GeometryTypeId (..) + , GeometryType (..) , Reader () , Writer () , area - , coordinateSequence , envelope - , exteriorRing - , geometryTypeId + , geomTypeId + , getCoordSeq + , getExteriorRing , getGeometry , getNumGeometries , getSize @@ -73,7 +73,7 @@ data CoordinateSequence = CoordinateSequence ContextStateRef GEOSCoordSequencePtr -- |Represents a geometry type ID -data GeometryTypeId = +data GeometryType = Point | LineString | LinearRing | @@ -104,51 +104,67 @@ value <- peek valuePtr return $ Just (realToFrac value) --- |Returns a 'CoordinateSequence' from the supplied 'Geometry' -coordinateSequence :: Geometry -> IO CoordinateSequence -coordinateSequence (Geometry sr hGeometry) = do +checkAndDoNotTrack :: ContextStateRef -> (GEOSContextHandle_t -> IO GEOSGeometryPtr) -> IO (Maybe Geometry) +checkAndDoNotTrack sr f = do ContextState{..} <- readIORef sr - h <- c_GEOSGeom_getCoordSeq_r hCtx hGeometry - -- Do not track - --modifyIORef' sr (\p@ContextState{..} -> p { hCoordinateSequences = h : hCoordinateSequences }) - return $ CoordinateSequence sr h + h <- f hCtx + return $ if isNullPtr h + then Nothing + else Just $ Geometry sr h -doNotTrack :: ContextStateRef -> (GEOSContextHandle_t -> IO GEOSGeometryPtr) -> IO Geometry -doNotTrack sr f = do +checkAndTrack :: ContextStateRef -> (GEOSContextHandle_t -> IO GEOSGeometryPtr) -> IO (Maybe Geometry) +checkAndTrack sr f = do ContextState{..} <- readIORef sr h <- f hCtx - return $ Geometry sr h + if isNullPtr h + then return Nothing + else do + modifyIORef' sr $ (\p@ContextState{..} -> p { hGeometries = h : hGeometries }) + return $ Just (Geometry sr h) -- |Returns a 'Geometry' instance representing the envelope of the supplied -- 'Geometry' -envelope :: Geometry -> IO Geometry +envelope :: Geometry -> IO (Maybe Geometry) envelope (Geometry sr h) = - track sr (\hCtx -> c_GEOSEnvelope_r hCtx h) - --- |Returns a 'Geometry' instance representing the exterior ring of the --- supplied 'Geometry' -exteriorRing :: Geometry -> IO Geometry -exteriorRing (Geometry sr h) = - doNotTrack sr (\hCtx -> c_GEOSGetExteriorRing_r hCtx h) + checkAndTrack sr (\hCtx -> c_GEOSEnvelope_r hCtx h) --- |Returns type ID of a 'Geometry' instance -geometryTypeId :: Geometry -> IO GeometryTypeId -geometryTypeId (Geometry sr h) = do +-- |Returns type of a 'Geometry' instance +geomTypeId :: Geometry -> IO (Maybe GeometryType) +geomTypeId (Geometry sr h) = do ContextState{..} <- readIORef sr value <- c_GEOSGeomTypeId_r hCtx h - return $ toEnum (fromIntegral value) + return $ if value == -1 + then Nothing + else Just $ toEnum (fromIntegral value) +-- |Returns a 'CoordinateSequence' from the supplied 'Geometry' +getCoordSeq :: Geometry -> IO (Maybe CoordinateSequence) +getCoordSeq (Geometry sr hGeometry) = do + ContextState{..} <- readIORef sr + h <- c_GEOSGeom_getCoordSeq_r hCtx hGeometry + return $ if isNullPtr h + then Nothing + else Just $ CoordinateSequence sr h + +-- |Returns a 'Geometry' instance representing the exterior ring of the +-- supplied 'Geometry' +getExteriorRing :: Geometry -> IO (Maybe Geometry) +getExteriorRing (Geometry sr h) = + checkAndDoNotTrack sr (\hCtx -> c_GEOSGetExteriorRing_r hCtx h) + -- |Returns child 'Geometry' at given index -getGeometry :: Geometry -> Int -> IO Geometry +getGeometry :: Geometry -> Int -> IO (Maybe Geometry) getGeometry (Geometry sr h) index = - doNotTrack sr (\hCtx -> c_GEOSGetGeometryN_r hCtx h (fromIntegral index)) + checkAndDoNotTrack sr (\hCtx -> c_GEOSGetGeometryN_r hCtx h (fromIntegral index)) -- |Gets the number of geometries in a 'Geometry' instance -getNumGeometries :: Geometry -> IO Int +getNumGeometries :: Geometry -> IO (Maybe Int) getNumGeometries (Geometry sr h) = do ContextState{..} <- readIORef sr value <- c_GEOSGetNumGeometries_r hCtx h - return $ fromIntegral value + return $ if value == -1 + then Nothing + else Just $ fromIntegral value getOrdinate :: (GEOSContextHandle_t -> GEOSCoordSequencePtr -> CUInt -> Ptr CDouble -> IO CInt) -> CoordinateSequence -> Word -> IO (Maybe Double) @@ -188,16 +204,19 @@ -- |Returns a 'Geometry' instance representing the intersection of the two -- supplied 'Geometry' instances: -intersection :: Geometry -> Geometry -> IO Geometry +intersection :: Geometry -> Geometry -> IO (Maybe Geometry) intersection (Geometry sr0 h0) (Geometry sr1 h1) = - track sr0 (\hCtx -> c_GEOSIntersection_r hCtx h0 h1) + checkAndTrack sr0 (\hCtx -> c_GEOSIntersection_r hCtx h0 h1) -- |Returns value indicating if specified 'Geometry' instance is empty -isEmpty :: Geometry -> IO Bool +isEmpty :: Geometry -> IO (Maybe Bool) isEmpty (Geometry sr h) = do ContextState{..} <- readIORef sr value <- c_GEOSisEmpty_r hCtx h - return $ value /= 0 + return $ case value of + 0 -> Just False + 1 -> Just True + _ -> Nothing mkContext :: IO Context mkContext = do @@ -207,27 +226,33 @@ -- |Creates a reader used to deserialize 'Geometry' instances from -- -format text: -mkReader :: Context -> IO Reader +mkReader :: Context -> IO (Maybe Reader) mkReader (Context sr) = do ContextState{..} <- readIORef sr h <- c_GEOSWKTReader_create_r hCtx - modifyIORef' sr (\p@ContextState{..} -> p { hReaders = h : hReaders }) - return $ Reader sr h + if isNullPtr h + then return Nothing + else do + modifyIORef' sr (\p@ContextState{..} -> p { hReaders = h : hReaders }) + return $ Just (Reader sr h) -- |Creates a writer used to serialize 'Geometry' instances to -- -format text: -mkWriter :: Context -> IO Writer +mkWriter :: Context -> IO (Maybe Writer) mkWriter (Context sr) = do ContextState{..} <- readIORef sr h <- c_GEOSWKTWriter_create_r hCtx - modifyIORef' sr (\p@ContextState{..} -> p { hWriters = h : hWriters }) - return $ Writer sr h + if isNullPtr h + then return Nothing + else do + modifyIORef' sr (\p@ContextState{..} -> p { hWriters = h : hWriters }) + return $ Just (Writer sr h) -- |Deserializes a 'Geometry' instance from the given 'String' using the -- supplied 'Reader': -readGeometry :: Reader -> String -> IO Geometry +readGeometry :: Reader -> String -> IO (Maybe Geometry) readGeometry (Reader sr h) str = withCString str $ \cs -> do - track sr (\hCtx -> c_GEOSWKTReader_read_r hCtx h cs) + checkAndTrack sr (\hCtx -> c_GEOSWKTReader_read_r hCtx h cs) releaseContext :: Context -> IO () releaseContext (Context sr) = do @@ -238,13 +263,6 @@ mapM_ (c_GEOSWKTReader_destroy_r hCtx) hReaders c_finishGEOS_r hCtx -track :: ContextStateRef -> (GEOSContextHandle_t -> IO GEOSGeometryPtr) -> IO Geometry -track sr f = do - ContextState{..} <- readIORef sr - h <- f hCtx - modifyIORef' sr $ (\p@ContextState{..} -> p { hGeometries = h : hGeometries }) - return $ Geometry sr h - -- |Reports version of GEOS API version :: IO String version = c_GEOSversion >>= peekCString @@ -264,11 +282,10 @@ withGEOS = bracket mkContext releaseContext -- |Serializes a 'Geometry' instance to a 'String' using the supplied 'Writer': -writeGeometry :: Writer -> Geometry -> IO String +writeGeometry :: Writer -> Geometry -> IO (Maybe String) writeGeometry (Writer sr hWriter) (Geometry _ hGeometry) = do ContextState{..} <- readIORef sr - str <- bracket + bracket (c_GEOSWKTWriter_write_r hCtx hWriter hGeometry) (c_GEOSFree_r_CString hCtx) - peekCString - return str + (\cs -> if cs == nullPtr then return Nothing else Just <$> peekCString cs) diff --git a/src/lib/Data/Geolocation/GEOS/Imports.hs b/src/lib/Data/Geolocation/GEOS/Imports.hs --- a/src/lib/Data/Geolocation/GEOS/Imports.hs +++ b/src/lib/Data/Geolocation/GEOS/Imports.hs @@ -14,7 +14,7 @@ For the high-level API, see "Data.Geolocation.GEOS". - + -} module Data.Geolocation.GEOS.Imports @@ -23,6 +23,7 @@ , GEOSGeometryPtr () , GEOSWKTReaderPtr () , GEOSWKTWriterPtr () + , NullablePtr (isNullPtr) , c_GEOSArea_r , c_GEOSCoordSeq_destroy_r , c_GEOSCoordSeq_getSize_r @@ -53,20 +54,33 @@ import Foreign.C import Foreign.Ptr +-- |Determines if given pointer is null +class NullablePtr a where + -- |Returns @True@ if pointer is null, @False@ otherwise + isNullPtr :: a -> Bool + -- |Wraps @GEOSContextHandle_t@ newtype GEOSContextHandle_t = GEOSContextHandle_t (Ptr GEOSContextHandle_t) -- |Wraps @GEOSCoordSequence*@ newtype GEOSCoordSequencePtr = GEOSCoordSequencePtr (Ptr GEOSCoordSequencePtr) +instance NullablePtr GEOSCoordSequencePtr where + isNullPtr (GEOSCoordSequencePtr p) = p == nullPtr -- |Wraps @GEOSGeometry*@ newtype GEOSGeometryPtr = GEOSGeometryPtr (Ptr GEOSGeometryPtr) +instance NullablePtr GEOSGeometryPtr where + isNullPtr (GEOSGeometryPtr p) = p == nullPtr -- |Wraps @GEOSWKTReader*@ newtype GEOSWKTReaderPtr = GEOSWKTReaderPtr (Ptr GEOSWKTReaderPtr) +instance NullablePtr GEOSWKTReaderPtr where + isNullPtr (GEOSWKTReaderPtr p) = p == nullPtr -- |Wraps @GEOSWKTWriter*@ newtype GEOSWKTWriterPtr = GEOSWKTWriterPtr (Ptr GEOSWKTWriterPtr) +instance NullablePtr GEOSWKTWriterPtr where + isNullPtr (GEOSWKTWriterPtr p) = p == nullPtr -- |Wraps @GEOSArea_r@ foreign import ccall "GEOSArea_r" diff --git a/src/test/GEOSTest/Arith.hs b/src/test/GEOSTest/Arith.hs new file mode 100644 --- /dev/null +++ b/src/test/GEOSTest/Arith.hs @@ -0,0 +1,12 @@ +module GEOSTest.Arith + ( frange + , mfloor + ) where + +mfloor :: Double -> Double -> Double +mfloor m x = (fromInteger $ floor (x / m)) * m + +frange :: Double -> Double -> Double -> [Double] +frange lower upper step = + let count = (upper - lower) / step + in [lower + (fromIntegral i) * step | i <- [0..(round count - 1)]] diff --git a/src/test/GEOSTest/HighLevelAPI.hs b/src/test/GEOSTest/HighLevelAPI.hs new file mode 100644 --- /dev/null +++ b/src/test/GEOSTest/HighLevelAPI.hs @@ -0,0 +1,18 @@ +module GEOSTest.HighLevelAPI (demo) where + +import Data.Geolocation.GEOS + +-- Demonstrates use of high-level API +-- Lifetimes of GEOS objects are automatically managed by the context objects +-- which guarantees that they are released when the context goes out of scope +demo :: IO () +demo = do + withGEOS $ \ctx -> do + (Just reader) <- mkReader ctx + (Just g0) <- readGeometry reader "POLYGON (( 10 10, 10 20, 20 20, 20 10, 10 10 ))" + (Just g1) <- readGeometry reader "POLYGON (( 11 11, 11 12, 12 12, 12 11, 11 11 ))" + (Just g2) <- intersection g0 g1 + (Just writer) <- mkWriter ctx + (Just str) <- writeGeometry writer g2 + putStrLn str + putStrLn "HighLevelAPI.demo done" diff --git a/src/test/GEOSTest/LowLevelAPI.hs b/src/test/GEOSTest/LowLevelAPI.hs new file mode 100644 --- /dev/null +++ b/src/test/GEOSTest/LowLevelAPI.hs @@ -0,0 +1,39 @@ +module GEOSTest.LowLevelAPI (demo) where + +import Control.Exception +import Data.Geolocation.GEOS.Imports +import Foreign.C + +-- Demonstrates direct use of imports +-- Lifetimes of various GEOS objects, including readers, writers and +-- geometries, must be managed explicitly by clients using explicit calls to +-- various destroy/free functions +demo :: IO () +demo = do + wkt0 <- newCString "POLYGON (( 10 10, 10 20, 20 20, 20 10, 10 10 ))" + wkt1 <- newCString "POLYGON (( 11 11, 11 12, 12 12, 12 11, 11 11 ))" + + withGEOS $ \ctx -> do + withWKTReader ctx $ \reader -> do + withGeometry ctx reader wkt0 $ \g0 -> do + withGeometry ctx reader wkt1 $ \g1 -> do + bracket (c_GEOSIntersection_r ctx g0 g1) (c_GEOSGeom_destroy_r ctx) $ \g2 -> do + withWKTWriter ctx $ \writer -> do + str <- bracket + (c_GEOSWKTWriter_write_r ctx writer g2) + (c_GEOSFree_r_CString ctx) + peekCString + putStrLn str + putStrLn "LowLevelAPI.demo done" + where + withGEOS :: (GEOSContextHandle_t -> IO a) -> IO a + withGEOS = bracket c_initializeGEOSWithHandlers c_finishGEOS_r + withWKTReader :: GEOSContextHandle_t -> (GEOSWKTReaderPtr -> IO a) -> IO a + withWKTReader ctx = bracket (c_GEOSWKTReader_create_r ctx) (c_GEOSWKTReader_destroy_r ctx) + withWKTWriter :: GEOSContextHandle_t -> (GEOSWKTWriterPtr -> IO a) -> IO a + withWKTWriter ctx = bracket (c_GEOSWKTWriter_create_r ctx) (c_GEOSWKTWriter_destroy_r ctx) + withGeometry :: GEOSContextHandle_t -> GEOSWKTReaderPtr -> CString -> (GEOSGeometryPtr -> IO a) -> IO a + withGeometry ctx reader wkt = + bracket + (c_GEOSWKTReader_read_r ctx reader wkt) + (c_GEOSGeom_destroy_r ctx) diff --git a/src/test/GEOSTest/Sample.hs b/src/test/GEOSTest/Sample.hs new file mode 100644 --- /dev/null +++ b/src/test/GEOSTest/Sample.hs @@ -0,0 +1,149 @@ +{-# LANGUAGE RecordWildCards #-} + +module GEOSTest.Sample (demo) where + +import Control.Monad +import Data.Geolocation.GEOS +import Data.List +import Data.String.Utils +import GEOSTest.Arith +import Paths_hgeos +import Text.Printf + +type Coordinate = (Double, Double, Double) +type Latitude = Double +type Longitude = Double +type Resolution = Double + +printGeometry :: Writer -> Geometry -> IO () +printGeometry r g = do + maybeStr <- writeGeometry r g + case maybeStr of + Nothing -> putStrLn "(invalid)" + (Just str) -> putStrLn str + +getXYZs :: CoordinateSequence -> IO (Maybe [Coordinate]) +getXYZs coordSeq = do + maybeSize <- getSize coordSeq + case maybeSize of + Nothing -> return Nothing + Just size -> do + xyzs <- forM [0..(size - 1)] $ \i -> do + (Just x) <- getX coordSeq i + (Just y) <- getY coordSeq i + (Just z) <- getZ coordSeq i + return (x, y, z) + return $ Just xyzs + +data Extent = Extent + { minX :: Double + , maxX :: Double + , minY :: Double + , maxY :: Double + , minZ :: Double + , maxZ :: Double + } deriving Show + +extent :: [Coordinate] -> Extent +extent ((x, y, z) : cs) = + let (minX, maxX, minY, maxY, minZ, maxZ) = + foldr (\(x', y', z') (minX, maxX, minY, maxY, minZ, maxZ) -> + (if x' < minX then x' else minX, + if x' > maxX then x' else maxX, + if y' < minY then y' else minY, + if y' > maxY then y' else maxY, + if z' < minZ then z' else minZ, + if z' > maxX then z' else maxZ)) (x, x, y, y, z, z) cs + in Extent minX maxX minY maxY minZ maxZ + +mkSquare :: Reader -> Longitude -> Latitude -> Resolution -> IO Geometry +mkSquare reader longitude latitude resolution = do + let points = [ + (longitude, latitude), + (longitude + resolution, latitude), + (longitude + resolution, latitude + resolution), + (longitude, latitude + resolution), + (longitude, latitude) + ] + wkt = printf "POLYGON ((%s))" (intercalate "," (map (\(a, b) -> printf "%f %f" a b) points)) + (Just g) <- readGeometry reader wkt + return g + +getGeometries :: Geometry -> IO [Geometry] +getGeometries geometry = do + (Just count) <- getNumGeometries geometry + forM [0..(count - 1)] $ \i -> do + (Just g) <- getGeometry geometry i + return g + +findBiggestPolygon :: Geometry -> IO Geometry +findBiggestPolygon geometry = do + geometries@(g : gs) <- getGeometries geometry + a : as <- sequence $ map area geometries + let (g', _) = foldr (\p@(_, a') biggest@(_, aBiggest) -> if a' > aBiggest then p else biggest) (g, a) (zip gs as) + return g' + +resolution :: Double +resolution = 1.0 + +processPolygon :: String -> Resolution -> Longitude -> Latitude -> Geometry -> IO () +processPolygon tableName resolution longitude latitude p = do + let pId = polygonId resolution longitude latitude + (Just shell) <- getExteriorRing p + (Just coordSeq) <- getCoordSeq shell + (Just xyzs) <- getXYZs coordSeq + forM_ (zip [0..] xyzs) $ \(pointId, (x, y, _)) -> do + let s = printf + "INSERT INTO %s (polygon_id, point_id, longitude, latitude) VALUES ('%s', %d, %f, %f);\n" + tableName + pId + (pointId :: Int) + x + y + putStr s + +polygonId :: Resolution -> Longitude -> Latitude -> String +polygonId resolution longitude latitude = "id_" ++ formatValue longitude ++ "_" ++ formatValue latitude + where + formatValue :: Double -> String + formatValue = replace "-" "m" . replace "." "_" . (printf "%.2f") . mfloor resolution + +-- A more involved example of use of API +demo :: IO () +demo = do + fileName <- getDataFileName "data/namibia.wkt" + wkt <- readFile fileName + withGEOS $ \ctx -> do + (Just reader) <- mkReader ctx + (Just writer) <- mkWriter ctx + let p = printGeometry writer + + (Just country) <- readGeometry reader wkt + (Just env) <- envelope country + (Just shell) <- getExteriorRing env + (Just coordSeq) <- getCoordSeq shell + (Just xyzs) <- getXYZs coordSeq + forM_ xyzs $ \(x, y, z) -> print (x, y, z) + let Extent{..} = extent xyzs + mfloorRes = mfloor resolution + longitudeBegin = mfloorRes minX + longitudeEnd = mfloorRes maxX + resolution + latitudeBegin = mfloorRes minY + latitudeEnd = mfloorRes maxY + resolution + print longitudeBegin + print longitudeEnd + print latitudeBegin + print latitudeEnd + let longitudes = frange longitudeBegin longitudeEnd 1.0 + latitudes = frange latitudeBegin latitudeEnd 1.0 + forM_ [(i, j) | i <- longitudes, j <- latitudes] $ \(longitude, latitude) -> do + square <- mkSquare reader longitude latitude resolution + (Just overlap) <- intersection square country + (Just x) <- isEmpty overlap + unless x $ do + (Just t) <- geomTypeId overlap + polygon <- case t of + MultiPolygon -> findBiggestPolygon overlap + Polygon -> return overlap + processPolygon "foo" resolution longitude latitude polygon + putStrLn "Sample.demo done" diff --git a/src/test/Main.hs b/src/test/Main.hs new file mode 100644 --- /dev/null +++ b/src/test/Main.hs @@ -0,0 +1,16 @@ +{-# LANGUAGE RecordWildCards #-} + +module Main (main) where + +import Data.Geolocation.GEOS +import qualified GEOSTest.HighLevelAPI as HighLevelAPI +import qualified GEOSTest.LowLevelAPI as LowLevelAPI +import qualified GEOSTest.Sample as Sample + +main :: IO () +main = do + v <- version + putStrLn v + LowLevelAPI.demo + HighLevelAPI.demo + Sample.demo