diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,5 @@ +# hgis + +## 1.0.0.3 + + * Use `spherical` package diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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 diff --git a/Setup.hs b/Setup.hs deleted file mode 100644 --- a/Setup.hs +++ /dev/null @@ -1,2 +0,0 @@ -import Distribution.Simple -main = defaultMain diff --git a/hgis.cabal b/hgis.cabal --- a/hgis.cabal +++ b/hgis.cabal @@ -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 diff --git a/src/GIS/Graphics/Plot.hs b/src/GIS/Graphics/Plot.hs --- a/src/GIS/Graphics/Plot.hs +++ b/src/GIS/Graphics/Plot.hs @@ -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. diff --git a/src/GIS/Hylo.hs b/src/GIS/Hylo.hs --- a/src/GIS/Hylo.hs +++ b/src/GIS/Hylo.hs @@ -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] diff --git a/src/GIS/Math/Projections.hs b/src/GIS/Math/Projections.hs --- a/src/GIS/Math/Projections.hs +++ b/src/GIS/Math/Projections.hs @@ -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 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 @@ -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) diff --git a/src/GIS/Math/Utils.hs b/src/GIS/Math/Utils.hs deleted file mode 100644 --- a/src/GIS/Math/Utils.hs +++ /dev/null @@ -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 diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -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` ())