{-# LANGUAGE RankNTypes #-}
module GIS.Hylo where
import Control.Lens hiding (lens)
import Data.Default
import Data.List
import Data.Maybe
import Data.Semigroup
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, perimeter)
import GIS.Utils
import System.Directory
districtArea :: [District] -> String
districtArea districts = concat . intercalate (pure "\n") $ fmap (pure . show . distA) districts
where distA (District _ label _ area _) = (label, sum area)
districtPerimeter :: [District] -> String
districtPerimeter districts = concat . intercalate (pure "\n") $ fmap (pure . show . distP) districts
where distP (District _ label perimeter _ _) = (label, perimeter)
districtCompactness :: [District] -> String
districtCompactness districts = concat . 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 [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
let shapes = fmap (getPolygon . fromJust . shpRecContents) . shpRecs $ file
let perimeters = fmap (totalPerimeter . getPolygon . fromJust . shpRecContents) . shpRecs $ file
let areas = fmap (fmap areaPolygon . getPolygon . fromJust . shpRecContents) . shpRecs $ file
let compacticity = fmap (relativeCompactness . concat . getPolygon . fromJust . shpRecContents) . shpRecs $ file
pure $ zipWith5 District shapes districtLabels perimeters areas compacticity
getPolygon :: RecContents -> [Polygon]
getPolygon RecPolygon { recPolPoints = pt } = pt
getPolygon RecPolygonM { recPolMPoints = pt } = pt
getPolygon RecPolygonZ { recPolZPoints = pt } = pt
getPolygon _ = error "shapefile contents not supported. please ensure the shapefile contains a list of polygons."