diff --git a/hgis.cabal b/hgis.cabal --- a/hgis.cabal +++ b/hgis.cabal @@ -1,100 +1,88 @@ -name: hgis -version: 0.1.3.2 -cabal-version: >=1.10 -build-type: Simple -license: BSD3 -license-file: LICENSE -copyright: Copyright: (c) 2016 Vanessa McHale -maintainer: tmchale@wisc.edu -homepage: https://github.com/vmchale/hgis#readme -synopsis: Package and command-line for GIS with Haskell -description: - Package containting functions to make graphs, read - shapefiles, and compute areas/perimeters of - geographic features. -category: Web -author: Vanessa McHale -extra-source-files: - README.md - stack.yaml - docs/manual.tex - docs/manual.pdf - src/depends/readshp/LICENSE - LICENSE-readshp - -source-repository head - type: git - location: https://github.com/vmchale/hgis +name: hgis +version: 0.1.3.3 +synopsis: Package and command-line for GIS with Haskell +description: Package containting functions to make graphs, read + shapefiles, and compute areas/perimeters of + geographic features. +homepage: https://github.com/vmchale/hgis#readme +license: BSD3 +license-file: LICENSE +author: Vanessa McHale +maintainer: tmchale@wisc.edu +copyright: Copyright: (c) 2016 Vanessa McHale +category: Web +build-type: Simple +extra-source-files: README.md, stack.yaml, docs/manual.tex, docs/manual.pdf, src/depends/readshp/LICENSE, LICENSE-readshp +cabal-version: >=1.10 -flag llvm-fast - description: - Enable build with llvm backend - default: False +Flag llvm-fast { + Description: Enable build with llvm backend + Default: False +} library - exposed-modules: - GIS.Exe.Opt - GIS.Exe.OptCairo - GIS.Math.Spherical - GIS.Graphics.PlotPNG - GIS.Graphics.PlotSVG - GIS.Hylo - GIS.Math.Projections - build-depends: - base >=4.7 && <5, - optparse-applicative >=0.13.2.0 && <0.14, - Chart >=1.8.2 && <1.9, - Chart-cairo >=1.8.2 && <1.9, - Chart-diagrams >=1.8.2 && <1.9, - lens >=4.15.1 && <4.16, - composition >=1.0.2.1 && <1.1, - ansi-wl-pprint >=0.6.7.3 && <0.7, - transformers >=0.5.2.0 && <0.6, - directory >=1.3.0.0 && <1.4, - colour >=2.3.3 && <2.4, - data-default >=0.7.1.1 && <0.8, - binary >=0.8.3.0 && <0.9, - bytestring >=0.10.8.1 && <0.11, - data-binary-ieee754 >=0.4.4 && <0.5, - filepath >=1.4.1.1 && <1.5, - monad-loops >=0.4.3 && <0.5 - default-language: Haskell2010 - default-extensions: DeriveGeneric OverloadedStrings - hs-source-dirs: src src/depends/readshp - other-modules: - GIS.Math.Utils - GIS.Utils - GIS.Graphics.Plot - GIS.Graphics.Types - GIS.Types - GIS.Exe.Parser - Geometry.Shapefile - Geometry.Shapefile.MergeShpDbf - Geometry.Shapefile.ReadDbf - Geometry.Shapefile.ReadShp - Geometry.Shapefile.Types - Geometry.Shapefile.Internal + hs-source-dirs: src, src/depends/readshp + exposed-modules: GIS.Exe.Opt + , GIS.Exe.OptCairo + , GIS.Math.Spherical + , GIS.Graphics.PlotPNG + , GIS.Graphics.PlotSVG + , GIS.Hylo + , GIS.Math.Projections + other-modules: GIS.Math.Utils + , GIS.Utils + , GIS.Graphics.Plot + , GIS.Graphics.Types + , GIS.Types + , GIS.Exe.Parser + , Geometry.Shapefile + , Geometry.Shapefile.MergeShpDbf + , Geometry.Shapefile.ReadDbf + , Geometry.Shapefile.ReadShp + , Geometry.Shapefile.Types + , Geometry.Shapefile.Internal + build-depends: base >= 4.7 && < 5 + , optparse-applicative + , Chart + , Chart-cairo + , Chart-diagrams + , lens + , composition + , ansi-wl-pprint + , transformers + , directory + , colour + , data-default + , binary + , bytestring + , data-binary-ieee754 + , filepath + , monad-loops + default-language: Haskell2010 + default-extensions: DeriveGeneric + , OverloadedStrings executable hgis - - if flag(llvm-fast) - ghc-options: -threaded -rtsopts -with-rtsopts=-N -fllvm -optlo-O3 -O3 - else - ghc-options: -threaded -rtsopts -with-rtsopts=-N - main-is: MainPng.hs - build-depends: - base >=4.9.1.0 && <4.10, - hgis >=0.1.3.2 && <0.2 - default-language: Haskell2010 - hs-source-dirs: app + hs-source-dirs: app + main-is: MainPng.hs + if flag(llvm-fast) + ghc-options: -threaded -rtsopts -with-rtsopts=-N -fllvm -optlo-O3 -O3 + else + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: base + , hgis + default-language: Haskell2010 test-suite hgis-test - type: exitcode-stdio-1.0 - main-is: Spec.hs - build-depends: - base >=4.9.1.0 && <4.10, - hgis >=0.1.3.2 && <0.2, - hspec >=2.4.2 && <2.5 - default-language: Haskell2010 - hs-source-dirs: test - ghc-options: -threaded -rtsopts -with-rtsopts=-N + type: exitcode-stdio-1.0 + hs-source-dirs: test + main-is: Spec.hs + build-depends: base + , hgis + , hspec + ghc-options: -threaded -rtsopts -with-rtsopts=-N + default-language: Haskell2010 + +source-repository head + type: git + location: https://github.com/vmchale/hgis diff --git a/src/GIS/Hylo.hs b/src/GIS/Hylo.hs --- a/src/GIS/Hylo.hs +++ b/src/GIS/Hylo.hs @@ -22,18 +22,18 @@ -- | Get the areas of various objects and return a string suitable for printing districtArea :: [District] -> String districtArea districts = concat . intercalate (pure "\n") $ map (pure . show . distA) districts - where distA (District _ label _ area) = (label, sum area) -- TODO figure out which one is the correct one + where distA (District _ label _ area _) = (label, sum area) -- TODO figure out which one is the correct one -- should mention it's km -- | Get the perimeters of various objects and return a string suitable for printing districtPerimeter :: [District] -> String districtPerimeter districts = concat . intercalate (pure "\n") $ map (pure . show . distP) districts - where distP (District _ label perimeter _) = (label, perimeter) + where distP (District _ label perimeter _ _) = (label, perimeter) --- FIXME use the other function? put compactness in +-- | Label with relative compactness districtCompactness :: [District] -> String districtCompactness districts = concat . intercalate (pure "\n") $ map (pure . show . distC) districts - where distC (District _ label perimeter area) = (label, (sum area)/perimeter^2) + where distC (District _ label _ _ compacticity) = (label, compacticity) -- | Given a projection and lists of districts, draw a map. districtToMapP :: Projection -> [District] -> Map @@ -51,8 +51,7 @@ -- | Given a list of districts, return a list of maps. districtToMapFiles :: [District] -> [Map] -districtToMapFiles = map (\(District polygons label _ area) -> title .~ label ++ "-" ++ (show . sum $ area) $ labelledDistricts .~ (zip polygons (nullLabel polygons)) $ def) ---districtToMapFiles = map (\(District polygons label _ _) -> title .~ label $ labelledDistricts .~ (zip polygons (nullLabel polygons)) $ def) +districtToMapFiles = map (\(District polygons label _ area _) -> title .~ label ++ "-" ++ (show . sum $ area) $ labelledDistricts .~ (zip polygons (nullLabel polygons)) $ def) where nullLabel polys = map (const "") [1..(length polys)] -- | Given the path to a shapefile, return a list of districts @@ -64,8 +63,10 @@ let shapes = (map (getPolygon . fromJust . shpRecContents)) . shpRecs $ file let perimeters = map (totalPerimeter . getPolygon . fromJust . shpRecContents) . shpRecs $ file let areas = map (fmap areaPolygon . getPolygon . fromJust . shpRecContents) . shpRecs $ file - let tuple = zip4 shapes districtLabels perimeters areas - pure $ fmap (\(a,b,c,d) -> District a b c d) tuple + let compacticity = map (relativeCompactness . concat . getPolygon . fromJust . shpRecContents) . shpRecs $ file + pure $ zipWith5 (\a b c d e -> District a b c d e) shapes districtLabels perimeters areas compacticity + --let tuple = zip5 shapes districtLabels perimeters areas compactness + --pure $ fmap (\(a,b,c,d,e) -> District a b c d e) tuple -- | Helper function for extracting from shapefiles. getPolygon :: RecContents -> [Polygon] diff --git a/src/GIS/Math/Spherical.hs b/src/GIS/Math/Spherical.hs --- a/src/GIS/Math/Spherical.hs +++ b/src/GIS/Math/Spherical.hs @@ -19,6 +19,7 @@ avg list = sum list / (fromIntegral . length $ list) -- | Compute the area of a triangle using L'Huillier's formula +areaTriangle :: Point -> Point -> Point -> Double areaTriangle x1 x2 x3 = r^2 * e where r = 6371 e = 4 * atan(sqrt(tan(s/2) * tan((s - a)/2) * tan((s - b)/2) * tan((s - c)/2))) @@ -28,20 +29,25 @@ c = distanceRad x2 x3 distanceRad = on centralAngle toRadians --- mandelbrot/fractal dimension? +-- TODO mandelbrot/fractal dimension? -- consider "area of largest circumscribable circle" as well. -compactness1 :: Polygon -> Double -compactness1 p = (areaPolygon p)/(perimeterPolygon p) +-- | Relative compactness, i.e. compactness divided by the compactness of a Euclidean circle +relativeCompactness :: Polygon -> Double +relativeCompactness = (*scale) . compactness1 + where scale = (1/4*pi) +-- | Take the area of the polygon and divide by the perimeter squared. This is a dimensionless measurement. +compactness1 :: Polygon -> Double +compactness1 p = (areaPolygon p)/(perimeterPolygon p)^2 + -- | Compute the area of a convex polygon on the surface of a sphere. areaConvex :: Polygon -> Double areaConvex (base1:base2:pts) = (view _1) $ foldr stepArea (0,base2) pts where stepArea point (sum, base) = (sum + (areaTriangle base1 base point), point) -- | Uses areal projection; then finds area of the polygon. --- WARNING: not yet to scale. --- (need to compute jacobian; ideally km^2/steridian +-- Result is in km^2 areaPolygon :: Polygon -> Double areaPolygon = (*factor) . areaPolyRectangular . fmap (bonne . toRadians) where factor = 1717856/4.219690791828533e-2 diff --git a/src/GIS/Types.hs b/src/GIS/Types.hs --- a/src/GIS/Types.hs +++ b/src/GIS/Types.hs @@ -18,6 +18,7 @@ , _districtLabel :: String , _perimeter :: Double , _area :: [Double] + , _compactness :: Double } deriving (Generic, Show) makeLenses ''District