diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/app/MainPng.hs b/app/MainPng.hs new file mode 100644 --- /dev/null +++ b/app/MainPng.hs @@ -0,0 +1,6 @@ +module Main where + +import GIS.Exe.OptCairo + +main :: IO () +main = exec diff --git a/hgis.cabal b/hgis.cabal new file mode 100644 --- /dev/null +++ b/hgis.cabal @@ -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 diff --git a/src/GIS/Exe/Opt.hs b/src/GIS/Exe/Opt.hs new file mode 100644 --- /dev/null +++ b/src/GIS/Exe/Opt.hs @@ -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 diff --git a/src/GIS/Exe/OptCairo.hs b/src/GIS/Exe/OptCairo.hs new file mode 100644 --- /dev/null +++ b/src/GIS/Exe/OptCairo.hs @@ -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 diff --git a/src/GIS/Exe/Parser.hs b/src/GIS/Exe/Parser.hs new file mode 100644 --- /dev/null +++ b/src/GIS/Exe/Parser.hs @@ -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 diff --git a/src/GIS/Graphics/Plot.hs b/src/GIS/Graphics/Plot.hs new file mode 100644 --- /dev/null +++ b/src/GIS/Graphics/Plot.hs @@ -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 diff --git a/src/GIS/Graphics/PlotPNG.hs b/src/GIS/Graphics/PlotPNG.hs new file mode 100644 --- /dev/null +++ b/src/GIS/Graphics/PlotPNG.hs @@ -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 +--} diff --git a/src/GIS/Graphics/PlotSVG.hs b/src/GIS/Graphics/PlotSVG.hs new file mode 100644 --- /dev/null +++ b/src/GIS/Graphics/PlotSVG.hs @@ -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 +--} diff --git a/src/GIS/Graphics/Types.hs b/src/GIS/Graphics/Types.hs new file mode 100644 --- /dev/null +++ b/src/GIS/Graphics/Types.hs @@ -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 } diff --git a/src/GIS/Hylo.hs b/src/GIS/Hylo.hs new file mode 100644 --- /dev/null +++ b/src/GIS/Hylo.hs @@ -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 diff --git a/src/GIS/Math/Projections.hs b/src/GIS/Math/Projections.hs new file mode 100644 --- /dev/null +++ b/src/GIS/Math/Projections.hs @@ -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))) diff --git a/src/GIS/Math/Spherical.hs b/src/GIS/Math/Spherical.hs new file mode 100644 --- /dev/null +++ b/src/GIS/Math/Spherical.hs @@ -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)) diff --git a/src/GIS/Math/Utils.hs b/src/GIS/Math/Utils.hs new file mode 100644 --- /dev/null +++ b/src/GIS/Math/Utils.hs @@ -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) diff --git a/src/GIS/Types.hs b/src/GIS/Types.hs new file mode 100644 --- /dev/null +++ b/src/GIS/Types.hs @@ -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 "") diff --git a/src/GIS/Utils.hs b/src/GIS/Utils.hs new file mode 100644 --- /dev/null +++ b/src/GIS/Utils.hs @@ -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) diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -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` ())