hgis 1.0.0.2 → 1.0.0.3

raw patch · 10 files changed

+41/−234 lines, 10 filesdep +sphericaldep ~Chartdep ~Chart-cairodep ~Chart-diagramssetup-changedPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: spherical

Dependency ranges changed: Chart, Chart-cairo, Chart-diagrams, ansi-wl-pprint, binary, bytestring, colour, composition-prelude, data-binary-ieee754, data-default, directory, filepath, hspec, lens, monad-loops

API changes (from Hackage documentation)

Files

• CHANGELOG.md +5/−0
• README.md +2/−39
• Setup.hs +0/−2
• hgis.cabal +6/−4
• src/GIS/Graphics/Plot.hs +9/−2
• src/GIS/Hylo.hs +7/−5
• src/GIS/Math/Projections.hs +2/−77
• src/GIS/Math/Spherical.hs +3/−69
• src/GIS/Math/Utils.hs +0/−29
• test/Spec.hs +7/−7

@@ -0,0 +1,5 @@+# hgis++## 1.0.0.3++  * Use `spherical` package

@@ -1,51 +1,14 @@-# HGIS (aichgis)+# HGIS  ## About -This is a library and command-line tool for GIS.+This is a library for GIS.  Currently, it generates maps and computes perimters, areas, and compactness. You can use it to look at your congressional district for evidence of gerrymandering. -## Installation--hgis is best installed with [stack](https://haskellstack.org/). Try--```-stack install hgis-```--or--```-stack install hgis --flag hgis:-cairo-```--to install with bindings [cairo](https://cairographics.org/) library, which will-enable outputs as PNGs.--You can also clone this repository and type `stack build` if you want to start hacking.- ## Usage--### Command-line--From the command line, type--```-hgis map --output OUTPUT.svg FILE.shp-```--to read a shapefile and write an svg.--You can use e.g.--```-hgis map --output OUTPUT.svg --projection bonne FILE.shp-```--to force the use of the Bonne projection.  ### Library 

@@ -1,2 +0,0 @@-import Distribution.Simple-main = defaultMain

@@ -1,6 +1,6 @@ cabal-version: 2.0 name: hgis-version: 1.0.0.2+version: 1.0.0.3 license: BSD3 license-file: LICENSE copyright: Copyright: (c) 2016-2018 Vanessa McHale@@ -16,6 +16,7 @@ extra-source-files:     cabal.project.local extra-doc-files: README.md+                 CHANGELOG.md                  docs/manual.tex                  dbf/LICENSE @@ -45,7 +46,6 @@     hs-source-dirs: src     other-modules:         GIS.Graphics.Types-        GIS.Math.Utils         GIS.Utils         GIS.Graphics.Plot         GIS.Types@@ -63,7 +63,8 @@         directory -any,         colour -any,         data-default -any,-        hgis-readshp -any+        hgis-readshp -any,+        spherical -any          if flag(cairo)         exposed-modules:@@ -113,7 +114,8 @@     build-depends:         base -any,         hgis -any,-        hspec -any+        hspec -any,+        spherical >=0.1.2.0          if flag(development)         ghc-options: -Werror

@@ -15,11 +15,18 @@ import           Data.Colour.Names import           Data.Foldable                 (fold) import           GIS.Graphics.Types            hiding (title)-import           GIS.Math.Utils-import           GIS.Types+import           GIS.Types                     as GIS import           GIS.Utils import           Graphics.Rendering.Chart import           Graphics.Rendering.Chart.Easy hiding (lens, points)++-- | Averages the coördinates of a polygon, returning a point.+shittyCentroid :: Polygon -> GIS.Point+shittyCentroid poly = (avg $ fmap fst poly, avg $ fmap snd poly)++-- | Average over a foldable container+avg :: (RealFrac a, Foldable t) => t a -> a+avg list = sum list / (fromIntegral . length $ list)  -- | Given a map, return a `Renderable ()` for use with the -- """Graphics.Rendering.Char""" module.

@@ -27,6 +27,7 @@ import           GIS.Math.Spherical import           GIS.Types                      hiding (area) import           GIS.Utils+import           Math.Geometry.Spherical import           System.Directory  -- | Get the areas of various objects and return a string suitable for printing@@ -70,19 +71,20 @@ districtToMapFiles = fmap g where     g (District polygons label _ area' _) = title .~ label <> "-" <> (show . sum $ area') $ labelledDistricts .~ pure (polygons,"") $ def --- | Given the path to a shapefile, return a list of districts-getDistricts :: FilePath -> IO [District]+-- | Given the path to a shapefile, return a list of districts, if it contains+-- them.+getDistricts :: FilePath -> IO (Maybe [District]) getDistricts filepath = do         dbfExists <- doesFileExist (stripExt filepath <> ".dbf")         file <- if dbfExists then readShpWithDbf filepath else readShpFile filepath         -- FIXME throw a real error.         let extract f = fmap (f . getPolygon . fromJust . shpRecContents) . shpRecs $ file-            districtLabels = fromJust $ fmap labels $ traverse shpRecLabel . shpRecs $ file+            districtLabels = fmap labels $ traverse shpRecLabel . shpRecs $ file             shapes = extract id             perimeters = extract totalPerimeter-            areas = extract (fmap areaPolygon)+            areas = extract (fmap (areaPolygon 6371))             compacticity = extract (relativeCompactness . fold)-        pure $ zipWith5 District shapes districtLabels perimeters areas compacticity+        pure $ zipWith5 District shapes <$> districtLabels <*> pure perimeters <*> pure areas <*> pure compacticity  -- | Helper function for extracting from shapefiles. getPolygon :: RecContents -> [Polygon]

