Hylo.hs

{-# 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."