From 6888879712792c8c0bc88fb0611b92419972a6d0 Mon Sep 17 00:00:00 2001 From: Frank Staals <991345+noinia@users.noreply.github.com> Date: Fri, 16 Aug 2024 08:11:43 +0200 Subject: [PATCH] partial integration/support for geojson (#237) * some instances for the points * exporting functionality + cleaning up imports * some work on the examples * bunch of instances for the polygon type * some fiddling * more instances :) * uncommented the partial instance for now * qualified import * missing exe --- .../src/HGeometry/Foldable/Util.hs | 3 + hgeometry-examples/geojson/Main.hs | 43 +++ hgeometry-examples/hgeometry-examples.cabal | 16 + hgeometry-examples/skia/Main.hs | 3 +- hgeometry/geojson/src/HGeometry/GeoJSON.hs | 329 +++++++++++++++++- hgeometry/hgeometry.cabal | 7 +- 6 files changed, 397 insertions(+), 4 deletions(-) create mode 100644 hgeometry-examples/geojson/Main.hs diff --git a/hgeometry-combinatorial/src/HGeometry/Foldable/Util.hs b/hgeometry-combinatorial/src/HGeometry/Foldable/Util.hs index 27154027d..6ed4d65bd 100644 --- a/hgeometry-combinatorial/src/HGeometry/Foldable/Util.hs +++ b/hgeometry-combinatorial/src/HGeometry/Foldable/Util.hs @@ -61,6 +61,9 @@ instance Monoid c => HasFromFoldable (Const c) where instance HasFromFoldable1 NonEmpty where fromNonEmpty = id +instance HasFromFoldable1 Seq.Seq where + fromNonEmpty = Seq.fromList . F.toList + instance HasFromFoldable Vector.Vector where fromList = Vector.fromList diff --git a/hgeometry-examples/geojson/Main.hs b/hgeometry-examples/geojson/Main.hs new file mode 100644 index 000000000..ab29fdb75 --- /dev/null +++ b/hgeometry-examples/geojson/Main.hs @@ -0,0 +1,43 @@ +{-# LANGUAGE QuasiQuotes #-} +module Main(main) where + +import Control.Lens +import Data.Aeson +import Data.Geospatial +import Data.Maybe +import HGeometry.GeoJSON +import HGeometry.Point +import HGeometry.Polygon.Class +import HGeometry.Polygon.Simple +-- import HGeometry.Vector () +import Ipe +import Paths_hgeometry_examples +import qualified System.File.OsPath as File +import System.OsPath +-------------------------------------------------------------------------------- + +type R = Double + +parseGeoJSONFile :: OsPath -> IO (Either String (GeoFeatureCollection Value)) +parseGeoJSONFile = fmap eitherDecode . File.readFile + + +-- toPolygon :: IpeOut GeoPolygon Ipe.Path R +-- toPolygon = ipePolygon . fromMaybe (error "failed") . toPolygon' + +-- toPolygon' :: GeoPolygon -> Maybe (SimplePolygon (Point 2 R)) +-- toPolygon' = fromPoints . view (vertices.asPoint) + + +main :: IO () +main = do + res <- parseGeoJSONFile [osp|data/ne_110m_admin_1_states_provinces_shp.geojson|] + case res of + Left err -> print err + Right fCollection -> do + mapM_ print $ fCollection^..geofeatures.traverse.geometry._Polygon + -- let outFp = [osp|foo.ipe|] + -- out = [ iO $ toPolygon pg + -- | pg <- fCollection^..geofeatures.traverse.geometry._Polygon + -- ] + -- writeIpeFile outFp . singlePageFromContent $ out diff --git a/hgeometry-examples/hgeometry-examples.cabal b/hgeometry-examples/hgeometry-examples.cabal index 933f756e2..07d67e006 100644 --- a/hgeometry-examples/hgeometry-examples.cabal +++ b/hgeometry-examples/hgeometry-examples.cabal @@ -19,6 +19,7 @@ extra-source-files: data-files: data/**/*.in data/**/*.out + data/**/**.geojson tested-with: GHC == 9.6.5 @@ -48,6 +49,7 @@ common setup , hgeometry:ipe , hgeometry:svg , hgeometry:miso + , hgeometry:geojson , hiraffe >= 0.1 , containers >= 0.6 , vector >= 0.13 @@ -71,6 +73,9 @@ common setup , transformers >= 0.6.0.0 , infinite-list >= 0.1.1 && < 0.2 , text >= 2.0 + , file-io + , filepath + , geojson >= 4.1.1 -- , dependent-map >= 0.4 -- , dependent-sum >= 0.7.1 @@ -188,6 +193,17 @@ executable hgeometry-polyLineDrawing -- other-modules: -- Miso.Event.Extra +-------------------------------------------------------------------------------- +-- * GeoJSON example + +executable hgeometry-geojson + import: setup, miso-setup + hs-source-dirs: geojson + main-is: Main.hs + other-modules: + Paths_hgeometry_examples + -- Miso.Event.Extra + -------------------------------------------------------------------------------- -- * Polyline Drawing diff --git a/hgeometry-examples/skia/Main.hs b/hgeometry-examples/skia/Main.hs index 786da021d..6f7a572c3 100644 --- a/hgeometry-examples/skia/Main.hs +++ b/hgeometry-examples/skia/Main.hs @@ -52,9 +52,10 @@ import SkiaCanvas ( mouseCoordinates, dimensions, canvasKitRefs, surfa , Canvas ) -import SkiaCanvas.CanvasKit.Image +import HGeometry.GeoJSON import SkiaCanvas.CanvasKit hiding (Style(..)) import SkiaCanvas.CanvasKit.GeomPrims (ltrbRect) +import SkiaCanvas.CanvasKit.Image import SkiaCanvas.CanvasKit.Paint (SkPaintRef) import SkiaCanvas.CanvasKit.Picture (serialize, withPicture, drawPicture) import SkiaCanvas.CanvasKit.PictureRecorder (recordAsPicture) diff --git a/hgeometry/geojson/src/HGeometry/GeoJSON.hs b/hgeometry/geojson/src/HGeometry/GeoJSON.hs index 6fb620e02..e37edc958 100644 --- a/hgeometry/geojson/src/HGeometry/GeoJSON.hs +++ b/hgeometry/geojson/src/HGeometry/GeoJSON.hs @@ -1,10 +1,335 @@ +{-# OPTIONS_GHC -Wno-orphans #-} module HGeometry.GeoJSON - ( - + ( GeoPositionWithoutCRS' + , RestGeoPosition(..) + , _GeoPositionWithoutCRS ) where +import Control.Lens +import Data.Bifunctor +import Data.Coerce +import qualified Data.Foldable as F +import Data.Foldable1 +import Data.Functor.Apply (Apply, (<.*>), MaybeApply(..)) +import Data.Geospatial +import Data.LinearRing +import Data.List.NonEmpty (NonEmpty(..)) +import Data.Maybe (fromMaybe) +import Data.Semigroup.Traversable +import Data.Sequence (Seq(..)) +import qualified Data.Sequence as Seq +import GHC.Generics (Generic) +import HGeometry.Cyclic +import HGeometry.Ext +import HGeometry.Foldable.Util +import HGeometry.Point.Class +import HGeometry.Polygon.Class +import HGeometry.Polygon.Simple +import HGeometry.Properties +import HGeometry.Vector + +-------------------------------------------------------------------------------- +-- * PointXY + +type instance NumType PointXY = Double +type instance Dimension PointXY = 2 + +instance HasVector PointXY PointXY where + vector = lens (\(PointXY x y) -> Vector2 x y) (\_ (Vector2 x y) -> PointXY x y) + +instance HasCoordinates PointXY PointXY +instance Affine_ PointXY 2 Double +instance Point_ PointXY 2 Double + +-------------------------------------------------------------------------------- +-- * PointXYZ + +type instance NumType PointXYZ = Double +type instance Dimension PointXYZ = 3 + +instance HasVector PointXYZ PointXYZ where + vector = lens (\(PointXYZ x y z) -> Vector3 x y z) (\_ (Vector3 x y z) -> PointXYZ x y z) + +instance HasCoordinates PointXYZ PointXYZ +instance Affine_ PointXYZ 3 Double +instance Point_ PointXYZ 3 Double + +-------------------------------------------------------------------------------- +-- * PointXYZM + +type instance NumType PointXYZM = Double +type instance Dimension PointXYZM = 4 + +instance HasVector PointXYZM PointXYZM where + vector = lens (\(PointXYZM x y z m) -> Vector4 x y z m) + (\_ (Vector4 x y z m) -> PointXYZM x y z m) + +instance HasCoordinates PointXYZM PointXYZM +instance Affine_ PointXYZM 4 Double +instance Point_ PointXYZM 4 Double + +-------------------------------------------------------------------------------- +-- * GeoPositionWithoutCRS + +type GeoPositionWithoutCRS' = PointXY :+ Maybe RestGeoPosition + +type instance NumType GeoPositionWithoutCRS = Double + +data RestGeoPosition = Z {-#UNPACK#-}!Double + | ZM {-#UNPACK#-}!Double {-#UNPACK#-}!Double + deriving (Show,Eq,Ord) + +-- | Convert between a GeoPositionWithoutCRS and a PointXY +_GeoPositionWithoutCRS :: Prism' GeoPositionWithoutCRS (PointXY :+ Maybe RestGeoPosition) +_GeoPositionWithoutCRS = prism' toGeoP fromGeoP + where + toGeoP (p@(PointXY x y) :+ mr) = case mr of + Nothing -> GeoPointXY p + Just (Z z) -> GeoPointXYZ $ PointXYZ x y z + Just (ZM z m) -> GeoPointXYZM $ PointXYZM x y z m + fromGeoP = \case + GeoEmpty -> Nothing + GeoPointXY p -> Just $ p :+ Nothing + GeoPointXYZ (PointXYZ x y z) -> Just $ (PointXY x y) :+ Just (Z z) + GeoPointXYZM (PointXYZM x y z m) -> Just $ (PointXY x y) :+ Just (ZM z m) + + +-------------------------------------------------------------------------------- +-- * Polygon + +type instance NumType GeoPolygon = Double + +type instance Dimension GeoPolygon = 2 + +instance Wrapped GeoPolygon +instance Rewrapped GeoPolygon GeoPolygon + +-------------------------------------------------------------------------------- + +newtype Seq1 a = Seq1 (Seq a) + deriving newtype (Show,Eq,Ord,Functor,Foldable) + deriving Generic + +_Seq1Seq :: Iso (Seq1 a) (Seq1 b) (Seq a) (Seq b) +_Seq1Seq = coerced + +instance Traversable Seq1 where + traverse f (Seq1 s) = Seq1 <$> traverse f s + +type instance Index (Seq1 a) = Int -- Index (Seq a) +type instance IxValue (Seq1 a) = a -- IxValue (Seq a) + +instance Ixed (Seq1 a) where + ix i = _Seq1Seq . ix i + +instance Wrapped (Seq1 a) -- where + -- type Unwrapped (Seq1 a) = Seq a + -- _Wrapped' = _Seq1Seq + +instance Rewrapped (Seq1 a) (Seq1 b) + +instance FunctorWithIndex Int Seq1 where + imap f (Seq1 s) = Seq1 $ imap f s + +instance FoldableWithIndex Int Seq1 where + ifoldMap f (Seq1 s) = ifoldMap f s + +instance TraversableWithIndex Int Seq1 where + itraverse f (Seq1 s) = Seq1 <$> itraverse f s + +instance Foldable1 Seq1 where + foldMap1 = foldMap1Default +instance Traversable1 Seq1 where + traverse1 f (Seq1 s) = Seq1 <$> traverse1Seq f s + +instance HasFromFoldable1 Seq1 where + fromNonEmpty = Seq1 . fromNonEmpty + +-- | Traverse a non-empty sequence +traverse1Seq :: Apply f => (a -> f b) -> Seq a -> f (Seq b) +traverse1Seq f = \case + Seq.Empty -> error "traverse1Seq: precondition violated" + (x :<| s) -> (:<|) <$> f x <.*> traverse1Maybe f s + +-- | Indexed traversal of a non-empty Sequence +traversed1Seq :: IndexedTraversal1 Int (Seq a) (Seq b) a b +traversed1Seq = conjoined traverse1Seq (indexing traverse1Seq) + + + + + + + +-- type instance Index (LinearRing a) = Int +-- type instance IxValue (LinearRing a) = a + +-------------------------------------------------------------------------------- + +-- | pre: the sequence has at leat 3 elements +_RingSeq :: (Eq b, Show b) => Iso (LinearRing a) (LinearRing b) (Seq a) (Seq b) +_RingSeq = iso ringToSeq (fromMaybe err . seqToRing) + where + err = error "_RingSeq: failed" + +ringToSeq :: LinearRing a -> Seq a +ringToSeq r = let (s :|> _) = toSeq r in s + +seqToRing :: (Eq b, Show b) => Seq b -> Maybe (LinearRing b) +seqToRing s@(x :<| _) = either (const Nothing) Just . fromSeq $ s |> x + +-- | pre: the sequence has at least three elements +_RingSeq1 :: (Eq b, Show b) => Iso (LinearRing a) (LinearRing b) (Seq1 a) (Seq1 b) +_RingSeq1 = _RingSeq . from _Seq1Seq + + +-- traverse1Ring :: Apply f => (a -> f b) -> LinearRing a -> f (LinearRing b) +-- traverse1Ring f + +-- | convert a ring into a polygon +-- +-- pre: the ring should be in CCW order. +-- +-- O(1) +uncheckedRingAsPolygon :: LinearRing GeoPositionWithoutCRS + -> SimplePolygonF Seq GeoPositionWithoutCRS +uncheckedRingAsPolygon = MkSimplePolygon . ringToSeq + +-- uncheckedRingAsPolygon' :: LinearRing GeoPositionWithoutCRS +-- -> SimplePolygonF Seq GeoPositionWithoutCRS' +-- uncheckedRingAsPolygon' = uncheckedRingAsPolygon + +-- MkSimplePolygon . ringToSeq + +_RingAsSimplePolygon :: Iso' (LinearRing GeoPositionWithoutCRS) + (SimplePolygonF Seq1 GeoPositionWithoutCRS) +_RingAsSimplePolygon = _RingSeq1 . from _SimplePolygonContainer + +_SimplePolygonContainer :: Iso (SimplePolygonF Seq1 point) (SimplePolygonF Seq1 point') + (Seq1 point) (Seq1 point') +_SimplePolygonContainer = coerced + +_RingAsSimplePolygon' :: Iso' (LinearRing GeoPositionWithoutCRS) + (SimplePolygonF Seq1 GeoPositionWithoutCRS') +_RingAsSimplePolygon' = _RingSeq1 . go . from _SimplePolygonContainer + where + go = iso (fmap (view $ singular _GeoPositionWithoutCRS)) + (fmap (view $ singular (re _GeoPositionWithoutCRS))) + +-- instance HasDirectedTraversals Seq where + +instance HasDirectedTraversals Seq1 where + traverseRightFrom i = conjoined traverse' (itraverse' . indexed) + where + combine sb' sa' = Seq1 (sa' <> sb') + + traverse' :: Apply f => (a -> f b) -> Seq1 a -> f (Seq1 b) + traverse' f (Seq1 s) = let (sa,sb) = Seq.splitAt i s + in combine <$> traverse1Seq f sb + <.*> traverse1Maybe f sa + itraverse' :: Apply f => (Int -> a -> f b) -> Seq1 a -> f (Seq1 b) + itraverse' f (Seq1 s) = let (sa,sb) = Seq.splitAt i s + in combine <$> itraverse1Seq (\j x -> f (i+j) x) sb + <.*> itraverseMaybe f sa + + traverseLeftFrom i = conjoined traverse' (itraverse' . indexed) + where + combine sb' sa' = Seq1 (Seq.reverse sa' <> Seq.reverse sb') + + traverse' :: Apply f => (a -> f b) -> Seq1 a -> f (Seq1 b) + traverse' f (Seq1 s) = let (sa,sb) = Seq.splitAt i s + in combine <$> traverse1Seq f (Seq.reverse sa) + <.*> traverse1Maybe f (Seq.reverse sb) + itraverse' :: Apply f => (Int -> a -> f b) -> Seq1 a -> f (Seq1 b) + itraverse' f (Seq1 s) = let (sa,sb) = Seq.splitAt (i+1) s + in combine <$> itraverse1Seq (\j x -> f (i+j) x) (Seq.reverse sa) + <.*> itraverseMaybe f (Seq.reverse sb) + + + +itraverse1Seq :: Apply f => (Int -> a -> f b) -> Seq a -> f (Seq b) +itraverse1Seq = itraverseOf traversed1Seq + +itraverseMaybe :: (Apply f, TraversableWithIndex i t) + => (i -> a -> f b) -> t a -> MaybeApply f (t b) +itraverseMaybe f = itraverse (\i -> MaybeApply . Left . f i) + + +-------------------------------------------------------------------------------- + +instance HasVertices' GeoPolygon where + type VertexIx GeoPolygon = (Int,Int) -- first one is the ringId + type Vertex GeoPolygon = GeoPositionWithoutCRS' + + vertexAt (i,j) = _Wrapped .> iix i <.> _RingAsSimplePolygon .> vertexAt j + <. _GeoPositionWithoutCRS + + numVertices (GeoPolygon ss) = F.foldl' (\a s -> a + ringLength s - 1) 0 ss + -- -- the ring incluces a copy of the first element, so it overestimates the length by 1 + +instance HasVertices GeoPolygon GeoPolygon where + vertices = _Wrapped .> traversed1Seq <.> _RingAsSimplePolygon .> vertices + <. singular _GeoPositionWithoutCRS +-- TODO: the internal ones should be reversed + +withRing :: Int -> IndexedLens' Int GeoPolygon (SimplePolygonF Seq1 GeoPositionWithoutCRS) +withRing i = _Wrapped .> singular (iix i) <. _RingAsSimplePolygon + +firstRing :: IndexedLens' Int GeoPolygon (SimplePolygonF Seq1 GeoPositionWithoutCRS) +firstRing = withRing 0 + +-- | unsafe conversion wrapping an edge +edge' :: Iso' (GeoPositionWithoutCRS, GeoPositionWithoutCRS) + (GeoPositionWithoutCRS', GeoPositionWithoutCRS') +edge' = iso (\(u,v) -> (u^?!_GeoPositionWithoutCRS, v^?!_GeoPositionWithoutCRS)) + (\(u,v) -> (u^?!re _GeoPositionWithoutCRS, v^?!re _GeoPositionWithoutCRS)) + +reIndexEdge :: Indexable ( (Int,Int) , (Int,Int)) p + => (Indexed (Int,(Int,Int)) a b -> r) -> p a b -> r +reIndexEdge = reindexed (\(i,(u,v)) -> ((i,u), (i,v))) + +instance HasOuterBoundary GeoPolygon where + outerBoundaryVertexAt (i,j) + | i == 0 = firstRing <.> singular (vertexAt j) <. singular _GeoPositionWithoutCRS + | otherwise = error "outerBoundaryVertex: not on first ring" + -- the first ring is outer boundary apparently + + ccwOuterBoundaryFrom (i,j) + | i == 0 = firstRing <.> ccwOuterBoundaryFrom j <. singular _GeoPositionWithoutCRS + | otherwise = error "ccwOuterBoundaryFrom: not on first ring" + + outerBoundary = firstRing <.> vertices <. singular _GeoPositionWithoutCRS + + cwOuterBoundaryFrom (i,j) + | i == 0 = firstRing <.> cwOuterBoundaryFrom j <. singular _GeoPositionWithoutCRS + | otherwise = error "cwOuterBoundaryFrom: not on first ring" + + outerBoundaryEdges = reIndexEdge + $ firstRing <.> outerBoundaryEdges <. edge' + + + outerBoundaryEdgeAt (_,j) = reIndexEdge + $ firstRing <.> outerBoundaryEdgeAt j <. edge' + +-- instance Polygon_ GeoPolygon GeoPositionWithoutCRS' Double where +-- area pg = case toNonEmptyOf (_Wrapped.from _Seq1Seq.traverse1._RingAsSimplePolygon') pg of +-- outer :| inners -> area outer - sum (map area inners) + + -- ccwPredecessorOf (i,j) = withRing i <.> ccwPredecessorOf j <. singular _GeoPositionWithoutCRS + -- ccwSuccessorOf (i,j) = withRing i <.> ccwSuccessorOf j <. singular _GeoPositionWithoutCRS + +{- + +instance SimplePolygon_ GeoPolygon GeoPositionWithoutCRS' Double where + uncheckedFromCCWPoints = undefined + fromPoints = undefined -- featureCollection :: -- just make instances + +-} + +-- parseGeoJSONFile :: OsPath -> IO () +-- parseGeoJSONFile path = diff --git a/hgeometry/hgeometry.cabal b/hgeometry/hgeometry.cabal index d645a5c79..cfb850b6a 100644 --- a/hgeometry/hgeometry.cabal +++ b/hgeometry/hgeometry.cabal @@ -535,7 +535,12 @@ library geojson visibility: public hs-source-dirs: geojson/src build-depends: - geojson >= 4.1.1 && < 5.0 + geojson >= 4.1.1 && < 5.0 + -- , hgeometry-combinatorial + , hgeometry:vector + , hgeometry:point + , hgeometry:kernel + , hgeometry exposed-modules: HGeometry.GeoJSON