diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,11 @@ # Change log for [naqsha] +## [0.2.0.0] - 22 July, 2017 + +* Overall change in module structure. + - Naqsha.Position becomes Naqsha.Geometry.Coordinate +* Support for geohash + ## [0.1.0.0] - 6th May, 2017 Very first release of naqsha. @@ -10,4 +16,5 @@ [naqsha]: "Naqsha library" +[0.2.0.0]: "Release 0.2.0.0" [0.1.0.0]: "Release 0.1.0.0" diff --git a/Naqsha.hs b/Naqsha.hs deleted file mode 100644 --- a/Naqsha.hs +++ /dev/null @@ -1,9 +0,0 @@ --- | Naqsha is a library to work with geo data. This module exposes --- the basic types and some utility functions. -module Naqsha - ( module Naqsha.Position - , module Naqsha.Geometry.Angle - ) where - -import Naqsha.Position -import Naqsha.Geometry.Angle diff --git a/Naqsha/Geometry.hs b/Naqsha/Geometry.hs new file mode 100644 --- /dev/null +++ b/Naqsha/Geometry.hs @@ -0,0 +1,62 @@ +-- | The geometric types and values exposed by naqsha. +module Naqsha.Geometry + ( module Naqsha.Geometry.Coordinate + + -- ** Angles and angular quantities. + , module Naqsha.Geometry.Angle + + -- * Geometric hashing. + -- $geohashing$ + + -- * Distance calculation. + -- + -- $distance$ + -- + + -- * Internal details + -- $internals$ + ) where + +import Naqsha.Geometry.Angle +import Naqsha.Geometry.Coordinate + +-- Nothing imported here. Only for docs. +import Naqsha.Geometry.Spherical() + + +-- $geohashing$ +-- +-- Geometric hashing is a technique of converting geometric +-- coordinates into 1-dimension strings. Often these hashes ensures +-- that string with large common prefix are close by (although not the +-- converse). Hence, these hashes can be used to stored geo-cordinates +-- in database and build into it a sense of location awareness. We support +-- the following geometric hashing: +-- +-- ["Naqsha.Geometry.Coordinate.GeoHash":] The geohash standard +-- (). +-- +-- None of these modules are imported by default the user may import +-- the one that is most desirable. + +-- $distance$ +-- +-- Calculating quantities like distance, bearing etc depends on the +-- model of the globe that we choose. Even in a given model we might +-- have different algorithms to compute the distance depending on +-- speed-accuracy trade-offs. Choosing the correct model and +-- algorithms is application dependent and hence we do not expose any +-- default ones. The following modules can be imported depending on the need +-- +-- ["Naqsha.Geometry.Spherical": ] Assume a spherical model of the +-- globe. Distance is calculated using the haversine formula. + + +-- $internals$ +-- +-- The basic types like `Latitude` or `Longitude` are exposed as +-- opaque types from this module. For type safety, we encourage the +-- users to use this module mostly when dealing with those times. For +-- the rare case when some non-trivial operations need to be defined, +-- we expose the internal module "Naqsha.Geometry.Internal". However, +-- use this interface with caution. diff --git a/Naqsha/Geometry/Angle.hs b/Naqsha/Geometry/Angle.hs --- a/Naqsha/Geometry/Angle.hs +++ b/Naqsha/Geometry/Angle.hs @@ -1,8 +1,3 @@ -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE CPP #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE TypeFamilies #-} - -- | Basic types associated with geometry. module Naqsha.Geometry.Angle ( Angle @@ -12,85 +7,7 @@ , Angular(..) ) where - -import Control.Monad ( liftM ) -import Data.Default -import Data.Group -import Data.Int -import GHC.Real -import Data.Vector.Unboxed ( MVector(..), Vector, Unbox) -import qualified Data.Vector.Generic as GV -import qualified Data.Vector.Generic.Mutable as GVM - - ------------------------------ Angles and Angular quantities ----------------------- - --- | An abstract angle. Internally, angles are represented as a 64-bit --- integer with each unit contribute 1/2^64 fraction of a complete --- circle. This means that angles are accurate up to a resolution of 2 --- π / 2^64 radians. Angles form a group under the angular addition --- and the fact that these are represented as integers means one can --- expect high speed accurate angle arithmetic. --- --- When expressing angles one can use a more convenient notation: --- --- > myAngle = degree 21.71167 --- > yourAngle = degree 21 <> minute 42 <> second 42 --- -newtype Angle = Angle {unAngle :: Int64} deriving (Enum, Eq, Ord, Unbox, Show, Read) - --- | Express angle in degrees. -degree :: Rational -> Angle -degree = Angle . fromInteger . round . (*scale) - where scale = (2^(64:: Int)) % 360 - --- | Express angle in minutes. -minute :: Rational -> Angle -minute = degree . (*scale) - where scale = 1 % 60 - --- | Express angle in seconds. -second :: Rational -> Angle -second = degree . (*scale) - where scale = 1 % 3600 - --- | Express angle in radians -radian :: Double -> Angle -radian = Angle . round . (*scale) - where scale = (2^(63:: Int)) / pi - - ----------------------- Decimal representation of angle ---------------------------------- - --- | Measure angle in degrees. This conversion may lead to loss of --- precision. -toDegree :: Fractional r => Angle -> r -toDegree = fromRational . (*conv) . toRational . unAngle - where conv = 360 % (2^(64 :: Int)) - --- | Measure angle in radians. This conversion may lead to loss of --- precision. -toRadian :: Angle -> Double -toRadian = (*conv) . fromIntegral . unAngle - where conv = pi / (2^(63:: Int)) - -instance Default Angle where - def = Angle 0 - -instance Angular Angle where - toAngle = id - -instance Monoid Angle where - mempty = Angle 0 - mappend (Angle x) (Angle y) = Angle $ x + y - mconcat = Angle . sum . map unAngle - -instance Group Angle where - invert (Angle x) = Angle $ negate x - -instance Bounded Angle where - maxBound = Angle maxBound - minBound = Angle minBound +import Naqsha.Geometry.Internal ------------------------------ The angular class ------------------------ @@ -98,55 +15,12 @@ class Angular a where toAngle :: a -> Angle -------------------- Making stuff suitable for unboxed vector. -------------------------- - -newtype instance MVector s Angle = MAngV (MVector s Int64) -newtype instance Vector Angle = AngV (Vector Int64) +instance Angular Angle where + toAngle = id -instance GVM.MVector MVector Angle where - {-# INLINE basicLength #-} - {-# INLINE basicUnsafeSlice #-} - {-# INLINE basicOverlaps #-} - {-# INLINE basicUnsafeNew #-} - {-# INLINE basicUnsafeReplicate #-} - {-# INLINE basicUnsafeRead #-} - {-# INLINE basicUnsafeWrite #-} - {-# INLINE basicClear #-} - {-# INLINE basicSet #-} - {-# INLINE basicUnsafeCopy #-} - {-# INLINE basicUnsafeGrow #-} - basicLength (MAngV v) = GVM.basicLength v - basicUnsafeSlice i n (MAngV v) = MAngV $ GVM.basicUnsafeSlice i n v - basicOverlaps (MAngV v1) (MAngV v2) = GVM.basicOverlaps v1 v2 - - basicUnsafeRead (MAngV v) i = Angle `liftM` GVM.basicUnsafeRead v i - basicUnsafeWrite (MAngV v) i (Angle x) = GVM.basicUnsafeWrite v i x - - basicClear (MAngV v) = GVM.basicClear v - basicSet (MAngV v) (Angle x) = GVM.basicSet v x - - basicUnsafeNew n = MAngV `liftM` GVM.basicUnsafeNew n - basicUnsafeReplicate n (Angle x) = MAngV `liftM` GVM.basicUnsafeReplicate n x - basicUnsafeCopy (MAngV v1) (MAngV v2) = GVM.basicUnsafeCopy v1 v2 - basicUnsafeGrow (MAngV v) n = MAngV `liftM` GVM.basicUnsafeGrow v n - -#if MIN_VERSION_vector(0,11,0) - basicInitialize (MAngV v) = GVM.basicInitialize v -#endif - -instance GV.Vector Vector Angle where - {-# INLINE basicUnsafeFreeze #-} - {-# INLINE basicUnsafeThaw #-} - {-# INLINE basicLength #-} - {-# INLINE basicUnsafeSlice #-} - {-# INLINE basicUnsafeIndexM #-} - {-# INLINE elemseq #-} - basicUnsafeFreeze (MAngV v) = AngV `liftM` GV.basicUnsafeFreeze v - basicUnsafeThaw (AngV v) = MAngV `liftM` GV.basicUnsafeThaw v - basicLength (AngV v) = GV.basicLength v - basicUnsafeSlice i n (AngV v) = AngV $ GV.basicUnsafeSlice i n v - basicUnsafeIndexM (AngV v) i = Angle `liftM` GV.basicUnsafeIndexM v i +instance Angular Latitude where + toAngle = unLat - basicUnsafeCopy (MAngV mv) (AngV v) = GV.basicUnsafeCopy mv v - elemseq _ (Angle x) = GV.elemseq (undefined :: Vector a) x +instance Angular Longitude where + toAngle = unLong diff --git a/Naqsha/Geometry/Coordinate.hs b/Naqsha/Geometry/Coordinate.hs new file mode 100644 --- /dev/null +++ b/Naqsha/Geometry/Coordinate.hs @@ -0,0 +1,195 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE Rank2Types #-} +-- | This module captures position of a point on the globe. +module Naqsha.Geometry.Coordinate + ( -- * Basics + -- $latandlong$ + Geo(..) + , northPole, southPole + -- ** Latitudes + , Latitude + , lat, north, south + , equator + , tropicOfCancer + , tropicOfCapricon + -- ** Longitudes. + , Longitude + , lon, east, west + , greenwich + ) where + +import Control.Monad ( liftM ) +import Data.Group +import Data.Vector.Unboxed ( MVector(..), Vector) +import qualified Data.Vector.Generic as GV +import qualified Data.Vector.Generic.Mutable as GVM + + +import Prelude -- To avoid redundunt import warnings. + +import Naqsha.Geometry.Angle +import Naqsha.Geometry.Internal + + + +-- $latandlong$ +-- +-- A point on the globe is specified by giving its geo coordinates +-- represented by the type `Geo`. It is essentially a pair of the +-- `Latitude` and `Longitude` of the point. +-- +-- == Examples +-- +-- > kanpurLatitude :: Latitude +-- > kanpurLatitude = lat $ degree 26.4477777 +-- > kanpurLongitude :: Longitude +-- > kanpurLongitude = lon $ degree 80.3461111 +-- > kanpurGeo :: Geo +-- > kanpurGeo = Geo kanpurLatitude kanpurLongitude +-- +-- You can also specify the latitude and longitude in units of degree, +-- minute and seconds. +-- +-- > kanpurLatitude = lat $ degree 26 <> minute 26 <> second 52 +-- > kanpurLongitude = lon $ degree 80 <> minute 20 <> second 46 +-- +-- The show and read instance of the `Latitude` and `Longitude` types +-- uses degrees for displaying and reading respectively. Show and Read +-- instances can express these quantities up to Nano degree precision. +-- +-- == Convention on sign. +-- +-- For latitudes, positive means north of the equator and negative +-- means south. In the case of longitudes, positive means east of the +-- longitude zero and negative means west. However, if you find these +-- conventions confusing you can use the combinators `north`, `south`, +-- `east`, and `west` when constructing latitudes or longitudes. +-- + + +-- | Convert an angle to a northern latitude +-- +-- > tropicOfCancer = north $ degree 23.5 +-- +north :: Angle -> Latitude +north = lat + +-- | Convert an angle to a southern latitude. +-- +-- > tropicOfCapricon = south $ degree 23.5 +-- +south :: Angle -> Latitude +south = lat . invert + + +-- | The latitude of equator. +equator :: Latitude +equator = lat $ degree 0 + +-- | The latitude corresponding to the Tropic of Cancer. +tropicOfCancer :: Latitude +tropicOfCancer = north $ degree 23.5 + +-- | The latitude corresponding to the Tropic of Capricon +tropicOfCapricon :: Latitude +tropicOfCapricon = south $ degree 23.5 + +-- | Convert angle to an eastern longitude. +-- +-- > kanpurLongitude = east $ degree 80.3461 +-- +east :: Angle -> Longitude +east = lon + +-- | Convert angle to a western longitude +-- +-- > newyorkLongitude = west $ degree 74.0059 +-- +west :: Angle -> Longitude +west = lon . invert + + +-- | The zero longitude. +greenwich :: Longitude +greenwich = lon $ degree 0 + + +------------------- The geometric coordinates. ----------------- + +-- | The coordinates of a point on the earth's surface. +data Geo = Geo {-# UNPACK #-} !Latitude + {-# UNPACK #-} !Longitude + deriving Show + + +-- | The North pole +northPole :: Geo +northPole = Geo maxBound $ lon $ degree 0 + +-- | The South pole +southPole :: Geo +southPole = Geo minBound $ lon $ degree 0 + +instance Eq Geo where + (==) (Geo xlat xlong) (Geo ylat ylong) + | xlat == maxBound = ylat == maxBound -- longitude irrelevant for north pole + | xlat == minBound = ylat == minBound -- longitude irrelevant for south pole + | otherwise = xlat == ylat && xlong == ylong + +----------------------------- Vector Instance for Geo --------------------------------------------- + + +newtype instance MVector s Geo = MGeoV (MVector s (Angle,Angle)) +newtype instance Vector Geo = GeoV (Vector (Angle,Angle)) + + +instance GVM.MVector MVector Geo where + {-# INLINE basicLength #-} + {-# INLINE basicUnsafeSlice #-} + {-# INLINE basicOverlaps #-} + {-# INLINE basicUnsafeNew #-} + {-# INLINE basicUnsafeReplicate #-} + {-# INLINE basicUnsafeRead #-} + {-# INLINE basicUnsafeWrite #-} + {-# INLINE basicClear #-} + {-# INLINE basicSet #-} + {-# INLINE basicUnsafeCopy #-} + {-# INLINE basicUnsafeGrow #-} + basicLength (MGeoV v) = GVM.basicLength v + basicUnsafeSlice i n (MGeoV v) = MGeoV $ GVM.basicUnsafeSlice i n v + basicOverlaps (MGeoV v1) (MGeoV v2) = GVM.basicOverlaps v1 v2 + + basicUnsafeRead (MGeoV v) i = do (x,y) <- GVM.basicUnsafeRead v i + return $ Geo (Latitude x) $ Longitude y + basicUnsafeWrite (MGeoV v) i (Geo x y) = GVM.basicUnsafeWrite v i (unLat x, unLong y) + + basicClear (MGeoV v) = GVM.basicClear v + basicSet (MGeoV v) (Geo x y) = GVM.basicSet v (unLat x, unLong y) + + basicUnsafeNew n = MGeoV `liftM` GVM.basicUnsafeNew n + basicUnsafeReplicate n (Geo x y) = MGeoV `liftM` GVM.basicUnsafeReplicate n (unLat x, unLong y) + basicUnsafeCopy (MGeoV v1) (MGeoV v2) = GVM.basicUnsafeCopy v1 v2 + basicUnsafeGrow (MGeoV v) n = MGeoV `liftM` GVM.basicUnsafeGrow v n + +#if MIN_VERSION_vector(0,11,0) + basicInitialize (MGeoV v) = GVM.basicInitialize v +#endif + +instance GV.Vector Vector Geo where + {-# INLINE basicUnsafeFreeze #-} + {-# INLINE basicUnsafeThaw #-} + {-# INLINE basicLength #-} + {-# INLINE basicUnsafeSlice #-} + {-# INLINE basicUnsafeIndexM #-} + {-# INLINE elemseq #-} + basicUnsafeFreeze (MGeoV v) = GeoV `liftM` GV.basicUnsafeFreeze v + basicUnsafeThaw (GeoV v) = MGeoV `liftM` GV.basicUnsafeThaw v + basicLength (GeoV v) = GV.basicLength v + basicUnsafeSlice i n (GeoV v) = GeoV $ GV.basicUnsafeSlice i n v + basicUnsafeIndexM (GeoV v) i =do (x,y) <- GV.basicUnsafeIndexM v i + return $ Geo (Latitude x) $ Longitude y + + basicUnsafeCopy (MGeoV mv) (GeoV v) = GV.basicUnsafeCopy mv v + elemseq _ (Geo x y) = GV.elemseq (undefined :: Vector a) (unLat x, unLong y) diff --git a/Naqsha/Geometry/Coordinate/GeoHash.hs b/Naqsha/Geometry/Coordinate/GeoHash.hs new file mode 100644 --- /dev/null +++ b/Naqsha/Geometry/Coordinate/GeoHash.hs @@ -0,0 +1,214 @@ +-- | This module implements the geohash encoding of geo-locations. +-- https://en.wikipedia.org/wiki/Geohash. To try out geohash encoding +-- on web visit http://geohash.org + +module Naqsha.Geometry.Coordinate.GeoHash + ( GeoHash, encode, decode, accuracy, toByteString + ) where + + +import Data.Bits +import qualified Data.ByteString as B +import Data.ByteString.Internal ( c2w, w2c ) +import Data.Char ( ord ) +import Data.String +import Data.Monoid ( (<>) ) +import Data.Word ( Word8 ) + + + +import Naqsha.Geometry.Internal +import Naqsha.Geometry.Coordinate ( Geo(..) ) + + + +-- | Precision of encoding measured in base32 digits. +accuracyBase32 :: Int +accuracyBase32 = 12 + +-- | Precision of encoding measured in bits. +accuracy :: Int +accuracy = accuracyBase32 * 5 + +-- | The length of the output. +outputLength :: Int +outputLength = 2 * accuracyBase32 + +-- | The encoding of geo-coordinates as a geohash string. Currently, +-- the encoding supports 24 base32 digits of geo hash value which +-- means we loose about 4-bits of accuracy w.r.t the representation of +-- angles in the library. However, this loss is rather theoretical as +-- the angular error that results from such loss is so insignificant +-- that for all practical purposes, this accuracy is good enough --- +-- GPS devices will have much greater errors. The quantity `accuracy` +-- gives the number of bits of precision supported by the geohash +-- implementation exposed here. As expected GeoHash implementations +-- here will have problems at regions close to the poles. +newtype GeoHash = GeoHash B.ByteString deriving (Eq, Ord) + +instance Show GeoHash where + show (GeoHash x) = map b32ToChar $ B.unpack x + +instance IsString GeoHash where + fromString = GeoHash . B.pack . map cToB32 . take 24 + +------------------------------------------ Base 32 encoding used by geohash -------------------------- + +-- The digit ranges are +-- 0-9, b-h, jk, mn, p-z +-- +-- b - 10 +-- c - 11 +-- d - 12 +-- e - 13 +-- f - 14 +-- g - 15 +-- h - 16 +--------- Broken range --- +-- j - 17 +-- k - 18 +--------- Broken range ---- +-- m - 19 +-- n - 20 +---------- Broken range --- +-- p - 21 +-- q - 22 +-- r - 23 +-- s - 24 +-- t - 25 +-- u - 26 +-- v - 27 +-- w - 28 +-- x - 29 +-- y - 30 +-- x - 31 + +cToB32 :: Char -> Word8 +cToB32 x + | '0' <= x && x <= '9' = toEnum $ ord x - ord '0' + | 'b' <= x && x <= 'h' = toEnum $ ord x - ord 'b' + 10 + | 'p' <= x && x <= 'z' = toEnum $ ord x - ord 'p' + 21 + | x == 'j' = 17 + | x == 'k' = 18 + | x == 'm' = 19 + | x == 'n' = 20 + | otherwise = error $ "geohash: bad character " ++ show x + +b32ToChar8 :: Word8 -> Word8 +b32ToChar8 b32 + | 0 <= w && w <= 9 = c2w '0' + w + | 10 <= w && w <= 16 = c2w 'b' + w - 10 + | 21 <= w && w <= 32 = c2w 'p' + w - 21 + | w == 17 = c2w 'j' + | w == 18 = c2w 'k' + | w == 19 = c2w 'm' + | w == 20 = c2w 'n' + | otherwise = error "geohash: fatal this should never happen" + where w = b32 .&. 0x1F + +b32ToChar :: Word8 -> Char +b32ToChar = w2c . b32ToChar8 + +-- Geohash encoding +-- --------------- +-- +-- Notice that the bits of the geohash encoding is essentially got by +-- iterleaving the bits of the longitude and the latitude. However, +-- the first bit is 0 for negative angles 1 for positive +-- angles. Therefore we need to complement the sign bit. Since +-- longitudes vary over the entire range of angles, this is all we +-- need to do for adjustment before interleaving the bits. +-- +-- However, the latitudes like in the range -90 to +90. If we ignore the +-- +90 angle then we have the following property of its bit pattern +-- +-- 1. Every positive angle (other than +90) is of the form 00xxxxxxx. +-- +-- 2. Every negative angle is of the form 11xxxxxxx. +-- +-- Therefore, to get the actual bits that need to be interleaved, we +-- need to shift left the bits left by 1 and complement the +-- bit. During decoding, we need to do the reverse, i.e. complement +-- the bit and shift right by 1 with sign extension. +-- +-- For the +90 case while encoding we treat the bit stream as all 1's. +-- While decoding we will never get an angle of 90 but can be pretty +-- close. + +-- | Adjust the latitude for encoding. +adjustEncodeLat :: Latitude -> Angle +adjustEncodeLat lt + | testBit lt 63 = clearBit a 63 -- negative angle (starting bit = 0) + | testBit a 63 = complement zeroBits -- +90 + | otherwise = setBit a 63 -- positive angle (starting bit = 1) + where a = unsafeShiftL (unLat lt) 1 + +-- | Adjust the angle while decoding. +adjustDecodeLat :: Angle -> Latitude +adjustDecodeLat a = sgnBit .|. unsafeShiftR lt 1 + where lt = Latitude $ complementBit a 63 + sgnBit = bit 63 .&. lt + +-- | Adjusting longitude while encoding. Just nee +adjustEncodeLon :: Longitude -> Angle +adjustEncodeLon = flip complementBit 63 . unLong + +-- | Adjusting longitude while decoding +adjustDecodeLon :: Angle -> Longitude +adjustDecodeLon = Longitude . flip complementBit 63 + + +-- | Convert the geo hash to bytestring. +toByteString :: GeoHash -> B.ByteString +toByteString (GeoHash x) = B.map b32ToChar8 x + +--------------- Interleaved base32 encoding ------ + +-- | The @interleaveAndMerge (x,y)@ merges 5-bits, 3 from @x@ and 2 +-- from @y@ into a word and returns it. An appropriate swap is done so +-- that the next bytes are taken from y and x respectively. +interleaveAndMerge :: (Angle, Angle) -> (Word8, (Angle, Angle)) +interleaveAndMerge (x,y) = (w, (yp, xp)) + where xp = rotateL x 3 -- Take the top 3 bits + yp = rotateL y 2 -- Take the top 2 bits + wx = fromIntegral $ unAngle xp + wy = fromIntegral $ unAngle yp + w = unsafeShiftL (wx .&. 4) 2 -- x2 -> w4 + .|. unsafeShiftL (wx .&. 2) 1 -- x1 -> w2 + .|. (wx .&. 1) -- x0 -> w0 + .|. unsafeShiftL (wy .&. 2) 2 -- y1 -> w3 + .|. unsafeShiftL (wy .&. 1) 1 -- y0 -> w1 + + +-- | Encode a geo-location into its GeoHash string. +encode :: Geo -> GeoHash +encode (Geo lt lng) = GeoHash $ fst $ B.unfoldrN outputLength fld (adjustEncodeLon lng , adjustEncodeLat lt) + where fld = Just . interleaveAndMerge + +-------------------------- Decoding -------------------------------- + +-- | This function distributes the bits of the Word8 argument +-- (actually only 5-bits matter) to x and y in an interleaved fashion. +-- x gets 3-bits and y gets 2. The arguments are switched so that for +-- the next byte is distributed to y and x respectively. +splitAndDistribute :: (Angle, Angle) -> Word8 -> (Angle , Angle) +splitAndDistribute (x,y) w = (yp,xp) + where xp = unsafeShiftL x 3 + .|. (4 `bitTo` 2) + .|. (2 `bitTo` 1) + .|. (0 `bitTo` 0) + yp = unsafeShiftL y 2 + .|. (3 `bitTo` 1) + .|. (1 `bitTo` 0) + bitTo i j = Angle $ fromIntegral $ unsafeShiftL (unsafeShiftR w i .&. 1) j + + +-- | Decode the geo-location from its GeoHash string. +decode :: GeoHash -> Geo +decode (GeoHash hsh) = Geo lt ln + where lt = adjustDecodeLat $ unsafeShiftL y 4 + ln = adjustDecodeLon $ unsafeShiftL x 4 + (x,y) = B.foldl splitAndDistribute (Angle 0,Angle 0) strP + hshLen = B.length hsh + strP = if hshLen > outputLength then B.take outputLength hsh + else hsh <> B.replicate (outputLength - hshLen) 0 diff --git a/Naqsha/Geometry/Internal.hs b/Naqsha/Geometry/Internal.hs new file mode 100644 --- /dev/null +++ b/Naqsha/Geometry/Internal.hs @@ -0,0 +1,308 @@ +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE TypeFamilies #-} + +-- | The internal module that exposes the basic geometric types in +-- naqsha. This interface is subject to change and hence use with +-- caution. +module Naqsha.Geometry.Internal + ( Angle(..) + , degree , minute, second + , radian + , toDegree, toRadian + , Latitude(..), Longitude(..), lat, lon + ) where + +import Control.Monad ( liftM ) +import Data.Bits ( Bits ) +import Data.Fixed +import Data.Group +import Data.Int +import Data.Monoid +import GHC.Real +import Data.Vector.Unboxed ( MVector(..), Vector, Unbox) +import qualified Data.Vector.Generic as GV +import qualified Data.Vector.Generic.Mutable as GVM +import Text.Read + +----------------------------- Angles and Angular quantities ----------------------- + +-- | An abstract angle. Internally, angles are represented as a 64-bit +-- integer with each unit contribute 1/2^64 fraction of a complete +-- circle. This means that angles are accurate up to a resolution of 2 +-- π / 2^64 radians. Angles form a group under the angular addition +-- and the fact that these are represented as integers means one can +-- expect high speed accurate angle arithmetic. +-- +-- When expressing angles one can use a more convenient notation: +-- +-- > myAngle = degree 21.71167 +-- > yourAngle = degree 21 <> minute 42 <> second 42 +-- +newtype Angle = Angle {unAngle :: Int64} deriving (Enum, Eq, Ord, Unbox, Show, Read, Bits) + +instance Monoid Angle where + mempty = Angle 0 + mappend (Angle x) (Angle y) = Angle $ x + y + mconcat = Angle . sum . map unAngle + +instance Group Angle where + invert (Angle x) = Angle $ negate x + +instance Bounded Angle where + maxBound = Angle maxBound + minBound = Angle minBound + +-- | Express angle in degrees. +degree :: Rational -> Angle +degree = Angle . fromInteger . round . (*scale) + where scale = (2^(64:: Int)) % 360 + +-- | Express angle in minutes. +minute :: Rational -> Angle +minute = degree . (*scale) + where scale = 1 % 60 + +-- | Express angle in seconds. +second :: Rational -> Angle +second = degree . (*scale) + where scale = 1 % 3600 + +-- | Express angle in radians +radian :: Double -> Angle +radian = Angle . round . (*scale) + where scale = (2^(63:: Int)) / pi + + +---------------------- Decimal representation of angle ---------------------------------- + +-- | Measure angle in degrees. This conversion may lead to loss of +-- precision. +toDegree :: Fractional r => Angle -> r +toDegree = fromRational . (*conv) . toRational . unAngle + where conv = 360 % (2^(64 :: Int)) + +-- | Measure angle in radians. This conversion may lead to loss of +-- precision. +toRadian :: Angle -> Double +toRadian = (*conv) . fromIntegral . unAngle + where conv = pi / (2^(63:: Int)) + +------------------- Making stuff suitable for unboxed vector. -------------------------- + +newtype instance MVector s Angle = MAngV (MVector s Int64) +newtype instance Vector Angle = AngV (Vector Int64) + + +instance GVM.MVector MVector Angle where + {-# INLINE basicLength #-} + {-# INLINE basicUnsafeSlice #-} + {-# INLINE basicOverlaps #-} + {-# INLINE basicUnsafeNew #-} + {-# INLINE basicUnsafeReplicate #-} + {-# INLINE basicUnsafeRead #-} + {-# INLINE basicUnsafeWrite #-} + {-# INLINE basicClear #-} + {-# INLINE basicSet #-} + {-# INLINE basicUnsafeCopy #-} + {-# INLINE basicUnsafeGrow #-} + basicLength (MAngV v) = GVM.basicLength v + basicUnsafeSlice i n (MAngV v) = MAngV $ GVM.basicUnsafeSlice i n v + basicOverlaps (MAngV v1) (MAngV v2) = GVM.basicOverlaps v1 v2 + + basicUnsafeRead (MAngV v) i = Angle `liftM` GVM.basicUnsafeRead v i + basicUnsafeWrite (MAngV v) i (Angle x) = GVM.basicUnsafeWrite v i x + + basicClear (MAngV v) = GVM.basicClear v + basicSet (MAngV v) (Angle x) = GVM.basicSet v x + + basicUnsafeNew n = MAngV `liftM` GVM.basicUnsafeNew n + basicUnsafeReplicate n (Angle x) = MAngV `liftM` GVM.basicUnsafeReplicate n x + basicUnsafeCopy (MAngV v1) (MAngV v2) = GVM.basicUnsafeCopy v1 v2 + basicUnsafeGrow (MAngV v) n = MAngV `liftM` GVM.basicUnsafeGrow v n + +#if MIN_VERSION_vector(0,11,0) + basicInitialize (MAngV v) = GVM.basicInitialize v +#endif + +instance GV.Vector Vector Angle where + {-# INLINE basicUnsafeFreeze #-} + {-# INLINE basicUnsafeThaw #-} + {-# INLINE basicLength #-} + {-# INLINE basicUnsafeSlice #-} + {-# INLINE basicUnsafeIndexM #-} + {-# INLINE elemseq #-} + basicUnsafeFreeze (MAngV v) = AngV `liftM` GV.basicUnsafeFreeze v + basicUnsafeThaw (AngV v) = MAngV `liftM` GV.basicUnsafeThaw v + basicLength (AngV v) = GV.basicLength v + basicUnsafeSlice i n (AngV v) = AngV $ GV.basicUnsafeSlice i n v + basicUnsafeIndexM (AngV v) i = Angle `liftM` GV.basicUnsafeIndexM v i + + basicUnsafeCopy (MAngV mv) (AngV v) = GV.basicUnsafeCopy mv v + elemseq _ (Angle x) = GV.elemseq (undefined :: Vector a) x + +------------------------------------- Latitude and Longitude --------------------------------- + + +-- | The latitude of a point. Positive denotes North of Equator where +-- as negative South. +newtype Latitude = Latitude { unLat :: Angle } deriving (Eq, Ord, Bits) + + +instance Show Latitude where + show = show . (toDegree :: Angle -> Nano) . unLat + +instance Read Latitude where + readPrec = conv <$> readPrec + where conv = lat . degree . (toRational :: Nano -> Rational) + +instance Bounded Latitude where + maxBound = lat $ degree 90 + minBound = lat $ degree (-90) + + +-- | Construct latitude out of an angle. +lat :: Angle -> Latitude +lat = Latitude . normLat + + +-- | normalise latitude values. +normLat :: Angle -> Angle +normLat ang | degree (-90) <= ang && ang < degree 90 = ang + | ang > degree 90 = succ (maxBound <> invert ang) + | otherwise = minBound <> invert ang + + +-------------------------- Longitude ------------------------------------------ + +-- | The longitude of a point. Positive denotes East of the Greenwich +-- meridian where as negative denotes West. +newtype Longitude = Longitude { unLong :: Angle } + deriving (Eq, Bounded, Ord, Monoid, Group, Bits) + +-- | Convert angles to longitude. +lon :: Angle -> Longitude +lon = Longitude + + +instance Show Longitude where + show = show . (toDegree :: Angle -> Nano) . unLong + +instance Read Longitude where + readPrec = conv <$> readPrec + where conv = lon . degree . (toRational :: Nano -> Rational) + + +--------------------------- Internal helper functions ------------------------ + + +newtype instance MVector s Latitude = MLatV (MVector s Angle) +newtype instance Vector Latitude = LatV (Vector Angle) + + +newtype instance MVector s Longitude = MLongV (MVector s Angle) +newtype instance Vector Longitude = LongV (Vector Angle) + + +-------------------- Instance for Angle -------------------------------------------- + + +-------------------- Instance for latitude -------------------------------------------- + +instance GVM.MVector MVector Latitude where + {-# INLINE basicLength #-} + {-# INLINE basicUnsafeSlice #-} + {-# INLINE basicOverlaps #-} + {-# INLINE basicUnsafeNew #-} + {-# INLINE basicUnsafeReplicate #-} + {-# INLINE basicUnsafeRead #-} + {-# INLINE basicUnsafeWrite #-} + {-# INLINE basicClear #-} + {-# INLINE basicSet #-} + {-# INLINE basicUnsafeCopy #-} + {-# INLINE basicUnsafeGrow #-} + basicLength (MLatV v) = GVM.basicLength v + basicUnsafeSlice i n (MLatV v) = MLatV $ GVM.basicUnsafeSlice i n v + basicOverlaps (MLatV v1) (MLatV v2) = GVM.basicOverlaps v1 v2 + + basicUnsafeRead (MLatV v) i = Latitude `liftM` GVM.basicUnsafeRead v i + basicUnsafeWrite (MLatV v) i (Latitude x) = GVM.basicUnsafeWrite v i x + + basicClear (MLatV v) = GVM.basicClear v + basicSet (MLatV v) (Latitude x) = GVM.basicSet v x + + basicUnsafeNew n = MLatV `liftM` GVM.basicUnsafeNew n + basicUnsafeReplicate n (Latitude x) = MLatV `liftM` GVM.basicUnsafeReplicate n x + basicUnsafeCopy (MLatV v1) (MLatV v2) = GVM.basicUnsafeCopy v1 v2 + basicUnsafeGrow (MLatV v) n = MLatV `liftM` GVM.basicUnsafeGrow v n + +#if MIN_VERSION_vector(0,11,0) + basicInitialize (MLatV v) = GVM.basicInitialize v +#endif + +instance GV.Vector Vector Latitude where + {-# INLINE basicUnsafeFreeze #-} + {-# INLINE basicUnsafeThaw #-} + {-# INLINE basicLength #-} + {-# INLINE basicUnsafeSlice #-} + {-# INLINE basicUnsafeIndexM #-} + {-# INLINE elemseq #-} + basicUnsafeFreeze (MLatV v) = LatV `liftM` GV.basicUnsafeFreeze v + basicUnsafeThaw (LatV v) = MLatV `liftM` GV.basicUnsafeThaw v + basicLength (LatV v) = GV.basicLength v + basicUnsafeSlice i n (LatV v) = LatV $ GV.basicUnsafeSlice i n v + basicUnsafeIndexM (LatV v) i = Latitude `liftM` GV.basicUnsafeIndexM v i + + basicUnsafeCopy (MLatV mv) (LatV v) = GV.basicUnsafeCopy mv v + elemseq _ (Latitude x) = GV.elemseq (undefined :: Vector a) x + + +-------------------------------- Instance for Longitude ----------------------------------- + +instance GVM.MVector MVector Longitude where + {-# INLINE basicLength #-} + {-# INLINE basicUnsafeSlice #-} + {-# INLINE basicOverlaps #-} + {-# INLINE basicUnsafeNew #-} + {-# INLINE basicUnsafeReplicate #-} + {-# INLINE basicUnsafeRead #-} + {-# INLINE basicUnsafeWrite #-} + {-# INLINE basicClear #-} + {-# INLINE basicSet #-} + {-# INLINE basicUnsafeCopy #-} + {-# INLINE basicUnsafeGrow #-} + basicLength (MLongV v) = GVM.basicLength v + basicUnsafeSlice i n (MLongV v) = MLongV $ GVM.basicUnsafeSlice i n v + basicOverlaps (MLongV v1) (MLongV v2) = GVM.basicOverlaps v1 v2 + + basicUnsafeRead (MLongV v) i = Longitude `liftM` GVM.basicUnsafeRead v i + basicUnsafeWrite (MLongV v) i (Longitude x) = GVM.basicUnsafeWrite v i x + + basicClear (MLongV v) = GVM.basicClear v + basicSet (MLongV v) (Longitude x) = GVM.basicSet v x + + basicUnsafeNew n = MLongV `liftM` GVM.basicUnsafeNew n + basicUnsafeReplicate n (Longitude x) = MLongV `liftM` GVM.basicUnsafeReplicate n x + basicUnsafeCopy (MLongV v1) (MLongV v2) = GVM.basicUnsafeCopy v1 v2 + basicUnsafeGrow (MLongV v) n = MLongV `liftM` GVM.basicUnsafeGrow v n + +#if MIN_VERSION_vector(0,11,0) + basicInitialize (MLongV v) = GVM.basicInitialize v +#endif + +instance GV.Vector Vector Longitude where + {-# INLINE basicUnsafeFreeze #-} + {-# INLINE basicUnsafeThaw #-} + {-# INLINE basicLength #-} + {-# INLINE basicUnsafeSlice #-} + {-# INLINE basicUnsafeIndexM #-} + {-# INLINE elemseq #-} + basicUnsafeFreeze (MLongV v) = LongV `liftM` GV.basicUnsafeFreeze v + basicUnsafeThaw (LongV v) = MLongV `liftM` GV.basicUnsafeThaw v + basicLength (LongV v) = GV.basicLength v + basicUnsafeSlice i n (LongV v) = LongV $ GV.basicUnsafeSlice i n v + basicUnsafeIndexM (LongV v) i = Longitude `liftM` GV.basicUnsafeIndexM v i + + basicUnsafeCopy (MLongV mv) (LongV v) = GV.basicUnsafeCopy mv v + elemseq _ (Longitude x) = GV.elemseq (undefined :: Vector a) x diff --git a/Naqsha/Geometry/Spherical.hs b/Naqsha/Geometry/Spherical.hs --- a/Naqsha/Geometry/Spherical.hs +++ b/Naqsha/Geometry/Spherical.hs @@ -1,5 +1,9 @@ -- | Geometric operations on earth surface assuming that earth is a -- sphere of radius 6371008 m. +-- +-- TODO: Port some calculations from +-- http://www.movable-type.co.uk/scripts/latlong.html +-- module Naqsha.Geometry.Spherical ( -- * Distance calculation. @@ -9,7 +13,7 @@ import Data.Monoid import Data.Group -import Naqsha.Position +import Naqsha.Geometry.Coordinate import Naqsha.Geometry.Angle --------------------- Distance calculation ------------------------------------- @@ -29,7 +33,7 @@ distance = distance' rMean --- | A generalisation of `dHvS` that takes the radius as +-- | A generalisation of `distance` that takes the radius as -- argument. Will work on Mars for example once we set up a latitude -- longitude system there. For this function units does not matter --- -- the computed distance is in the same unit as the input radius. We diff --git a/Naqsha/Position.hs b/Naqsha/Position.hs deleted file mode 100644 --- a/Naqsha/Position.hs +++ /dev/null @@ -1,361 +0,0 @@ -{-# LANGUAGE CPP #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE Rank2Types #-} --- | This module captures position of a point on the globe. -module Naqsha.Position - ( -- * Basics - -- $latandlong$ - Geo(..) - , northPole, southPole - -- ** Latitudes - , Latitude - , north, south, lat - , equator - , tropicOfCancer - , tropicOfCapricon - -- ** Longitudes. - , Longitude - , east, west, lon - , greenwich - ) where - -import Control.Monad ( liftM ) -import Data.Default -import Data.Fixed -import Data.Monoid -import Data.Group -import Data.Vector.Unboxed ( MVector(..), Vector) -import qualified Data.Vector.Generic as GV -import qualified Data.Vector.Generic.Mutable as GVM -import Text.Read - -import Prelude -- To avoid redundunt import warnings. - -import Naqsha.Geometry.Angle - --- $latandlong$ --- --- A point on the globe is specified by giving its geo coordinates --- captures by the type `Geo`. It is essentially a pair of the --- `Latitude` and `Longitude` of the point. --- --- == Examples --- --- > kanpurLatitude :: Latitude --- > kanpurLatitude = lat $ degree 26.4477777 --- > kanpurLongitude :: Longitude --- > kanpurLongitude = lon $ degree 80.3461111 --- --- --- > kanpurLatitude = lat $ degree 26 <> minute 26 <> second 52 --- > kanpurLongitude = lon $ degree 80 <> minute 20 <> second 46 --- --- The show and read instance of the `Latitude` and `Longitude` types --- uses degrees for displaying and reading respectively. Show and Read --- instances can express these quantities up to Nano degree precision. --- --- == Convention on sign. --- --- For latitudes, positive means north and negative means south. For --- longitudes, positive means east and negative means west. However, --- if you find these conventions confusing you can use the combinators --- `north`, `south`, `east`, and `west` when constructing latitudes or --- longitudes. --- - - ------------------------------ Lattitude ---------------------------------- - --- | The latitude of a point. Positive denotes North of Equator where --- as negative South. -newtype Latitude = Latitude { unLat :: Angle } deriving (Eq, Ord) - --- | Construct latitude out of an angle. -lat :: Angle -> Latitude -lat = Latitude . normLat - --- | Convert an angle to a northern latitude --- --- > tropicOfCancer = north $ degree 23.5 --- -north :: Angle -> Latitude -north = lat - --- | Convert an angle to a southern latitude. --- --- > tropicOfCapricon = south $ degree 23.5 --- -south :: Angle -> Latitude -south = lat . invert - - -instance Angular Latitude where - toAngle = unLat - -instance Show Latitude where - show = show . (toDegree :: Angle -> Nano) . unLat - -instance Read Latitude where - readPrec = conv <$> readPrec - where conv = lat . degree . (toRational :: Nano -> Rational) - - -instance Default Latitude where - def = equator - --- | The latitude of equator. -equator :: Latitude -equator = lat $ degree 0 - --- | The latitude corresponding to the Tropic of Cancer. -tropicOfCancer :: Latitude -tropicOfCancer = north $ degree 23.5 - --- | The latitude corresponding to the Tropic of Capricon -tropicOfCapricon :: Latitude -tropicOfCapricon = south $ degree 23.5 - - -instance Bounded Latitude where - maxBound = lat $ degree 90 - minBound = lat $ degree (-90) - - --------------------------- Longitude ------------------------------------------ - --- | The longitude of a point. Positive denotes East of the Greenwich --- meridian where as negative denotes West. -newtype Longitude = Longitude { unLong :: Angle } - deriving (Eq, Bounded, Default, Angular, Ord, Monoid, Group) - -instance Show Longitude where - show = show . (toDegree :: Angle -> Nano) . unLong - -instance Read Longitude where - readPrec = conv <$> readPrec - where conv = lon . degree . (toRational :: Nano -> Rational) - --- | Convert angles to longitude. -lon :: Angle -> Longitude -lon = Longitude - --- | Convert angle to an eastern longitude. --- --- > kanpurLongitude = east $ degree 80.3461 --- -east :: Angle -> Longitude -east = lon - --- | Convert angle to a western longitude --- --- > newyorkLongitude = west $ degree 74.0059 --- -west :: Angle -> Longitude -west = lon . invert - - - --- | The zero longitude. -greenwich :: Longitude -greenwich = lon $ degree 0 - --- | The coordinates of a point on the earth's surface. -data Geo = Geo {-# UNPACK #-} !Latitude - {-# UNPACK #-} !Longitude - deriving Show - -instance Default Geo where - def = Geo def def - - --- | The North pole -northPole :: Geo -northPole = Geo maxBound $ lon $ degree 0 - --- | The South pole -southPole :: Geo -southPole = Geo minBound $ lon $ degree 0 - -instance Eq Geo where - (==) (Geo xlat xlong) (Geo ylat ylong) - | xlat == maxBound = ylat == maxBound -- longitude irrelevant for north pole - | xlat == minBound = ylat == minBound -- longitude irrelevant for south pole - | otherwise = xlat == ylat && xlong == ylong - --- | normalise latitude values. -normLat :: Angle -> Angle -normLat ang | degree (-90) <= ang && ang < degree 90 = ang - | ang > degree 90 = succ (maxBound <> invert ang) - | otherwise = minBound <> invert ang - - ---------------------------- Internal helper functions ------------------------ - - -newtype instance MVector s Latitude = MLatV (MVector s Angle) -newtype instance Vector Latitude = LatV (Vector Angle) - - -newtype instance MVector s Longitude = MLongV (MVector s Angle) -newtype instance Vector Longitude = LongV (Vector Angle) - - -newtype instance MVector s Geo = MGeoV (MVector s (Angle,Angle)) -newtype instance Vector Geo = GeoV (Vector (Angle,Angle)) - - --------------------- Instance for Angle -------------------------------------------- - - --------------------- Instance for latitude -------------------------------------------- - -instance GVM.MVector MVector Latitude where - {-# INLINE basicLength #-} - {-# INLINE basicUnsafeSlice #-} - {-# INLINE basicOverlaps #-} - {-# INLINE basicUnsafeNew #-} - {-# INLINE basicUnsafeReplicate #-} - {-# INLINE basicUnsafeRead #-} - {-# INLINE basicUnsafeWrite #-} - {-# INLINE basicClear #-} - {-# INLINE basicSet #-} - {-# INLINE basicUnsafeCopy #-} - {-# INLINE basicUnsafeGrow #-} - basicLength (MLatV v) = GVM.basicLength v - basicUnsafeSlice i n (MLatV v) = MLatV $ GVM.basicUnsafeSlice i n v - basicOverlaps (MLatV v1) (MLatV v2) = GVM.basicOverlaps v1 v2 - - basicUnsafeRead (MLatV v) i = Latitude `liftM` GVM.basicUnsafeRead v i - basicUnsafeWrite (MLatV v) i (Latitude x) = GVM.basicUnsafeWrite v i x - - basicClear (MLatV v) = GVM.basicClear v - basicSet (MLatV v) (Latitude x) = GVM.basicSet v x - - basicUnsafeNew n = MLatV `liftM` GVM.basicUnsafeNew n - basicUnsafeReplicate n (Latitude x) = MLatV `liftM` GVM.basicUnsafeReplicate n x - basicUnsafeCopy (MLatV v1) (MLatV v2) = GVM.basicUnsafeCopy v1 v2 - basicUnsafeGrow (MLatV v) n = MLatV `liftM` GVM.basicUnsafeGrow v n - -#if MIN_VERSION_vector(0,11,0) - basicInitialize (MLatV v) = GVM.basicInitialize v -#endif - -instance GV.Vector Vector Latitude where - {-# INLINE basicUnsafeFreeze #-} - {-# INLINE basicUnsafeThaw #-} - {-# INLINE basicLength #-} - {-# INLINE basicUnsafeSlice #-} - {-# INLINE basicUnsafeIndexM #-} - {-# INLINE elemseq #-} - basicUnsafeFreeze (MLatV v) = LatV `liftM` GV.basicUnsafeFreeze v - basicUnsafeThaw (LatV v) = MLatV `liftM` GV.basicUnsafeThaw v - basicLength (LatV v) = GV.basicLength v - basicUnsafeSlice i n (LatV v) = LatV $ GV.basicUnsafeSlice i n v - basicUnsafeIndexM (LatV v) i = Latitude `liftM` GV.basicUnsafeIndexM v i - - basicUnsafeCopy (MLatV mv) (LatV v) = GV.basicUnsafeCopy mv v - elemseq _ (Latitude x) = GV.elemseq (undefined :: Vector a) x - - --------------------------------- Instance for Longitude ----------------------------------- - -instance GVM.MVector MVector Longitude where - {-# INLINE basicLength #-} - {-# INLINE basicUnsafeSlice #-} - {-# INLINE basicOverlaps #-} - {-# INLINE basicUnsafeNew #-} - {-# INLINE basicUnsafeReplicate #-} - {-# INLINE basicUnsafeRead #-} - {-# INLINE basicUnsafeWrite #-} - {-# INLINE basicClear #-} - {-# INLINE basicSet #-} - {-# INLINE basicUnsafeCopy #-} - {-# INLINE basicUnsafeGrow #-} - basicLength (MLongV v) = GVM.basicLength v - basicUnsafeSlice i n (MLongV v) = MLongV $ GVM.basicUnsafeSlice i n v - basicOverlaps (MLongV v1) (MLongV v2) = GVM.basicOverlaps v1 v2 - - basicUnsafeRead (MLongV v) i = Longitude `liftM` GVM.basicUnsafeRead v i - basicUnsafeWrite (MLongV v) i (Longitude x) = GVM.basicUnsafeWrite v i x - - basicClear (MLongV v) = GVM.basicClear v - basicSet (MLongV v) (Longitude x) = GVM.basicSet v x - - basicUnsafeNew n = MLongV `liftM` GVM.basicUnsafeNew n - basicUnsafeReplicate n (Longitude x) = MLongV `liftM` GVM.basicUnsafeReplicate n x - basicUnsafeCopy (MLongV v1) (MLongV v2) = GVM.basicUnsafeCopy v1 v2 - basicUnsafeGrow (MLongV v) n = MLongV `liftM` GVM.basicUnsafeGrow v n - -#if MIN_VERSION_vector(0,11,0) - basicInitialize (MLongV v) = GVM.basicInitialize v -#endif - -instance GV.Vector Vector Longitude where - {-# INLINE basicUnsafeFreeze #-} - {-# INLINE basicUnsafeThaw #-} - {-# INLINE basicLength #-} - {-# INLINE basicUnsafeSlice #-} - {-# INLINE basicUnsafeIndexM #-} - {-# INLINE elemseq #-} - basicUnsafeFreeze (MLongV v) = LongV `liftM` GV.basicUnsafeFreeze v - basicUnsafeThaw (LongV v) = MLongV `liftM` GV.basicUnsafeThaw v - basicLength (LongV v) = GV.basicLength v - basicUnsafeSlice i n (LongV v) = LongV $ GV.basicUnsafeSlice i n v - basicUnsafeIndexM (LongV v) i = Longitude `liftM` GV.basicUnsafeIndexM v i - - basicUnsafeCopy (MLongV mv) (LongV v) = GV.basicUnsafeCopy mv v - elemseq _ (Longitude x) = GV.elemseq (undefined :: Vector a) x - - ------------------------------ Instance for Geo --------------------------------------------- - -instance GVM.MVector MVector Geo where - {-# INLINE basicLength #-} - {-# INLINE basicUnsafeSlice #-} - {-# INLINE basicOverlaps #-} - {-# INLINE basicUnsafeNew #-} - {-# INLINE basicUnsafeReplicate #-} - {-# INLINE basicUnsafeRead #-} - {-# INLINE basicUnsafeWrite #-} - {-# INLINE basicClear #-} - {-# INLINE basicSet #-} - {-# INLINE basicUnsafeCopy #-} - {-# INLINE basicUnsafeGrow #-} - basicLength (MGeoV v) = GVM.basicLength v - basicUnsafeSlice i n (MGeoV v) = MGeoV $ GVM.basicUnsafeSlice i n v - basicOverlaps (MGeoV v1) (MGeoV v2) = GVM.basicOverlaps v1 v2 - - basicUnsafeRead (MGeoV v) i = do (x,y) <- GVM.basicUnsafeRead v i - return $ Geo (Latitude x) $ Longitude y - basicUnsafeWrite (MGeoV v) i (Geo x y) = GVM.basicUnsafeWrite v i (unLat x, unLong y) - - basicClear (MGeoV v) = GVM.basicClear v - basicSet (MGeoV v) (Geo x y) = GVM.basicSet v (unLat x, unLong y) - - basicUnsafeNew n = MGeoV `liftM` GVM.basicUnsafeNew n - basicUnsafeReplicate n (Geo x y) = MGeoV `liftM` GVM.basicUnsafeReplicate n (unLat x, unLong y) - basicUnsafeCopy (MGeoV v1) (MGeoV v2) = GVM.basicUnsafeCopy v1 v2 - basicUnsafeGrow (MGeoV v) n = MGeoV `liftM` GVM.basicUnsafeGrow v n - -#if MIN_VERSION_vector(0,11,0) - basicInitialize (MGeoV v) = GVM.basicInitialize v -#endif - -instance GV.Vector Vector Geo where - {-# INLINE basicUnsafeFreeze #-} - {-# INLINE basicUnsafeThaw #-} - {-# INLINE basicLength #-} - {-# INLINE basicUnsafeSlice #-} - {-# INLINE basicUnsafeIndexM #-} - {-# INLINE elemseq #-} - basicUnsafeFreeze (MGeoV v) = GeoV `liftM` GV.basicUnsafeFreeze v - basicUnsafeThaw (GeoV v) = MGeoV `liftM` GV.basicUnsafeThaw v - basicLength (GeoV v) = GV.basicLength v - basicUnsafeSlice i n (GeoV v) = GeoV $ GV.basicUnsafeSlice i n v - basicUnsafeIndexM (GeoV v) i =do (x,y) <- GV.basicUnsafeIndexM v i - return $ Geo (Latitude x) $ Longitude y - - basicUnsafeCopy (MGeoV mv) (GeoV v) = GV.basicUnsafeCopy mv v - elemseq _ (Geo x y) = GV.elemseq (undefined :: Vector a) (unLat x, unLong y) diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,34 @@ +Naqsha +------ + +[![Build Staus][travis-status]][travis-naqsha] +[![Hackage][hackage-badge]][hackage] +[![Hackage Dependencies][hackage-deps-badge]][hackage-deps] +[![Stackage LTS][stackage-lts-badge]][stackage-lts] +[![Stackage Nightly][stackage-nightly-badge]][stackage-nightly] + + +Naqsha is a Haskell library to work with geospatial data types. The +goal of this library is to provide fast and high level access to +various operations in the geospatial setting like distance +calculations, azimuth etc. Due to the very nature of the geometry on +the surface of the globe, many of these basic tasks have multiple +algorithms based on criteria like speed of computation and +accuracy. We would like to provide all such algorithms in this +package. + +The word naqsha, or more accurately naqshA (नक़्शा), means a map or a +sketch. + +[travis-status]: "Build status" +[travis-naqsha]: + +[hackage]: +[hackage-badge]: +[hackage-deps-badge]: +[hackage-deps]: + +[stackage-lts]: +[stackage-nightly]: +[stackage-lts-badge]: +[stackage-nightly-badge]: diff --git a/naqsha.cabal b/naqsha.cabal --- a/naqsha.cabal +++ b/naqsha.cabal @@ -1,5 +1,5 @@ name: naqsha -version: 0.1.0.0 +version: 0.2.0.0 synopsis: A library for working with geospatial data types. description: Naqsha is a library to work with geospatial data types like latitudes and longitudes. It provides @@ -18,7 +18,7 @@ build-type: Simple cabal-version: >=1.10 extra-source-files: CHANGELOG.md - + , README.md bug-reports: https://github.com/naqsha/naqsha/issues @@ -29,16 +29,19 @@ library ghc-options: -Wall build-depends: base >= 4.6 && < 4.11 - , data-default + , bytestring >= 0.9 && < 0.11 , groups , vector >= 0.7.1 && < 0.13 - exposed-modules: Naqsha - , Naqsha.Position - , Naqsha.Geometry.Angle + exposed-modules: Naqsha.Geometry + , Naqsha.Geometry.Coordinate.GeoHash + , Naqsha.Geometry.Internal , Naqsha.Geometry.Spherical , Naqsha.Version other-modules: Paths_naqsha + , Naqsha.Geometry.Angle + , Naqsha.Geometry.Coordinate + default-language: Haskell2010 test-Suite test @@ -47,9 +50,11 @@ hs-source-dirs: tests main-is: Main.hs ghc-options: -Wall - other-modules: Naqsha.PositionSpec + other-modules: Naqsha.Arbitrary , Naqsha.Geometry.AngleSpec - , Naqsha.Arbitrary + , Naqsha.Geometry.CoordinateSpec + , Naqsha.Geometry.Coordinate.GeoHashSpec + build-depends: base , HUnit >= 1.2 , QuickCheck >= 2.4 @@ -58,4 +63,4 @@ -- -- This package -- - , naqsha == 0.1.0.0 + , naqsha diff --git a/tests/Naqsha/Arbitrary.hs b/tests/Naqsha/Arbitrary.hs --- a/tests/Naqsha/Arbitrary.hs +++ b/tests/Naqsha/Arbitrary.hs @@ -11,7 +11,8 @@ import Test.QuickCheck -import Naqsha +import Naqsha.Geometry +import Naqsha.Geometry.Coordinate.GeoHash instance Arbitrary Angle where @@ -26,3 +27,10 @@ instance Arbitrary Geo where arbitrary = Geo <$> arbitrary <*> arbitrary + + +geoHashRange :: Gen Geo +geoHashRange = suchThat arbitrary ( \ g -> g /= northPole && g /= southPole) + +instance Arbitrary GeoHash where + arbitrary = encode <$> geoHashRange diff --git a/tests/Naqsha/Geometry/AngleSpec.hs b/tests/Naqsha/Geometry/AngleSpec.hs --- a/tests/Naqsha/Geometry/AngleSpec.hs +++ b/tests/Naqsha/Geometry/AngleSpec.hs @@ -6,7 +6,7 @@ import Test.Hspec import Test.Hspec.QuickCheck -import Naqsha +import Naqsha.Geometry import Naqsha.Arbitrary () spec :: Spec diff --git a/tests/Naqsha/Geometry/Coordinate/GeoHashSpec.hs b/tests/Naqsha/Geometry/Coordinate/GeoHashSpec.hs new file mode 100644 --- /dev/null +++ b/tests/Naqsha/Geometry/Coordinate/GeoHashSpec.hs @@ -0,0 +1,27 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE Rank2Types #-} +module Naqsha.Geometry.Coordinate.GeoHashSpec where + +import Data.Bits +import Data.Group +import Data.Monoid +import Data.String + +import Test.Hspec +import Test.Hspec.QuickCheck + +import Naqsha.Geometry +import Naqsha.Geometry.Coordinate.GeoHash as GeoHash +import Naqsha.Arbitrary() + +approxEq :: Geo -> Geo -> Bool +approxEq (Geo x1 y1) (Geo x2 y2) = abs dx <= err && abs dy <= err + where dx = fromEnum $ toAngle x1 <> invert (toAngle x2) + dy = fromEnum $ toAngle y1 <> invert (toAngle y2) + err = bit $ 64 - accuracy + +spec :: Spec +spec = do + prop "fromString . show = id" $ \ (g :: GeoHash) -> (fromString $ show g) `shouldBe` g + prop "encode . decode = id" $ \ (g :: GeoHash) -> (encode $ decode g) `shouldBe` g + prop "decode . encode = id (approx)" $ \ (g :: Geo) -> (decode $ encode g) `approxEq` g diff --git a/tests/Naqsha/Geometry/CoordinateSpec.hs b/tests/Naqsha/Geometry/CoordinateSpec.hs new file mode 100644 --- /dev/null +++ b/tests/Naqsha/Geometry/CoordinateSpec.hs @@ -0,0 +1,69 @@ +{-# LANGUAGE ScopedTypeVariables #-} +module Naqsha.Geometry.CoordinateSpec where +import Data.Monoid +import Data.Fixed +import Test.QuickCheck +import Test.Hspec +import Test.Hspec.QuickCheck + +import Naqsha.Geometry +import Naqsha.Arbitrary () + +inRange :: Testable prop + => (Angle, Angle) -- ^ Range + -> String -- ^ description + -> ((Angle,Angle) -> prop) + -> Spec +inRange (mi,mx) descr prop_test = it msg $ forAll pair prop_test + where pair = (,) <$> gen <*> gen + msg = "in range " ++ show (toNano mi,toNano mx) ++ ": " ++ descr + gen = toEnum <$> choose (fromEnum mi, fromEnum mx) + + toNano :: Angle -> Nano + toNano = toDegree + +isIncreasing :: (Angle, Angle) -> Bool +isIncreasing (x , y) + | x == y = xA == yA + | x < y = xA < yA + | otherwise = xA > yA + where xA = lat x + yA = lat y + + +isDecreasing :: (Angle, Angle) -> Bool +isDecreasing (x, y) + | x == y = xA == yA + | x > y = xA < yA + | otherwise = xA > yA + where xA = lat x + yA = lat y + +shouldBeBounded :: (Arbitrary a, Ord a, Show a, Bounded a) => a -> Spec +shouldBeBounded a = prop msg $ \ x -> x >= mi && x <= mx + where msg = unwords [ "should lie between" + , show mi + , "and" + , show mx + ] + mi = minBound `asTypeOf` a + mx = maxBound `asTypeOf` a + + +spec :: Spec +spec = do + + describe "latitudes" $ do + + inRange (degree (-90), degree 90) "increases monotonically" isIncreasing + inRange (degree 90 , maxBound) "decreases monotonically" isDecreasing + inRange (minBound , degree (-90)) "decreases monotonically" isDecreasing + + shouldBeBounded (undefined :: Latitude) + + describe "longitudes" $ do + + prop "should have a period of 360 deg" $ + \ (x :: Longitude) -> x <> lon (degree 360) `shouldBe` x + + shouldBeBounded (undefined ::Longitude) diff --git a/tests/Naqsha/PositionSpec.hs b/tests/Naqsha/PositionSpec.hs deleted file mode 100644 --- a/tests/Naqsha/PositionSpec.hs +++ /dev/null @@ -1,70 +0,0 @@ -{-# LANGUAGE ScopedTypeVariables #-} -module Naqsha.PositionSpec where -import Data.Monoid -import Data.Fixed -import Test.QuickCheck -import Test.Hspec -import Test.Hspec.QuickCheck - -import Naqsha.Position -import Naqsha.Geometry.Angle -import Naqsha.Arbitrary () - -inRange :: Testable prop - => (Angle, Angle) -- ^ Range - -> String -- ^ description - -> ((Angle,Angle) -> prop) - -> Spec -inRange (mi,mx) descr prop_test = it msg $ forAll pair prop_test - where pair = (,) <$> gen <*> gen - msg = "in range " ++ show (toNano mi,toNano mx) ++ ": " ++ descr - gen = toEnum <$> choose (fromEnum mi, fromEnum mx) - - toNano :: Angle -> Nano - toNano = toDegree - -isIncreasing :: (Angle, Angle) -> Bool -isIncreasing (x , y) - | x == y = xA == yA - | x < y = xA < yA - | otherwise = xA > yA - where xA = lat x - yA = lat y - - -isDecreasing :: (Angle, Angle) -> Bool -isDecreasing (x, y) - | x == y = xA == yA - | x > y = xA < yA - | otherwise = xA > yA - where xA = lat x - yA = lat y - -shouldBeBounded :: (Arbitrary a, Ord a, Show a, Bounded a) => a -> Spec -shouldBeBounded a = prop msg $ \ x -> x >= mi && x <= mx - where msg = unwords [ "should lie between" - , show mi - , "and" - , show mx - ] - mi = minBound `asTypeOf` a - mx = maxBound `asTypeOf` a - - -spec :: Spec -spec = do - - describe "latitudes" $ do - - inRange (degree (-90), degree 90) "increases monotonically" isIncreasing - inRange (degree 90 , maxBound) "decreases monotonically" isDecreasing - inRange (minBound , degree (-90)) "decreases monotonically" isDecreasing - - shouldBeBounded (undefined :: Latitude) - - describe "longitudes" $ do - - prop "should have a period of 360 deg" $ - \ (x :: Longitude) -> x <> lon (degree 360) `shouldBe` x - - shouldBeBounded (undefined ::Longitude)