diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,5 @@
+# Changlog for persistent-spatial
+
+## 0.1.0.0
+
+Initial release, start of versioning.
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright © 2018-2019 Satsuma Labs
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,4 @@
+# persistent-spatial
+
+This package implements a type for storing and indexing geographic coordinates which can be used with any database which supports Word64. 
+See the haddocks for details.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/persistent-spatial.cabal b/persistent-spatial.cabal
new file mode 100644
--- /dev/null
+++ b/persistent-spatial.cabal
@@ -0,0 +1,53 @@
+name:                persistent-spatial
+version:             0.1.0.0
+synopsis:            Database agnostic, spatially indexed type for geographic points.
+description:         Defines type for storing geographic coordinates that can be spatially indexed by any database which supports Word64.
+                     This inxeding is implemented using a normal integer index on points repersented using a Morton Z-Order curve.
+                     Geographic regions are transformed into a covering set of tiles (contigious ranges) which can be used in a single query.
+homepage:            https://github.com/george-steel/persistent-spatial#readme
+license:             MIT
+license-file:        LICENSE
+author:              George Steel
+maintainer:          george.steel@gmail.com
+copyright:           2019 Satsuma Labs
+category:            Database, Geography
+build-type:          Simple
+cabal-version:       >=1.10
+extra-source-files:  README.md
+                   , CHANGELOG.md
+
+library
+  hs-source-dirs:      src
+  ghc-options:         -Wall -Wcompat -Wincomplete-record-updates -Wincomplete-uni-patterns -Wredundant-constraints
+  exposed-modules:     Data.Morton
+                       Data.LatLong
+  build-depends:       base >= 4.7 && < 5
+                     , integer-logarithms
+                     , aeson
+                     , lens
+                     , http-api-data
+                     , text
+                     , persistent
+  default-language:    Haskell2010
+
+
+test-suite tests
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      test
+  ghc-options:         -Wall -Wcompat -Wincomplete-record-updates -Wincomplete-uni-patterns -Wredundant-constraints -Wno-unused-imports
+  main-is:             Spec.hs
+  build-depends:       persistent-spatial
+                     , base >= 4.7 && < 5
+                     , aeson
+                     , http-api-data
+                     , text
+                     , persistent
+                     , hspec
+                     , QuickCheck
+  default-language:    Haskell2010
+
+
+
+source-repository head
+  type:     git
+  location: https://github.com/george-steel/persistent-spatial
diff --git a/src/Data/LatLong.hs b/src/Data/LatLong.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/LatLong.hs
@@ -0,0 +1,179 @@
+{-# 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
diff --git a/src/Data/Morton.hs b/src/Data/Morton.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Morton.hs
@@ -0,0 +1,224 @@
+{-# LANGUAGE BangPatterns, ScopedTypeVariables, GeneralizedNewtypeDeriving, DeriveFunctor, ViewPatterns, PatternSynonyms #-}
+
+{- |
+Module: Data.Morton
+Description: Morton reperesention of integer pairs
+Copyright: © 2018-2019 Satsuma labs
+
+Morton reperesentation of integer pairs (interleaved bits) used for creating spatial indexes.
+
+Bit interleaving code is originally from the documentation of Data.Sparse by Edward Kmett at
+<https://www.schoolofhaskell.com/user/edwardk/revisiting-matrix-multiplication/part-1>
+-}
+
+module Data.Morton (
+    Morton(Morton,MortonPair),
+    -- * Intervals
+    Interval(..), intersectInterval, intervalElem, intervalSize, intervalSizeMorton,
+    -- * Rectangles
+    MortonRect(MortonRect,MortonRectSides), mortonRectBounds, intersectMortonRect, mortonRectSize,
+    -- * Tiles
+    MortonTile(..), mortonTileBounds, mortonTileRect, enclosingMortonTile, splitMortonTile, trimMortonTile,
+    mortonTileCoverSized, mortonTileCover, mortonTileCoverTorus,
+) where
+
+import Data.Bits
+import Data.Word
+import Data.Monoid ((<>))
+import Data.Maybe
+import Data.Ord
+import Control.Monad
+import Numeric
+import Text.Read
+import Text.Read.Lex
+import Text.ParserCombinators.ReadPrec (readP_to_Prec)
+import Text.ParserCombinators.ReadP
+import Math.NumberTheory.Logarithms
+
+zeropad :: Int -> String -> String
+zeropad n s = replicate (n - length s) '0' ++ s
+
+-- | Type implementing a Morton Z-Order Curve.
+-- Stores two 'Word32' values with bits interleaved.
+-- This allows for spatial indexing by rectangular tiles which form contiguous intervals.
+newtype Morton = Morton Word64 deriving (Eq, Ord, Enum)
+-- Shows value in hex.
+instance (Show Morton) where
+    show (Morton m) = "Z" <> zeropad 16 (showHex m [])
+instance (Read Morton) where
+    readPrec = readP_to_Prec (const readMorton)
+
+readMorton :: ReadP Morton
+readMorton = char 'Z' >> fmap Morton readHexP
+
+
+-- interleaves the bits of two integers by performing AH,AL,BH,BL -> AH,BH,AL,BL then recursing on H and L portions in SIMD fashion
+interleaveM :: Word32 -> Word32 -> Morton
+interleaveM !x !y = Morton k5 where
+    k0 = unsafeShiftL (fromIntegral x) 32 .|. fromIntegral y
+    k1 = unsafeShiftL (k0 .&. 0x00000000FFFF0000) 16  .|. unsafeShiftR k0 16 .&. 0x00000000FFFF0000  .|. k0 .&. 0xFFFF00000000FFFF
+    k2 = unsafeShiftL (k1 .&. 0x0000FF000000FF00) 8   .|. unsafeShiftR k1 8  .&. 0x0000FF000000FF00  .|. k1 .&. 0xFF0000FFFF0000FF
+    k3 = unsafeShiftL (k2 .&. 0x00F000F000F000F0) 4   .|. unsafeShiftR k2 4  .&. 0x00F000F000F000F0  .|. k2 .&. 0xF00FF00FF00FF00F
+    k4 = unsafeShiftL (k3 .&. 0x0C0C0C0C0C0C0C0C) 2   .|. unsafeShiftR k3 2  .&. 0x0C0C0C0C0C0C0C0C  .|. k3 .&. 0xC3C3C3C3C3C3C3C3
+    k5 = unsafeShiftL (k4 .&. 0x2222222222222222) 1   .|. unsafeShiftR k4 1  .&. 0x2222222222222222  .|. k4 .&. 0x9999999999999999
+
+
+uninterleaveM :: Morton -> (Word32,Word32)
+uninterleaveM (Morton k0) = (fromIntegral (unsafeShiftR k5 32), fromIntegral (k5 .&. 0x00000000FFFFFFFF)) where
+    k5 = unsafeShiftL (k4 .&. 0x00000000FFFF0000) 16  .|. unsafeShiftR k4 16 .&. 0x00000000FFFF0000  .|. k4 .&. 0xFFFF00000000FFFF
+    k4 = unsafeShiftL (k3 .&. 0x0000FF000000FF00) 8   .|. unsafeShiftR k3 8  .&. 0x0000FF000000FF00  .|. k3 .&. 0xFF0000FFFF0000FF
+    k3 = unsafeShiftL (k2 .&. 0x00F000F000F000F0) 4   .|. unsafeShiftR k2 4  .&. 0x00F000F000F000F0  .|. k2 .&. 0xF00FF00FF00FF00F
+    k2 = unsafeShiftL (k1 .&. 0x0C0C0C0C0C0C0C0C) 2   .|. unsafeShiftR k1 2  .&. 0x0C0C0C0C0C0C0C0C  .|. k1 .&. 0xC3C3C3C3C3C3C3C3
+    k1 = unsafeShiftL (k0 .&. 0x2222222222222222) 1   .|. unsafeShiftR k0 1  .&. 0x2222222222222222  .|. k0 .&. 0x9999999999999999
+
+-- | Construct a Morton value from its two coordinates.
+pattern MortonPair :: Word32 -> Word32 -> Morton
+pattern MortonPair x y <- (uninterleaveM -> (x,y)) where
+    MortonPair x y = interleaveM x y
+{-# COMPLETE MortonPair #-}
+
+
+-- | Type for closed intervals. The second field should be greater than the first.
+data Interval a = Interval a a deriving (Eq, Show, Read, Functor)
+
+-- | Returns intersection of two intervals, or Nothing if they do not overlap.
+intersectInterval :: Ord a => Interval a -> Interval a -> Maybe (Interval a)
+intersectInterval (Interval a b) (Interval a' b') = do
+    guard (a <= b)
+    guard (a' <= b')
+    let x = max a a'
+        y = min b b'
+    guard (x <= y)
+    return $ Interval x y
+
+-- | Tests whether an element is contained within a given Interval.
+intervalElem :: (Ord a) => a -> Interval a -> Bool
+intervalElem x (Interval a b) = a <= x && x <= b
+
+-- | Returns the size of an integer interval.
+intervalSize :: (Integral a) => Interval a -> Integer
+intervalSize (Interval a b) = fromIntegral b - fromIntegral a + 1
+
+-- | Returns the size of a Morton interval.
+intervalSizeMorton :: Interval Morton -> Integer
+intervalSizeMorton (Interval (Morton a) (Morton b)) = fromIntegral b - fromIntegral a + 1
+
+
+-- | Type for retangles in Morton space reperesented by upper-left and lower-right corners
+data MortonRect = MortonRect {-# UNPACK #-} !Morton {-# UNPACK #-} !Morton deriving (Eq, Show, Read)
+
+-- | Returns x,y bounds of a rectangle
+mortonRectBounds :: MortonRect -> (Interval Word32, Interval Word32)
+mortonRectBounds (MortonRect (MortonPair x y) (MortonPair x' y')) = (Interval x x', Interval y y')
+
+-- | Construct/match rectangles by their sides.
+pattern MortonRectSides :: Interval Word32 -> Interval Word32 -> MortonRect
+pattern MortonRectSides xs ys <- (mortonRectBounds -> (xs,ys)) where
+    MortonRectSides (Interval x x') (Interval y y') = MortonRect (MortonPair x y) (MortonPair x' y')
+{-# COMPLETE MortonRectSides #-}
+
+-- | Rerurns intersection of two rectangles.
+intersectMortonRect :: MortonRect -> MortonRect -> Maybe MortonRect
+intersectMortonRect (MortonRect a1 b1) (MortonRect a2 b2) = mint where
+    selx (Morton m) = m .&. 0xAAAAAAAAAAAAAAAA
+    sely (Morton m) = m .&. 0x5555555555555555
+    -- interleaving bits with 0 preserves ordering so bit shifts are not necesary for comparisons
+    ax = max (selx a1) (selx a2)
+    ay = max (sely a1) (sely a2)
+    bx = min (selx b1) (selx b2)
+    by = min (sely b1) (sely b2)
+    mint = if (ax <= bx) && (ay <= by) then Just (MortonRect (Morton $ ax .|. ay) (Morton $ bx .|. by)) else Nothing
+
+-- | Returns the area of a rectangle
+mortonRectSize :: MortonRect -> Integer
+mortonRectSize (MortonRect (MortonPair ax ay) (MortonPair bx by)) =
+    (fromIntegral bx - fromIntegral ax + 1) * (fromIntegral by - fromIntegral ay + 1)
+
+
+-- | Type for a tile in Morton space, a special type of rectangle which is the set of all points sharing a common binary prefex.
+-- Reperesented as a point and mask length simillarly to a CIDR subnet.
+data MortonTile = MortonTile {-# UNPACK #-} !Morton {-# UNPACK #-} !Int
+instance (Show MortonTile) where
+    show t@(MortonTile _ n) = let Interval m _ = mortonTileBounds t in show m <> "/" <> show n
+instance (Read MortonTile) where
+    readPrec = readP_to_Prec . const $ do
+        m <- readMorton
+        char '/'
+        n <- readDecP
+        guard (n >= 0 && n <= 64)
+        return $ MortonTile m n
+
+-- | Values which reperesent the same tile compare equal even if the reperesentative points differ.
+instance (Eq MortonTile) where
+    a == b = mortonTileBounds a == mortonTileBounds b
+
+-- | A tile sorts immediately before its subtiles, i.e. x sorts before 0 and 1.
+instance (Ord MortonTile) where
+    compare = comparing (\x -> let Interval a b = mortonTileBounds x in (a, Down b))
+
+-- | Returns a tile as an 'Interval'.
+mortonTileBounds :: MortonTile -> Interval Morton
+mortonTileBounds (MortonTile (Morton x) n) = let
+    mask = shiftR 0xFFFFFFFFFFFFFFFF n
+    a = x .&. complement mask
+    b = x .|. mask
+    in Interval (Morton a) (Morton b)
+
+-- | Returns a tile as a 'MortonRect'.
+mortonTileRect :: MortonTile -> MortonRect
+mortonTileRect t = let
+    Interval a b = mortonTileBounds t
+    in MortonRect a b
+
+-- | Finds the smallest tile completely enclosing a rectangle.
+--  This can be arbitrarily large if the rectangle crosses a seam.
+enclosingMortonTile :: MortonRect -> MortonTile
+enclosingMortonTile (MortonRect (Morton a) (Morton b)) =
+    let n = countLeadingZeros $ a `xor` b in MortonTile (Morton a) n
+
+-- | Splits a 'MortonTile' in half. Does not split tiles containing a single value.
+splitMortonTile :: MortonTile -> [MortonTile]
+splitMortonTile t@(MortonTile _ 64) = [t]
+splitMortonTile (MortonTile (Morton m) n) = let
+    mask = shiftR 0xFFFFFFFFFFFFFFFF n
+    low = m .&. complement mask
+    high = m .|. mask
+    in [MortonTile (Morton low) (n+1), MortonTile (Morton high) (n+1)]
+
+-- | Trims a @MortonTile@ to its subtile overlapping a given rectangle.
+-- Returns 'Nothing' if the rectabgle and tile do not intersect.
+trimMortonTile :: MortonRect -> MortonTile -> Maybe MortonTile
+trimMortonTile rect t = let
+    Interval a b = mortonTileBounds t
+    trect = MortonRect a b
+    irect = intersectMortonRect rect trect
+    in fmap enclosingMortonTile irect
+
+-- | Covers a rectangle using tiles within a range of sizes (specified by their mask values).
+mortonTileCoverSized :: Int -> Maybe Int -> MortonRect -> [MortonTile]
+mortonTileCoverSized big small rect = let
+    expandTile tile@(MortonTile p m) = case small of
+        Just s | m > s -> MortonTile p s
+        _              -> tile
+    subdiv tile@(MortonTile _ m)
+        | m >= big = [expandTile tile]
+        | otherwise = (>>= subdiv) . mapMaybe (trimMortonTile rect) . splitMortonTile $ tile
+    in subdiv (enclosingMortonTile rect)
+
+-- | Covers a rectangle with tiles no larger then the area to be covered (no lower size limit).
+-- The total area coverd by these tiles bas a trivial upper bound of 8 tiles the rectangle's area plus the area of its enclosing square
+-- and the actual performance is usually significantly better (possibly always, although I have not proven so).
+mortonTileCover :: MortonRect -> [MortonTile]
+mortonTileCover rect = let
+    big = max 0 $ 64 - integerLog2 (mortonRectSize rect)
+    in mortonTileCoverSized big Nothing rect
+
+-- | Version of 'mortonTileCover' which allows the rectangle to wrap around the maximum x/y coordinates (as if the space were a torus).
+mortonTileCoverTorus :: Morton -> Morton -> [MortonTile]
+mortonTileCoverTorus (MortonPair ax ay) (MortonPair bx by) = let
+    rsize = (fromIntegral (bx - ax) + 1) * (fromIntegral (by - ay) + 1) -- integer overflows handle wraparound case
+    big = max 0 $ 64 - integerLog2 rsize
+    initrects = MortonRectSides
+        <$> (if bx >= ax then [Interval ax bx] else [Interval ax maxBound, Interval minBound bx])
+        <*> (if by >= ay then [Interval ay by] else [Interval ay maxBound, Interval minBound by])
+    in initrects >>= mortonTileCoverSized big Nothing
diff --git a/test/Spec.hs b/test/Spec.hs
new file mode 100644
--- /dev/null
+++ b/test/Spec.hs
@@ -0,0 +1,159 @@
+{-# OPTIONS_GHC -Wno-orphans #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Main where
+
+import Test.Hspec
+import Test.QuickCheck
+import Test.Hspec.QuickCheck
+import Data.Word
+import Data.List
+import Data.Maybe
+import Data.Aeson
+import Text.Read
+import Numeric
+import Web.HttpApiData
+import Database.Persist
+
+import Data.Morton
+import Data.LatLong
+
+instance (Arbitrary Morton) where
+    arbitrary = fmap Morton arbitrary
+instance (Arbitrary a, Ord a) => (Arbitrary (Interval a)) where
+    arbitrary = do
+        x <- arbitrary
+        y <- arbitrary
+        return $ Interval (min x y) (max x y)
+instance (Arbitrary MortonRect) where
+    arbitrary = MortonRectSides <$> arbitrary <*> arbitrary
+
+prop_morton_isom :: Word32 -> Word32 -> Bool
+prop_morton_isom x y = let
+    z = MortonPair x y
+    MortonPair a b = z
+    in (x == a) && (y == b)
+
+prop_morton_isom_rev :: Morton -> Bool
+prop_morton_isom_rev z = let
+    MortonPair x y = z
+    z' = MortonPair x y
+    in z == z'
+
+prop_morton_parse :: Morton -> Bool
+prop_morton_parse x = readMaybe (show x) == Just x
+
+prop_morton_intersect :: Interval Word32 -> Interval Word32 -> Interval Word32 -> Interval Word32 -> Property
+prop_morton_intersect a b c d = let
+    rx = MortonRectSides a b
+    ry = MortonRectSides c d
+    rz = intersectMortonRect rx ry
+    rz' = MortonRectSides <$> (intersectInterval a c) <*> (intersectInterval b d)
+    in classify (isJust rz) "overlapping" $ rz == rz'
+
+prop_morton_tile_parse :: Morton -> Property
+prop_morton_tile_parse m = forAll (choose (0,64)) $ \n -> let
+    t = MortonTile m n
+    mt = readMaybe (show t)
+    in mt == Just t
+
+prop_morton_tile_bounds :: MortonRect -> Bool
+prop_morton_tile_bounds rect = let
+    t@(MortonTile _ n) = enclosingMortonTile rect
+    s = intervalSizeMorton (mortonTileBounds t)
+    sr = mortonRectSize (mortonTileRect t)
+    s' = 2 ^ (64 - n)
+    in s == s' && sr == s'
+
+
+mortonRectAspect :: MortonRect -> Double
+mortonRectAspect (MortonRectSides ix iy) = let
+    dx = intervalSize ix
+    dy = intervalSize iy
+    in fromIntegral (max dx dy) / fromIntegral (min dx dy)
+
+mortonRectExpansion :: MortonRect -> Integer
+mortonRectExpansion rect = let
+    tiles = mortonTileCover rect
+    rsize = mortonRectSize rect
+    tsize = sum . fmap (intervalSizeMorton . mortonTileBounds) $ tiles
+    in (tsize + 1) `div` rsize
+
+prop_morton_split_8 :: MortonRect -> Property
+prop_morton_split_8 rect = let
+    ratio = mortonRectExpansion rect
+    in collect ratio $ mortonRectAspect rect < 8 ==> (ratio >= 1 && ratio <= 16)
+
+badSquare = MortonRect (Morton 0x0fffFfffFfffFfff) (Morton 0xc000000000000000)
+
+
+
+instance (Arbitrary LatLong) where
+    arbitrary = LatLong <$> choose (-90,90) <*> choose (-180,180)
+
+prop_latlong_json :: LatLong -> Bool
+prop_latlong_json p = let
+    mp = decode (encode p)
+    in mp == Just p
+
+prop_latlong_http :: LatLong -> Bool
+prop_latlong_http p = let
+    mp = parseUrlPiece (toUrlPiece p)
+    in mp == Right p
+
+prop_latlong_persist :: LatLong -> Bool
+prop_latlong_persist p = let
+    mp = fromPersistValue (toPersistValue p)
+    in mp == Right p
+
+prop_geo_triangle :: LatLong -> LatLong -> LatLong -> Bool
+prop_geo_triangle a b c = let
+    ttest x y z = geoDistance x y + geoDistance y z >= geoDistance x z
+    in ttest a b c && ttest b c a && ttest c a b
+
+pctError :: Double -> Double -> Double
+pctError ref samp = abs (samp - ref) / ref
+
+prop_square_corner_dist :: Property
+prop_square_corner_dist = let
+    gen = (,) <$> choose (10,200000) <*> (LatLong <$> choose (-70,70) <*> choose (-180,180))
+    test (r,p) = let
+        d = 2*r
+        (LatLong s w, LatLong n e) = geoSquare p r
+        laterr =  max (pctError d (geoDistance (LatLong n w) (LatLong s w)))
+                      (pctError d (geoDistance (LatLong n e) (LatLong s e)))
+        --laterrbin :: Int = max (-9) (ceiling (logBase 10 laterr))
+        longerr = max (pctError d (geoDistance (LatLong n w) (LatLong n e)))
+                      (pctError d (geoDistance (LatLong s w) (LatLong s e)))
+        longerrbin = minimum [x | x <- [100,10,5,2,1,0.5,0.1,0.01], x > longerr * 100]
+        in collect longerrbin . counterexample (show (laterr,longerr)) $ laterr < 1e-3 && longerr < 0.1
+    in forAll gen test
+
+prop_tiles_cover_points :: Property
+prop_tiles_cover_points = let
+    gen = (,) <$> choose (10,1000000) <*> (LatLong <$> choose (-70,70) <*> choose (-180,180))
+    test (r,p) = let
+        (sw@(LatLong s w), ne@(LatLong n e)) = geoSquare p r
+        tiles = latLongTileCover sw ne
+        pgen = LatLong <$> choose (s,n) <*> choose (w,e)
+        in forAll pgen (`tileSetElem` tiles)
+    in forAll gen test
+
+main :: IO ()
+main = hspec . modifyMaxSuccess (const 10000) $ do
+    describe "Data.Morton" $ do
+        prop "interleaving is reversible" prop_morton_isom
+        prop "deinterleaving is reversible" prop_morton_isom_rev
+        prop "point Read instance works" prop_morton_parse
+        prop "rectangle intersection mathes interval intersection" prop_morton_intersect
+        prop "tile Read instance works" prop_morton_tile_parse
+        prop "tile size (both definitions) matches mask value" prop_morton_tile_bounds
+        it   "pathological square expands at most 6-fold" $ mortonRectExpansion badSquare `shouldSatisfy` (< 6)
+        modifyMaxSuccess (const 100000) $ prop "rectangle expansion tile cover (eccentricity limit 8, floored expansion collected)" prop_morton_split_8
+    describe "Data.LatLong" $ do
+        prop "Aeson instances work" prop_latlong_json
+        prop "HttpApiData instances work" prop_latlong_http
+        prop "PersistField instance works" prop_latlong_persist
+        prop "geoDistance obeys triangle inequality" prop_geo_triangle
+        prop "geoSquare corners within tolerance (% error collected), 10m < r < 200km, 70S < lat < 70N" prop_square_corner_dist
+        prop "tile covers contain all points in square" prop_tiles_cover_points