module Geometry.Shapefile.ReadShp ( readShpFile,
readShpData,
) where
import Control.Monad (replicateM, replicateM_, when)
import Control.Monad.Loops (whileM)
import Data.Binary.Get
import Data.Binary.IEEE754 (getFloat64le)
import Data.Semigroup
import qualified Data.ByteString.Lazy as BL
import Geometry.Shapefile.Internal
import Geometry.Shapefile.Types
readShpFile :: String -> IO ShpData
readShpFile fp = readShpData <$> BL.readFile fp
readShpData :: BL.ByteString -> ShpData
readShpData = runGet $ do
shpH <- getShpHeader
shpRs <- whileM (not <$> isEmpty) getShpRec
pure ShpData { shpHeader = shpH,
dbfFieldDescs = Nothing,
shpRecs = shpRs }
getShpHeader :: Get ShpHeader
getShpHeader = do
fc <- getIntBE
when (fc /= 9994) (error "getShpHeader: Input file is not a SHP file")
_ <- replicateM_ 5 getWord32be
shpLen <- getIntBE
shpV <- getIntLE
shpT <- getShpType
shpbb <- getShpBBox shpT
pure ShpHeader { shpFileLength = shpLen,
shpVersion = shpV,
shpType = shpT,
shpBB = shpbb }
getShpType :: Get ShpType
getShpType = shpTypeFromId <$> getIntLE
getShpBBox :: ShpType -> Get ShpBBox
getShpBBox t = do
[xMin, yMin, xMax, yMax] <- replicateM 4 getFloat64le
[zMin, zMax] <- let ds = replicateM 2 getFloat64le
in if t `elem` zTypes
then fmap Just <$> ds
else ds >> pure [Nothing, Nothing]
[mMin, mMax] <- let ds = replicateM 2 getFloat64le
in if t `elem` mTypes
then fmap Just <$> ds
else ds >> pure [Nothing, Nothing]
pure ShpBBox { shpXMin = xMin,
shpXMax = xMax,
shpYMin = yMin,
shpYMax = yMax,
shpZMin = zMin,
shpZMax = zMax,
shpMMin = mMin,
shpMMax = mMax }
getShpRec :: Get ShpRec
getShpRec = do
recNum <- getIntBE
recLen <- getIntBE
recType <- getShpType
recContents <- if recType == ShpNull
then pure Nothing
else Just <$> getRecContents recType
pure ShpRec { shpRecNum = recNum,
shpRecLen = recLen,
shpRecContents = recContents,
shpRecLabel = Nothing,
shpRecType = recType }
getRecContents :: ShpType -> Get RecContents
getRecContents t = case t of
ShpPoint -> RecPoint <$> getPoint
ShpPointM -> do
p <- getPoint
m <- getFloat64le
pure RecPointM { pmPoint = p, pmM = m }
ShpPointZ -> do
p <- getPoint
z <- getFloat64le
m <- getFloat64le
pure RecPointZ { pzPoint = p, pzZ = z, pzM = m }
ShpMultiPoint -> do
(bb, nPoints, points) <- getPointsData
pure RecMultiPoint { recMPBBox = bb,
recMPNumPoints = nPoints,
recMPPoints = points }
ShpMultiPointM -> do
(bb, nPoints, points) <- getPointsData
(mMin, mMax, ms) <- getMData nPoints
pure RecMultiPointM { recMPMBBox = bb,
recMPMNumPoints = nPoints,
recMPMPoints = points,
recMPMMRange = (mMin, mMax),
recMPMMs = ms}
ShpMultiPointZ -> do
(bb, nPoints, points) <- getPointsData
(zMin, zMax, zs) <- getZData nPoints
(mMin, mMax, ms) <- getMData nPoints
pure RecMultiPointZ { recMPZBBox = bb,
recMPZNumPoints = nPoints,
recMPZPoints = points,
recMPZZRange = (zMin, zMax),
recMPZZs = zs,
recMPZMRange = (mMin, mMax),
recMPZMs = ms }
ShpPolyLine -> do
(bb, nParts, nPoints, partLengths, pntLists) <- getPolyData
pure RecPolyLine { recPolLBBox = bb,
recPolLNumParts = nParts,
recPolLNumPoints = nPoints,
recPolLPartLengths = partLengths,
recPolLPoints = pntLists }
ShpPolyLineM -> do
(bb, nParts, nPoints, partLengths, pntLists) <- getPolyData
(mMin, mMax, ms) <- getMData nPoints
pure RecPolyLineM { recPolLMBBox = bb,
recPolLMNumParts = nParts,
recPolLMNumPoints = nPoints,
recPolLMPartLengths = partLengths,
recPolLMPoints = pntLists,
recPolLMMRange = (mMin, mMax),
recPolLMMs = ms }
ShpPolyLineZ -> do
(bb, nParts, nPoints, partLengths, pntLists) <- getPolyData
(zMin, zMax, zs) <- getZData nPoints
(mMin, mMax, ms) <- getMData nPoints
pure RecPolyLineZ { recPolLZBBox = bb,
recPolLZNumParts = nParts,
recPolLZNumPoints = nPoints,
recPolLZPartLengths = partLengths,
recPolLZPoints = pntLists,
recPolLZZRange = (zMin, zMax),
recPolLZZs = zs,
recPolLZMRange = (mMin, mMax),
recPolLZMs = ms }
ShpPolygon -> do
(bb, nParts, nPoints, partLengths, pntLists) <- getPolyData
pure RecPolygon { recPolBBox = bb,
recPolNumParts = nParts,
recPolNumPoints = nPoints,
recPolPartLengths = partLengths,
recPolPoints = pntLists }
ShpPolygonM -> do
(bb, nParts, nPoints, partLengths, pntLists) <- getPolyData
(mMin, mMax, ms) <- getMData nPoints
pure RecPolygonM { recPolMBBox = bb,
recPolMNumParts = nParts,
recPolMNumPoints = nPoints,
recPolMPartLengths = partLengths,
recPolMPoints = pntLists,
recPolMMRange = (mMin, mMax),
recPolMMs = ms }
ShpPolygonZ -> do
(bb, nParts, nPoints, partLengths, pntLists) <- getPolyData
(zMin, zMax, zs) <- getZData nPoints
(mMin, mMax, ms) <- getMData nPoints
pure RecPolygonZ { recPolZBBox = bb,
recPolZNumParts = nParts,
recPolZNumPoints = nPoints,
recPolZPartLengths = partLengths,
recPolZPoints = pntLists,
recPolZZRange = (zMin, zMax),
recPolZZs = zs,
recPolZMRange = (mMin, mMax),
recPolZMs = ms }
ShpNull ->
pure RecNull
ShpMultiPatch ->
error "getShpRec: MultiPatch type is not supported, sorry!"
getPointsData :: Get (RecBBox, Int, [Point])
getPointsData = do
bb <- getRecBBox
nPoints <- getIntLE
points <- replicateM nPoints getPoint
pure (bb, nPoints, points)
getPolyData :: Get (RecBBox, Int, Int, [Int], [[Point]])
getPolyData = do
bb <- getRecBBox
nParts <- getIntLE
nPoints <- getIntLE
partIndices <- replicateM nParts getIntLE
let partLengths = steps $ partIndices <> [nPoints]
pntLists <- mapM getPointList partLengths
pure (bb, nParts, nPoints, partLengths, pntLists)
getMData :: Int -> Get (Double, Double, [Double])
getMData nPoints = do
[mMin, mMax] <- replicateM 2 getFloat64le
ms <- replicateM nPoints getFloat64le
pure (mMin, mMax, ms)
getZData :: Int -> Get (Double, Double, [Double])
getZData = getMData
getRecBBox :: Get RecBBox
getRecBBox = do
[xMin, yMin, xMax, yMax] <- replicateM 4 getFloat64le
pure RecBBox { recXMin = xMin,
recXMax = xMax,
recYMin = yMin,
recYMax = yMax }