packages feed

hgis (empty) → 0.1.0.0

raw patch · 19 files changed

+759/−0 lines, 19 filesdep +Chartdep +Chart-cairodep +Chart-diagramssetup-changed

Dependencies added: Chart, Chart-cairo, Chart-diagrams, ansi-wl-pprint, base, colour, composition, data-default, directory, hgis, hspec, lens, optparse-applicative, readshp, transformers

Files

+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Vanessa McHale (c) 2017++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Vanessa McHale nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,56 @@+# HGIS (aichgis)++## About++This is a library and command-line tool 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++## Documentation++### Haddock++Haddock documentation can be accessed through hackage.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ app/MainPng.hs view
@@ -0,0 +1,6 @@+module Main where++import GIS.Exe.OptCairo++main :: IO ()+main = exec
+ hgis.cabal view
@@ -0,0 +1,77 @@+name: hgis+version: 0.1.0.0+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++source-repository head+    type: git+    location: https://github.com/vmchale/hgis++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,+        readshp >=0.1.0.0 && <0.2,+        optparse-applicative >=0.13.1.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+    pkgconfig-depends: cairo >=0.13.3.1 && <0.14+    default-language: Haskell2010+    default-extensions: DeriveGeneric OverloadedStrings+    hs-source-dirs: src+    other-modules:+        GIS.Math.Utils+        GIS.Utils+        GIS.Graphics.Plot+        GIS.Graphics.Types+        GIS.Types+        GIS.Exe.Parser++executable hgis+    main-is: MainPng.hs+    build-depends:+        base >=4.9.1.0 && <4.10,+        hgis >=0.1.0.0 && <0.2+    default-language: Haskell2010+    hs-source-dirs: app+    ghc-options: -threaded -rtsopts -with-rtsopts=-N++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.0.0 && <0.2,+        hspec >=2.4.2 && <2.5+    default-language: Haskell2010+    hs-source-dirs: test+    ghc-options: -threaded -rtsopts -with-rtsopts=-N
+ src/GIS/Exe/Opt.hs view
@@ -0,0 +1,31 @@+module GIS.Exe.Opt where++import GIS.Hylo+import Control.Monad+import Control.Lens hiding (argument)+import Options.Applicative+import Control.Monad.IO.Class+import Data.Monoid hiding (getAll)+import System.Directory+--+import GIS.Utils+import GIS.Math.Projections+import GIS.Math.Spherical+import GIS.Graphics.PlotSVG+import GIS.Types+import GIS.Exe.Parser+import GIS.Graphics.Types++pick :: Program -> IO ()+pick (Program (MapMaker outfile False projection) infile) = let p = pickProjection projection in+    mkMapSVG outfile =<< districtToMapP p <$> getDistricts infile -- svg now idk?+pick (Program (MapMaker outfile True projection) infile) = let p = pickProjection projection in+    makeFoldersSVG =<< districtToMapFilesP p <$> getDistricts infile++exec :: IO ()+exec = pick =<< execParser helpDisplay++makeFoldersSVG :: [Map] -> IO ()+makeFoldersSVG maps = do+    createDirectoryIfMissing False "maps"+    mapM_ (\m -> mkMapSVG ("maps/" <> (view title m) <> ".svg") m) maps
+ src/GIS/Exe/OptCairo.hs view
@@ -0,0 +1,48 @@+-- | This contains the main program executable for the program when built with+-- the Cairo libary+module GIS.Exe.OptCairo+    ( -- * Data types+    Program (..)+    -- * Functions+    , exec+    ) where++import GIS.Hylo+import Control.Lens+import Control.Monad+import Options.Applicative+import Control.Monad.IO.Class+import Data.Monoid hiding (getAll)+import System.Directory+--+import GIS.Utils+import GIS.Math.Projections+import GIS.Math.Spherical+import GIS.Graphics.PlotSVG+import GIS.Graphics.PlotPNG+import GIS.Types+import GIS.Exe.Opt hiding (exec, pick)+import GIS.Exe.Parser+import GIS.Graphics.Types++-- | Main executable; parses command line options and runs program+exec :: IO ()+exec = pick =<< execParser helpDisplay++-- | Execute the program with parsed command-line input+pick :: Program -> IO ()+pick (Program (MapMaker outfile False projection) infile) = let p = pickProjection projection in+    mkMap outfile =<< districtToMapP p <$> getDistricts infile -- svg now idk?+pick (Program (MapMaker outfile True projection) infile) = let p = pickProjection projection in+    makeFoldersPng =<< districtToMapFilesP p <$> getDistricts infile+pick (Program (Computation comp Nothing) infile) = --slightly wrong but eh. +    case comp of+        "perimeter" -> putStrLn =<< districtPerimeter <$> getDistricts infile+        "area" -> putStrLn =<< districtArea <$> getDistricts infile+        --compactness too! ++-- | Make maps as png files.+makeFoldersPng :: [Map] -> IO ()+makeFoldersPng maps = do+    createDirectoryIfMissing False "maps"+    mapM_ (\m -> mkMapPng ("maps/" <> (view title m) <> ".png") m) maps
+ src/GIS/Exe/Parser.hs view
@@ -0,0 +1,81 @@+-- | This module contains the parser for the command-line options+module GIS.Exe.Parser where++import Options.Applicative+import GIS.Types+import GIS.Math.Projections+import Control.Lens hiding (argument)+import Data.Monoid++-- | Data type for one run of the program+data Program = Program { com :: Command +                       , file :: FilePath -- "Path to dbf file"+                       } ++-- | Data type for the appropriate subcommand+data Command = Computation { computation :: String , outputM :: Maybe FilePath } +    | MapMaker { output :: FilePath , generateAll :: Bool, projection :: Maybe String }+    -- also a "map and label command?"++-- | Parses the `Program` data type+program :: Parser Program+program = Program +    <$> hsubparser+        ( command "compute" (info computationP ( progDesc "Compute perimeter, area, etc. of map" ))+        <> command "map" (info mapMaker ( progDesc "Make a map from a shapefile database." )))+    <*> ( argument str ( metavar "SHAPEFILE" <> help "Path to .shp file" ) )++-- | Parses the `Command` datatype into a Computation+computationP :: Parser Command+computationP = Computation+    <$> argument str+        ( metavar "COMPUTATION"+        <> help "which computation to perform, e.g. area, perimeter, etc." )+    <*> ( optional $ strOption  -- because it's a maybe filepath+        ( long "output"+        <> short 'o'+        <> metavar "OUTPUT"+        <> help "Where to write output" ) )++-- | Parses the `Command` datatype into a map+mapMaker :: Parser Command+mapMaker = MapMaker+    <$> strOption+        (long "output"+        <> short 'o'+        <> metavar "OUTPUT"+        <> help "Where to write the image/map" )+    <*> switch+        (long "generate-all"+        <> short 'a'+        <> help "Whether to generate a separate file for each object in the shapefile" )+    <*> ( optional $ strOption+        (long "projection"+        <> short 'p'+        <> help "Which projection to use, e.g. mercator etc.") )++-- | display help if user enters e.g.+--+-- > hgis --help wrongcommand+helpDisplay = info ( program <**> helper)+    ( fullDesc+    <> progDesc "GIS for Haskell. Can make maps, compute areas/perimeters for shapefiles."+    <> header "hgis - GIS in Haskell" )++{--+pickComputation :: String -> District -> [Double]+pickComputation str dist = case str of+    "perimeter" -> pure $ view perimeter dist+    "area" -> view area dist+--}++-- | Parse a `Maybe String` into the appropriate `Projection`, doing nothing for+-- a `Nothing`+pickProjection :: Maybe String -> Projection+pickProjection str = case str of+    Just "mercator" -> mercator+    Just "bonne" -> bonne+    Just "albers" -> albers washingtonDC+    Just "winkel3" -> winkel3+    Just "littow" -> littow+    Nothing -> id
+ src/GIS/Graphics/Plot.hs view
@@ -0,0 +1,40 @@+-- | Module containing renderable objects and plots for our maps+module GIS.Graphics.Plot where++import GIS.Types+import Graphics.Rendering.Chart+import Graphics.Rendering.Chart.Easy+import GIS.Math.Spherical+import GIS.Utils+import Data.Colour.Names+import GIS.Graphics.Types+import GIS.Math.Projections++-- | Given a map, return a `Renderable ()` for use with the+-- """Graphics.Rendering.Char""" module. +mkMapR :: Map -> Renderable ()+mkMapR (Map { _projection = p , _title = tit , _labelEntities = le , _labelledDistricts = points }) = case le of+    True -> mkRenderableLabelled tit points+    False -> mkRenderableLabelled tit $ fmap (over _2 (const "")) points++mkRenderablePlots :: String -> [Plot Double Double] -> Renderable ()+mkRenderablePlots title plots = toRenderable $+    layout_title .~ title $ layout_plots .~ plots $ def++mkRenderableLabelled :: String -> [(Polygon, String)] -> Renderable ()+mkRenderableLabelled title points = mkRenderablePlots title [ plotDataPoints (map fst points), plotLabels points ]++-- | Helper function to plot data points as appropriate for a map, i.e. using+-- contiguous lines. +plotDataPoints :: [Polygon] -> Plot Double Double+plotDataPoints points = toPlot $+    plot_lines_values .~ points $ plot_lines_style . line_color .~ opaque blue $ plot_lines_title .~ "Border" $ def++-- | Helper function to plot labels on a map. The ordinate will be plotted at+-- the centroid of the abscissa, which may be outside the polygon if it is+-- concave. +plotLabels :: [(Polygon, String)] -> Plot Double Double+plotLabels points = toPlot texts+    where pairs = fmap (flatten . (over _1 shittyCentroid)) points+          fontStyle = font_size .~ 15 $ font_weight .~ FontWeightBold $ def+          texts = plot_annotation_values .~ pairs $ plot_annotation_style .~ fontStyle $ def
+ src/GIS/Graphics/PlotPNG.hs view
@@ -0,0 +1,47 @@+-- | Module to generate PNGs from shapefiles+module GIS.Graphics.PlotPNG where++import Graphics.Rendering.Chart+import Graphics.Rendering.Chart.Easy+import Graphics.Rendering.Chart.Backend.Cairo+import GIS.Types+import Data.Monoid+import Control.Monad+import GIS.Hylo+import GIS.Math.Projections+import Control.Lens+import Control.Lens.Tuple+--+import GIS.Math.Spherical+import GIS.Utils+import GIS.Graphics.Plot+import GIS.Graphics.Types+import GIS.Graphics.PlotSVG hiding (fileOptions)++-- | Default file options: PNG output and 1920x1080. To change the file size,+-- you can do e.g.+-- > fo_size .~ (640,480) $ fileOptions+fileOptions = def { _fo_size = (1920, 1080) , _fo_format = PNG }++-- | Given a `Map` write it to file, where the format is determined by the+-- extension.+mkMap :: FilePath -> Map -> IO ()+mkMap filepath map = case getExt filepath of+    "png" -> mkMapPng filepath map+    "svg" -> mkMapSVG filepath map++-- | Given a `Map`, write it to file. +mkMapPng :: FilePath -> Map -> IO ()+mkMapPng path map = do+    renderableToFile fileOptions path $ mkMapR map+    putStrLn ("...output written to " <> path)++{--+makeLabelledMapPng :: String -> FilePath -> [(Polygon, String)] -> IO ()+makeLabelledMapPng title filepath points = do+    renderableToFile fileOptions filepath $ mkRenderableLabelled title points+    putStrLn ("...output written to " <> filepath)++makeMapPng' :: String -> FilePath -> [Polygon] -> IO ()+makeMapPng' = (flip zip (forever $ "")) -.** makeLabelledMapPng+--}
+ src/GIS/Graphics/PlotSVG.hs view
@@ -0,0 +1,33 @@+-- | Module to generate SVGs from shapefiles.+module GIS.Graphics.PlotSVG where++import Graphics.Rendering.Chart.Easy+import Graphics.Rendering.Chart.Backend.Diagrams +import Control.Monad+import GIS.Hylo+import GIS.Types+import Data.Monoid+import GIS.Math.Projections+import GIS.Graphics.Plot+import GIS.Graphics.Types++-- | Default file options: PNG output and 1920x1080. To change the file size,+-- you can do e.g.+-- > fo_size .~ (640,480) $ fileOptions+fileOptions = def { _fo_size = (1920, 1080) , _fo_format = SVG }++-- | Given a `Map`, write it to file as an SVG.+mkMapSVG :: FilePath -> Map -> IO ()+mkMapSVG path map = do+    renderableToFile fileOptions path $ mkMapR map+    putStrLn ("...output written to " <> path)++{--+makeLabelledMapSVG :: String -> FilePath -> [(Polygon, String)] -> IO ()+makeLabelledMapSVG title filepath points = do+    renderableToFile fileOptions filepath $ mkRenderableLabelled title points+    putStrLn ("...output written to " <> filepath)++makeMapSVG :: String -> FilePath -> [Polygon] -> IO ()+makeMapSVG = (flip zip (forever "")) -.** makeLabelledMapSVG+--}
+ src/GIS/Graphics/Types.hs view
@@ -0,0 +1,22 @@+{-# LANGUAGE TemplateHaskell #-}++-- | Module for types associated with generating maps. Includes lenses. +module GIS.Graphics.Types where++import Data.Default+import GIS.Types+import Control.Lens++-- | Data type for a map+data Map = Map { _projection :: Projection+               , _title :: String+               , _labelEntities :: Bool -- whether to label districts+               , _labelledDistricts :: [(Polygon, String)] -- the data we actually want to map+               --technically just needs to be a +               --though for that matter then projection should also have a typeclass!+               }++makeLenses ''Map++instance Default Map where+    def = Map { _projection = id , _title = mempty , _labelEntities = False , _labelledDistricts = mempty }
+ src/GIS/Hylo.hs view
@@ -0,0 +1,73 @@+-- | Module to transform shapefiles.+module GIS.Hylo where++import Geometry.Shapefile.MergeShpDbf+import Geometry.Shapefile.ReadShp+import Geometry.Shapefile.Types+import Data.List+import Data.Maybe+import Data.Monoid+import GIS.Math.Projections+import GIS.Math.Spherical+import GIS.Utils+import Control.Lens.Tuple+import Control.Lens+import Control.Monad+import GIS.Types+import GIS.Graphics.Types+import Data.Default+import System.Directory++-- | 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)++-- 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)++-- | Given a projection and lists of districts, draw a map.+districtToMapP :: Projection -> [District] -> Map+districtToMapP p = projectMap p . districtToMap++-- | Given a list of districts, draw a map.+districtToMap :: [District] -> Map+districtToMap districts = labelledDistricts .~ dist $ def+    where dist = gc $ zip (fmap _shape districts) (fmap _districtLabel districts)+          gc = concatMap (\(a,b) -> zip a (replicate (length a) b))++-- | Given a projection and list of districts, return a list of maps.+districtToMapFilesP :: Projection -> [District] -> [Map]+districtToMapFilesP p = fmap (projectMap p) . districtToMapFiles++-- | Given a list of districts, return a list of maps.+districtToMapFiles :: [District] -> [Map]+districtToMapFiles = map (\(District polygons label _ _) -> title .~ label $ 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+getDistricts :: FilePath -> IO [District]+getDistricts filepath = do+        dbfExists <- doesFileExist (stripExt filepath <> ".dbf")+        file <- if dbfExists then readShpWithDbf filepath else readShpFile filepath+        let districtLabels = fromJust $ (fmap labels) $ mapM (shpRecLabel) . shpRecs $ file -- <$> for print? +        let shapes = (map (getPolygon . fromJust . shpRecContents)) . shpRecs $ file+        let perimeters = map (getPerimeter . getPolygon . fromJust . shpRecContents) . shpRecs $ file+        let areas = map (fmap areaPolyRectangular . getPolygon . fromJust . shpRecContents) . shpRecs $ file+        let tuple = zip4 shapes districtLabels perimeters areas+        pure $ fmap (\(a,b,c,d) -> District a b c d) tuple++-- | Given a list of polygons, return the total area.+getPerimeter :: [Polygon] -> Double+getPerimeter lines = sum $ fmap segmentLength lines+    where segmentLength [x1, x2]       = distance x1 x2+          segmentLength (x1:x2:points) = segmentLength (x2:points) + distance x1 x2++-- | Helper function for extracting from shapefiles.+getPolygon :: RecContents -> [Polygon]+getPolygon (RecPolygon { recPolPoints = pt }) = pt+getPolygon (RecPolygonM { recPolMPoints = pt }) = pt+getPolygon (RecPolygonZ { recPolZPoints = pt }) = pt
+ src/GIS/Math/Projections.hs view
@@ -0,0 +1,74 @@+-- | Module containing several useful projections for generating maps.+module GIS.Math.Projections where++import Control.Lens+import Control.Lens.Tuple+import GIS.Types+import GIS.Math.Utils+import GIS.Graphics.Types++-- | For use as a reference point in certain projections.+washingtonDC :: Point+washingtonDC = over _2 radians $ over _1 radians (38.9072, -77.0369)++-- | For use as a reference point with the Craig retroazimuthal projection+mecca :: Point+mecca = over _2 radians $ over _1 radians (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+-- (works on a subset of the world)+craig :: Point -> Projection+craig referencePoint (long, lat) = (long - referenceLong, y)+    where (referenceLong, referenceLat) = referencePoint+          y = if long - referenceLong == 0 then expr else (long - referenceLong)/(sin(long - referenceLong)) * expr+          expr = ((sin lat) * (cos (long - referenceLong)) - (tan referenceLat) * (cos lat)) ++-- | Winkel Tripel projection (standard for the National Geographic Society since +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. To make it usable you can+-- use+-- > 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 + 2 * n * (sin phi1)+          rho = sqrt (c - 2 * n * (sin lat)) / n+          rho' = sqrt (c - 2 * n * (sin referenceLat)) / n+          phi1 = (radians 20) -- standard parallels @ 20, 50 degrees+          phi2 = (radians 50)+          (referenceLong, referenceLat) = referencePoint ++-- | Helper to project given a `Polygon`.+project :: Projection -> Polygon -> Polygon+project f = fmap (f . toRadians)++-- | Helper to apply a projection given a `Map`.+projectMap :: Projection -> Map -> Map+projectMap p = over (labelledDistricts) (fmap (over _1 (project p)))
+ src/GIS/Math/Spherical.hs view
@@ -0,0 +1,56 @@+-- | Utilities to compute area, perimeter, etc. on the surface of a sphere.+module GIS.Math.Spherical where++import Control.Lens+import Control.Lens.Tuple+import Data.Function+import Data.List+import GIS.Types+import Data.Composition+import GIS.Math.Projections+import GIS.Math.Utils++-- | averages the coördinates of a polygon, returning a point.+shittyCentroid :: Polygon -> Point+shittyCentroid poly = (avg $ map fst poly, avg $ map snd poly)++-- | Average over a foldable container+avg :: (RealFrac a, Foldable t) => t a -> a+avg list = sum list / (fromIntegral . length $ list)++-- | Compute the area of a triangle using L'Huillier's formula+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)))+          s = (a + b + c) / 2+          a = distanceRad x1 x2+          b = distanceRad x1 x3+          c = distanceRad x2 x3+          distanceRad = on centralAngle toRadians++-- mandelbrot/fractal dimension? +-- consider "area of largest circumscribable circle" as well. ++-- | 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. +areaPolygon :: Polygon -> Double+areaPolygon = {--(*((6371*pi)^2)) . --}areaPolyRectangular . (fmap bonne)++-- | Find the area of a polygon with rectangular coördinates given. +areaPolyRectangular :: Polygon -> Double+areaPolyRectangular (pt:pts) = 0.5 * (fst (foldl' areaPolyCalc (0,pt) pts))+    where areaPolyCalc (sum,(x,y)) (xNext,yNext) = (sum + (x * yNext - xNext * y),(xNext,yNext))++-- | Distance in kilometers between two points given in degrees. +distance :: (Double, Double) -> (Double, Double) -> Double+distance = (*6371) .* (on centralAngle toRadians)++-- | Compute central angle from points given in radians+centralAngle :: (Double, Double) -> (Double, Double) -> Double+centralAngle (long1, lat1) (long2, lat2) = centralAngle+    where centralAngle = acos $ (sin lat1) * (sin lat2) + (cos lat1) * (cos lat2) * (cos (long1 - long2))
+ src/GIS/Math/Utils.hs view
@@ -0,0 +1,15 @@+-- | Miscellaneous utils for math+module GIS.Math.Utils where++import GIS.Types+import Control.Lens+import Control.Lens.Tuple++sinc x = sin x / x++-- | Convert a `Double` from degrees to radians.+radians = (*(pi/180))++-- | Convert both coördinates to radians.+toRadians :: Point -> Point+toRadians = (over _1 radians) . (over _2 radians)
+ src/GIS/Types.hs view
@@ -0,0 +1,30 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE TemplateHaskell   #-}++-- | Module containing types for various geometrical objects. Lenses included. +module GIS.Types where++import Text.PrettyPrint.ANSI.Leijen+import GHC.Generics+import Control.Lens++type Point = (Double, Double)+type Polygon = [Point]+type Projection = (Double, Double) -> (Double, Double)++-- | Data type for one record in a shape file, also capable of storing basic+-- information about the district. +data District = District { _shape         :: [Polygon]+                         , _districtLabel :: String+                         , _perimeter     :: Double+                         , _area          :: [Double]+                         } deriving (Generic, Show)++makeLenses ''District++-- possibly get rid of this idk+data DbfReadError = NotAPolygon | ShpNull++instance Show DbfReadError where+    show NotAPolygon = show $ red (text "Error: ") <> (text "Shape not a polygon! Are you sure you're opening a district?")+    show ShpNull = show $ red (text "Error: ") <> (text "")
+ src/GIS/Utils.hs view
@@ -0,0 +1,20 @@+-- | Miscellaneous utils and helper functions.+module GIS.Utils where++import Data.Char++-- | Strip extension from a filepath+stripExt :: String -> String+stripExt = reverse . drop 1 . dropWhile (/='.') . reverse++-- | Get extension given a filepath+getExt :: String -> String+getExt = fmap toLower . reverse . (takeWhile (/='.')) . reverse++-- | Make sure labels on shapefile objects are suitable for human reading.+labels :: (Show a) => [[a]] -> [String]+labels = map ((filter (/=' ')) . (take 35) . (drop 11) . show . (!!3))++-- | Flatten a nested tuple.+flatten :: ((a,b),c) -> (a,b,c)+flatten ((x,y),z) = (x,y,z)
+ test/Spec.hs view
@@ -0,0 +1,18 @@+import Test.Hspec+import GIS.Math.Spherical+import GIS.Hylo+import GIS.Graphics.PlotSVG+import GIS.Graphics.PlotPNG+import GIS.Math.Projections++main :: IO ()+main = hspec $ do+    describe "distance" $ do+        parallel $ it "computes distances along geodesics" $ do+            distance ((-74.0059), 40.7128) ((-118.2347), 34.0522) `shouldBe` 3934.997673163554+    describe "mercatorFullPng" $ do+        parallel $ it "Reads a map and writes a .png with the mercator projection" $ do+            (mkMapPng "test/testmap.png" =<< districtToMapP bonne <$> getDistricts "test/data/worldmap/TM_WORLD_BORDERS-0.3.shp") >>= (`shouldBe` ())+    describe "mercatorFullSVG" $ do+        parallel $ it "Reads a map and writes a .svg with the mercator projection" $ do+            (mkMapSVG "test/testmap.svg" =<< districtToMapP bonne <$> getDistricts "test/data/worldmap/TM_WORLD_BORDERS-0.3.shp") >>= (`shouldBe` ())