{-# LANGUAGE RankNTypes #-} {-# LANGUAGE RecordWildCards #-} -- | Module to transform shapefiles. module GIS.Hylo ( districtArea , districtPerimeter , districtCompactness , districtToMapP , districtToMapLensP , districtToMapLens , districtToMap , districtToMapFilesP , getDistricts , getPolygon ) where import Control.Lens hiding (lens) import Data.Default import Data.Foldable (fold) import Data.List import Data.Maybe import Geometry.Shapefile.MergeShpDbf import Geometry.Shapefile.ReadShp import Geometry.Shapefile.Types import GIS.Graphics.Types import GIS.Math.Projections import GIS.Math.Spherical import GIS.Types hiding (area) import GIS.Utils import System.Directory -- | Get the areas of various objects and return a string suitable for printing districtArea :: [District] -> String districtArea districts = fold . intercalate (pure "\n") $ fmap (pure . show . distA) districts where distA (District _ label _ area _) = (label, sum area) -- TODO figure out which one is the correct one -- | Get the perimeters of various objects and return a string suitable for printing districtPerimeter :: [District] -> String districtPerimeter districts = fold . intercalate (pure "\n") $ fmap (pure . show . distP) districts where distP (District _ label perimeter _ _) = (label, perimeter) -- | Label with relative compactness districtCompactness :: [District] -> String districtCompactness districts = fold . intercalate (pure "\n") $ fmap (pure . show . distC) districts where distC (District _ label _ _ compacticity) = (label, compacticity) -- | Given a projection and lists of districts, draw a map. districtToMapP :: Projection -> [District] -> Map districtToMapP p = projectMap p . districtToMap -- | Given a projection and lists of districts, draw a map, with labels determined by a lens districtToMapLensP :: (Show a) => Projection -> Lens' District a -> [District] -> Map districtToMapLensP p lens districts = projectMap p (districtToMapLens lens districts) -- | Given a list of districts, draw a map, with labels determined by a lens. districtToMapLens :: (Show a) => Lens' District a -> [District] -> Map districtToMapLens lens districts = labelledDistricts .~ dist $ def where dist = zip (fmap _shape districts) (fmap (show . view lens) districts) -- | Given a list of districts, draw a map. districtToMap :: [District] -> Map districtToMap = districtToMapLens districtLabel -- | 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 = 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] 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 shapes = extract id perimeters = extract totalPerimeter areas = extract (fmap areaPolygon) compacticity = extract (relativeCompactness . fold) pure $ zipWith5 District shapes districtLabels perimeters areas compacticity -- | Helper function for extracting from shapefiles. getPolygon :: RecContents -> [Polygon] getPolygon RecPolygon {..} = recPolPoints getPolygon RecPolygonM {..} = recPolMPoints getPolygon RecPolygonZ {..} = recPolZPoints getPolygon _ = error "shapefile contents not supported. please ensure the shapefile contains a list of polygons."