@@ -1,87 +1,12 @@ -- | Module containing several useful projections for generating maps.-module GIS.Math.Projections ( -- * Reference points-                              washingtonDC-                            , mecca-                            -- * Projections-                            , littow-                            , craig-                            , winkel3-                            , mercator-                            , bonne-                            , albers-                            -- * Helper functions-                            , project-                            , projectMap+module GIS.Math.Projections ( projectMap                             ) where  import           Control.Arrow import           Control.Lens import           GIS.Graphics.Types-import           GIS.Math.Utils import           GIS.Types---- | For use as a reference point-washingtonDC :: Point-washingtonDC = toRadians (38.9072, -77.0369)---- | For use as a reference point-mecca :: Point-mecca = toRadians (21.3891, 39.8579)---- | Littow retroazimuthal + conformal projection-littow :: Projection-littow (long, lat) = (sin(long - referenceLong)/cos lat, cos (long - referenceLong) * tan lat)-    where referenceLong = radians 0---- | Craig retroazimuthal projection-craig :: Point -> Projection-craig referencePoint (long, lat) = (long - referenceLong, y)-    where (referenceLong, referenceLat) = referencePoint-          expr = sin lat * cos (long - referenceLong) - tan referenceLat * cos lat-          y | long - referenceLong == 0 = expr-            | otherwise = (long - referenceLong) / sin (long - referenceLong) * expr---- | Winkel Tripel projection-winkel3 :: Projection-winkel3 (long, lat) = ((lambda * cos phi1 + (2 * cos lat * sin (lambda/2)/sinc alpha))/2, (lat + sin lat/sinc alpha)/2)-    where lambda = long - lambda0-          phi1 = acos $ 2 / pi-          alpha = acos $ cos lat * cos (lambda/2)-          lambda0 = radians (-77.0369)---- | Mercator projection.-mercator :: Projection-mercator (long, lat) = (long - meridian, asinh (tan lat))-    where meridian = radians (-98.5795)---- | Bonne projection with standard parallel at 45 N and central meridian--- centered at Washington DC.-bonne :: Projection-bonne (long, lat) = (rho * sin e, cot phi1 - rho * cos e)-    where rho = cot phi1 + phi1 - lat-          e = (long - meridian) * cos lat / rho-          phi1 = radians 45 -- standard parallel @ 45 N-          meridian = radians (-77.0369) -- central meridian @ dc-          cot = (1/) . tan---- | Albers projection for a given reference point.------ > ablers washingtonDC-albers :: Point -> Projection-albers referencePoint (long, lat) = (rho * sin theta, rho' - rho * cos theta)-    where n = (sin phi1 + sin phi2) / 2-          theta = n * (long - referenceLong)-          c = cos phi1^(2 :: Int) + 2 * n * sin phi1-          rho = sqrt (c - 2 * n * sin lat) / n-          rho' = sqrt (c - 2 * n * sin referenceLat) / n-          -- standard parallels @ 20, 50 degrees-          phi1 = radians 20-          phi2 = radians 50-          (referenceLong, referenceLat) = referencePoint---- | Project given a `Polygon`.-project :: Projection -> Polygon -> Polygon-project f = fmap (f . toRadians)+import           Math.Geometry.Spherical  -- | Apply a given projection a `Map`. projectMap :: Projection -> Map -> Map

@@ -1,76 +1,10 @@ -- | Utilities to compute area, perimeterPolygon, etc. on the surface of a sphere.-module GIS.Math.Spherical ( areaPolygon-                          , perimeterPolygon-                          , relativeCompactness-                          , totalPerimeter-                          -- * Spherical geometry-                          , areaTriangle-                          , areaConvex-                          -- * Internal functions-                          , distance+module GIS.Math.Spherical ( totalPerimeter                           ) where -import           Control.Composition-import           GIS.Math.Projections-import           GIS.Math.Utils import           GIS.Types---- | Compute the area of a triangle using L'Huillier's formula-areaTriangle :: Point -> Point -> Point -> Double-areaTriangle x1 x2 x3 = r^(2 :: Int) * e-    where r = 6371-          e = 4 * atan(sqrt(tan(s/2) * tan((s - a)/2) * tan((s - b)/2) * tan((s - c)/2)))-          s = (a + b + c) / 2-          a = distanceRad x1 x2-          b = distanceRad x1 x3-          c = distanceRad x2 x3-          distanceRad = on centralAngle toRadians---- TODO mandelbrot/fractal dimension?--- consider "area of largest circumscribable circle" as well.---- | Relative compactness. Dimensionless.-relativeCompactness :: Polygon -> Double-relativeCompactness = (*scale) . compactness1-    where scale = 1/4*pi---- | Take the area of the polygon and divide by the perimeter squared. Dimensionless.-compactness1 :: Polygon -> Double-compactness1 p = areaPolygon p/perimeterPolygon p^(2 :: Int)---- | Compute the area of a convex polygon on the surface of a sphere.-areaConvex :: Polygon -> Double-areaConvex (base1:base2:pts) = fst $ foldr stepArea (0,base2) pts-    where stepArea point (sum', base) = (sum' + areaTriangle base1 base point, point)-areaConvex _ = error "attempted to take area of polygon with < 3 points"---- | Uses areal projection; then finds area of the polygon.--- Result is in km^2-areaPolygon :: Polygon -> Double-areaPolygon = (*factor) . areaPolyRectangular . fmap (bonne . toRadians)-    where factor = 1717856/4.219690791828533e-2+import           Math.Geometry.Spherical  -- | Given a list of polygons, return the total perimeter. totalPerimeter :: [Polygon] -> Double-totalPerimeter = sum . fmap perimeterPolygon--perimeterPolygon :: Polygon -> Double-perimeterPolygon [x1, x2]       = distance x1 x2-perimeterPolygon (x1:x2:points) = perimeterPolygon (x2:points) + distance x1 x2-perimeterPolygon _              = error "Attempted to take area of polygon with no points"---- | Find the area of a polygon with rectangular coördinates given. This--- function will error if given a polygon with no points.-areaPolyRectangular :: Polygon -> Double-areaPolyRectangular (pt:pts) = abs . (*0.5) . fst $ foldr areaPolyCalc (0,pt) pts-    where areaPolyCalc (x2, y2) (sum',(x1,y1)) = (sum' + (x1 * y2 - x2 * y1),(x2,y2))-areaPolyRectangular _ = error "Attempted to take area of polygon with no points"---- | Distance in kilometers between two points given in degrees.-distance :: Point -> Point -> Double-distance = (*6371) .* on centralAngle toRadians---- | Compute central angle from points given in radians-centralAngle :: Point -> Point -> Double-centralAngle (long1, lat1) (long2, lat2) =-    acos $ sin lat1 * sin lat2 + cos lat1 * cos lat2 * cos (long1 - long2)+totalPerimeter = sum . fmap (perimeterPolygon 6371)

@@ -1,29 +0,0 @@--- | Miscellaneous utils for math-module GIS.Math.Utils ( radians-                      , toRadians-                      , sinc-                      , avg-                      , shittyCentroid-                      ) where--import           Control.Composition-import           GIS.Types---- | Averages the coördinates of a polygon, returning a point.-shittyCentroid :: Polygon -> Point-shittyCentroid poly = (avg $ fmap fst poly, avg $ fmap snd poly)---- | Average over a foldable container-avg :: (RealFrac a, Foldable t) => t a -> a-avg list = sum list / (fromIntegral . length $ list)--sinc :: Floating a => a -> a-sinc x = sin x / x---- | Convert a `Double` from degrees to radians.-radians :: Double -> Double-radians = (*(pi/180))---- | Convert both coördinates to radians.-toRadians :: Point -> Point-toRadians = both radians

@@ -1,20 +1,20 @@ module Main (main) where +import           Data.Maybe import           GIS.Graphics.PlotPNG import           GIS.Graphics.PlotSVG import           GIS.Hylo-import           GIS.Math.Projections-import           GIS.Math.Spherical+import           Math.Geometry.Spherical import           Test.Hspec  main :: IO () main = hspec $ do     describe "distance" $         it "computes distances along geodesics" $-            distance (-74.0059, 40.7128) (-118.2347, 34.0522) `shouldBe` 3934.997673163554+            distance (6371 :: Double) (-74.0059, 40.7128) (-118.2347, 34.0522) `shouldBe` 3934.997673163554     describe "mercatorFullPng" $-        it "Reads a map and writes a .png with the mercator projection" $-            (mkMapPng "test/testmap.png" =<< districtToMapP bonne <$> getDistricts "test/data/worldmap/TM_WORLD_BORDERS-0.3.shp") >>= (`shouldBe` ())+        it "Reads a map and writes a .png with the Bonne projection" $+            (mkMapPng "test/testmap.png" =<< districtToMapP (bonne (radians 45) (snd washingtonDC)) . fromJust <$> getDistricts "test/data/worldmap/TM_WORLD_BORDERS-0.3.shp") >>= (`shouldBe` ())     describe "mercatorFullSVG" $-        it "Reads a map and writes a .svg with the mercator projection" $-            (mkMapSVG "test/testmap.svg" =<< districtToMapP bonne <$> getDistricts "test/data/worldmap/TM_WORLD_BORDERS-0.3.shp") >>= (`shouldBe` ())+        it "Reads a map and writes a .svg with the Bonne projection" $+            (mkMapSVG "test/testmap.svg" =<< districtToMapP (bonne (radians 45) (snd washingtonDC)) . fromJust <$> getDistricts "test/data/worldmap/TM_WORLD_BORDERS-0.3.shp") >>= (`shouldBe` ())