{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE RecordWildCards #-}
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 Math.Geometry.Spherical
import System.Directory
districtArea :: [District] -> String
districtArea districts = fold . intercalate (pure "\n") $ fmap (pure . show . distA) districts
where distA (District _ label _ area _) = (label, sum area)
districtPerimeter :: [District] -> String
districtPerimeter districts = fold . intercalate (pure "\n") $ fmap (pure . show . distP) districts
where distP (District _ label perimeter _ _) = (label, perimeter)
districtCompactness :: [District] -> String
districtCompactness districts = fold . intercalate (pure "\n") $ fmap (pure . show . distC) districts
where distC (District _ label _ _ compacticity) = (label, compacticity)
districtToMapP :: Projection -> [District] -> Map
districtToMapP p = projectMap p . districtToMap
districtToMapLensP :: (Show a) => Projection -> Lens' District a -> [District] -> Map
districtToMapLensP p lens districts = projectMap p (districtToMapLens lens districts)
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)
districtToMap :: [District] -> Map
districtToMap = districtToMapLens districtLabel
districtToMapFilesP :: Projection -> [District] -> [Map]
districtToMapFilesP p = fmap (projectMap p) . districtToMapFiles
districtToMapFiles :: [District] -> [Map]
districtToMapFiles = fmap g where
g (District polygons label _ area' _) = title .~ label <> "-" <> (show . sum $ area') $ labelledDistricts .~ pure (polygons,"") $ def
getDistricts :: FilePath -> IO (Maybe [District])
getDistricts filepath = do
dbfExists <- doesFileExist (stripExt filepath <> ".dbf")
file <- if dbfExists then readShpWithDbf filepath else readShpFile filepath
let extract f = fmap (f . getPolygon . fromJust . shpRecContents) . shpRecs $ file
districtLabels = fmap labels $ traverse shpRecLabel . shpRecs $ file
shapes = extract id
perimeters = extract totalPerimeter
areas = extract (fmap (areaPolygon 6371))
compacticity = extract (relativeCompactness . fold)
pure $ zipWith5 District shapes <$> districtLabels <*> pure perimeters <*> pure areas <*> pure compacticity
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."