{-# LANGUAGE OverloadedStrings, ScopedTypeVariables, ViewPatterns, PatternSynonyms #-}
{-|
Module: Data.LatLong
Description: Spatially indexed type for geographic coordinates.
Copyright: © 2018-2019 Satsuma labs
Defines a type for georgraphic coordinates that can be spatially indexed by any database supporting 64 bit integer values.
This indexing works by reperesenting points using a 'Morton' Z-Order curce, with each coordinate reperesented as a 32-bit fixed-point value
which then have their bits interleaved into a 64-bit integer to the internal reperesentation.
Taking binary prefixes of these values divides the globe into a hierarchy of rectangular tiles (repereseteh here as 'LatLongTile' objects),
each of which is a contiguous interval when points are ordered according to their integer reperesentations.
As any geographic region can be covered by a small number of tiles of simillar size, this provides an easy to loop up data for specific reguions.
Instances and a filter for persistent are provided for this purpose.
-}
module Data.LatLong (
LatLong(LatLongZ, LatLong), lat, long,
earthRadius, geoDistance, geoSquare,
-- * Tiles
LatLongTile, latLongTileInterval,
latLongTileCover, latLongTileCoverSquare, tileSetElem, withinTileSet
) where
import Data.Morton
import Data.Aeson
import Data.Proxy
import qualified Data.Text as T
import Control.Monad
import Control.Lens (Lens')
import Data.Word
import Numeric
import Web.HttpApiData
import Database.Persist.Sql
-- | Type for storing geographic coordinates that can be spatially indexed ('Morton' ordering).
-- Each coordinate is reperesented as as 32-bit fixed point value and is also accessible as a Double through a pattern synonym.
-- Order follows a Morton Z-order curve which can be used to search a database by tiles.
-- This works with any database capable of storing and indexing 'Word64' (although this type only uses those values fitting in a 64 bit signed integer)
newtype LatLong =
-- | Underlying reperesentation and source of ordering for indexing
LatLongZ Morton
deriving (Eq, Ord)
two32f :: Double
two32f = 2 ^ (32 :: Int)
makeLatLong :: Double -> Double -> LatLong
makeLatLong theta phi = LatLongZ (MortonPair theta' phi') where
theta' = clampLat . floor $ (theta + 90) / 360 * two32f
phi' = wrapLong . floor $ (phi + 180) / 360 * two32f
clampLat (x::Int) = fromIntegral (max 0 (min 0x7fffffff x))
wrapLong (x::Int) = fromIntegral (x `mod` 0x100000000)
latLongCoords :: LatLong -> (Double, Double)
latLongCoords (LatLongZ (MortonPair theta' phi')) = (theta,phi) where
theta = (fromIntegral theta' * 360 / two32f) - 90
phi = (fromIntegral phi' * 360 / two32f) - 180
-- | Pattern for accessing latitide and longitude coordinates as 'Double' values.
-- This is not fully isomoprphic as latitude is clipped to ±90, longitude is wrapped mod 360 ±180,
-- and rounding error exists due to the internal fixed-point reperesentation.
pattern LatLong :: Double -> Double -> LatLong
pattern LatLong theta phi <- (latLongCoords -> (theta,phi)) where
LatLong theta phi = makeLatLong theta phi
{-# COMPLETE LatLong #-}
instance ToJSON LatLong where
toJSON (LatLong theta phi) = object ["lat" .= theta, "long" .= phi]
toEncoding (LatLong theta phi) = pairs $ "lat" .= theta <> "long" .= phi
instance FromJSON LatLong where
parseJSON = withObject "LatLong" $ \o -> do
theta <- o .: "lat"
phi <- o .: "long"
guard $ theta > -90 && theta < 90
guard $ phi >= -180 && phi < 180
return $ LatLong theta phi
instance FromHttpApiData LatLong where
parseUrlPiece s = maybe (Left "malformed coordinate pair") Right $ do
[theta',phi'] <- return $ T.splitOn "," s
Right theta <- return $ parseUrlPiece theta'
Right phi <- return $ parseUrlPiece phi'
guard $ theta >= -90 && theta <= 90
guard $ phi >= -180 && phi <= 180
return $ LatLong theta phi
instance ToHttpApiData LatLong where
toUrlPiece (LatLong theta phi) = toUrlPiece theta <> "," <> toUrlPiece phi
instance Show LatLong where
show (LatLong theta phi) = join
[ showFFloat (Just 5) (abs theta) []
, if theta >= 0 then " N " else " S "
, showFFloat (Just 5) (abs phi) []
, if phi >= 0 then " E" else " W" ]
instance PersistField LatLong where
toPersistValue (LatLongZ (Morton x)) = toPersistValue x
fromPersistValue = fmap (LatLongZ . Morton) . fromPersistValue
instance PersistFieldSql LatLong where
sqlType _ = sqlType (Proxy :: Proxy Word64)
-- | Lens for latitude.
lat :: Lens' LatLong Double
lat f (LatLong theta phi) = fmap (\theta' -> LatLong theta' phi) (f theta)
-- | Lens for longitude.
long :: Lens' LatLong Double
long f (LatLong theta phi) = fmap (\phi' -> LatLong theta phi') (f phi)
-- | Earth's average radius in meters
earthRadius :: Double
earthRadius = 6371.2e3
rads :: Double->Double
rads x = x / 180 * pi
degs :: Double->Double
degs x = x * 180 / pi
sindeg :: Double->Double
sindeg = sin . rads
cosdeg :: Double->Double
cosdeg = cos . rads
-- | Calculate distance between two points using the Haversine formula (up to 0.5% due to the assumption of a spherical Earth).
-- Distance is returned in meters.
geoDistance :: LatLong -> LatLong -> Double
geoDistance (LatLong theta1 phi1) (LatLong theta2 phi2) = earthRadius * sigma where
havdelta x y = sindeg ((x-y)/2) ^ (2::Int)
hav = havdelta theta1 theta2 + (cosdeg theta1 * cosdeg theta2 * havdelta phi1 phi2)
sigma = 2 * asin (sqrt hav)
-- | Calculates the corner coordinates of a square with a given center and radius (in meters).
-- Based on the Mercator projection thus has distortion near the poles
-- (within 5% for a radius at most 200km and latitude within ±70).
geoSquare :: LatLong -> Double -> (LatLong, LatLong)
geoSquare (LatLong theta phi) r = let
dtheta = degs (r / earthRadius)
dphi = degs (r / earthRadius / cosdeg theta)
se = LatLong (theta - dtheta) (phi - dphi)
nw = LatLong (theta + dtheta) (phi + dphi)
in (se,nw)
-- | Represents a LatLong tile, which is both a rectangle and a contoguous interval in the ordering.
newtype LatLongTile = LatLongTile MortonTile deriving (Eq, Read, Show)
-- | Gets the corners of a tile, which are also the bounds of its interval in sort order.
latLongTileInterval :: LatLongTile -> Interval LatLong
latLongTileInterval (LatLongTile t) = fmap LatLongZ (mortonTileBounds t)
-- | Covers a rectangle (defined by its corners) tiles of at most its size.
latLongTileCover :: LatLong -> LatLong -> [LatLongTile]
latLongTileCover se nw = let
LatLong s _ = se
LatLong n _ = nw
LatLongZ y = se
LatLongZ x = nw
in if s > n then [] else fmap LatLongTile (mortonTileCoverTorus x y)
-- | Covers a square (defined by its center and radius) by tiles.
latLongTileCoverSquare :: LatLong -> Double -> [LatLongTile]
latLongTileCoverSquare c r = uncurry latLongTileCover $ geoSquare c r
-- | Tests whether a point is contasined in a tile set.
tileSetElem :: LatLong -> [LatLongTile] -> Bool
tileSetElem p ts = or [intervalElem p (latLongTileInterval t) | t <- ts]
-- | Persistent filter producing the SQL equiveland ot 'tileSetElem'.
withinTileSet :: (EntityField row LatLong) -> [LatLongTile] -> Filter row
withinTileSet field tiles = let
tfilter tile = let Interval a b = latLongTileInterval tile in FilterAnd [field >=. a, field <=. b]
in FilterOr $ fmap tfilter tiles