diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,3 @@ +## 0.1.0 + +- Initial release. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright Colin Woodbury (c) 2017 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Author name here nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,149 @@ +# mapalgebra + +*Efficient, polymorphic Map Algebra for Haskell.* + +This library is an implementation of *Map Algebra* as described in the +book *GIS and Cartographic Modeling* by Dana Tomlin. The fundamental +primitive is the `Raster`, a rectangular grid of data that usually describes +some area on the earth. + +`mapalgebra` is built on top of [massiv](https://github.com/lehins/massiv), +a powerful Parallel Array library by Alexey Kuleshevich. + +## Usage + +Always compile with `-threaded`, `-O2`, and `-with-rtsopts=-N` for best +performance. + +### The `Raster` Type + +This library provides `Raster`s which are lazy, polymorphic, and typesafe. They +can hold any kind of data, and are aware of their projection and dimensions +at the type level. This means that imagery of different size or projection +are considered completely different types, which prevents an entire class +of potential bugs. + +`Raster`s have types signatures like this: + +```haskell +-- | A Raster of Ints backed by efficient byte-packed arrays, encoded +-- via the `Storable` typeclass. `P` (Prim), `U` (Unbox) and `B` (boxed) are also available. +-- +-- This is either a freshly read image, or the result of evaluating a "delayed" +-- (`D` or `DW`) Raster. +Raster S LatLng 256 256 Int + +-- | A "delayed" Raster of bytes. Likely the result of some Local Operation. +-- Waiting to be evaluated by the `strict` function. +Raster D WebMercator 512 512 Word8 + +-- | A "windowed" Raster of an ADT, the result of some Focal Operation. +-- Waiting to be evaluated by the `strict` function. +Raster DW p 1024 1024 (Maybe Double) +``` + +### Reading Imagery + +`mapalgebra` can currently read any image file of any value type, so long as +it is grayscale (singleband) or RGBA. True multiband rasters (like from LandSat) +are not yet supported. + +To read a Raster: + +```haskell +-- | You must know the image dimensions ahead of time. If you don't care +-- about the projection, then `p` can be left generic. +getRaster :: IO (Raster S p 512 512 Word8) +getRaster = do + erast <- fromGray "path/to/image.tif" + case erast of + Left err -> ... -- deal with the error. + Right r -> pure r +``` + +### Colouring and Viewing Imagery + +To quickly view a Raster you're working on, use the `display` function: + +```haskell +-- | Simplified type signature. +display :: Raster D p r c a -> IO () +``` + +This will automatically colour gray, evaluate, and display your Raster +using your OS's default image viewer. + +To colour a Raster gray yourself, use `grayscale`: + +```haskell +grayscale :: Functor (Raster u p r c) => Raster u p r c a -> Raster u p r c (Pixel Y a) +``` + +True colouring is done with the `classify` function and colour ramps inspired by +Gretchen N. Peterson's book *Cartographer's Toolkit*. + +```haskell +-- | Both `Raster D` and `Raster DW` are Functors, so this function works on +-- either of them. `Raster S`, etc., do not form Functors by design. +classify :: (Ord a, Functor f) => b -> Map a b -> f a -> f b + +-- | An invisible pixel (alpha channel set to 0) to be passed +-- to `classify` as a default. +invisible :: Pixel RGBA Word8 + +-- | Given a list of "breaks", forms a colour ramp to be passed +-- to `classify`. +greenRed :: Ord k => [k] -> Map k (Pixel RGBA Word8) +``` + +### Local Operations + +All Local Operations defined in *GIS and Cartographic Modeling* are available. +For the usual math ops, `Raster D` has a `Num` instance: + +```haskell +rast :: Raster D p 512 512 Int + +squared :: Raster D p 512 512 Int +squared = rast * rast -- Element-wise multiplication. +``` + +### Focal Operations + +Except for *Focal Ranking* and *Focal Insularity*, all Focal Operations of immedatiate +neighbourhoods are provided: + +```haskell +rast :: Raster S p 512 512 Double + +-- | `Raster DW` forms a Functor, so we can do simple unary transformations +-- (like colouring!) to it after Focal Ops. +averagedPlusAbit :: Raster S p 512 512 Double +averagedPlusAbit = strict S . fmap (+1) $ fmean rast +``` + +### Typesafe NoData Handling + +If it's known that your images have large areas of NoData, consider that `Maybe` +has a `Monoid` instance: + +```haskell +import Data.Monoid (Sum(..)) + +nodatafsum :: Raster S p r c Word8 -> Raster DW p r c 512 Word8 +nodatafsum = fmap (maybe 0 getSum) . fmonoid . strict B . fmap check . lazy + where check 0 = Nothing + check n = Just $ Sum n +``` + +In theory, one could construct special `newtype` wrappers with `Monoid` instances +that handle any Focal scenario imaginable. + +## Future Work + +- Projection handling at IO time +- Histograms for colour ramp generation +- Reprojections +- Extended neighbourhoods for Focal Ops +- Upsampling and Downsampling +- Improved NoData handling 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/bench/Bench.hs b/bench/Bench.hs new file mode 100644 --- /dev/null +++ b/bench/Bench.hs @@ -0,0 +1,250 @@ +{-# LANGUAGE DataKinds, TypeApplications #-} + +module Main ( main ) where + +import Criterion.Main +import Data.List.NonEmpty (NonEmpty(..)) +import Data.Massiv.Array as A hiding (zipWith) +import Data.Monoid ((<>), Sum(..)) +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as U +import GHC.TypeLits +import Geography.MapAlgebra +import Graphics.ColorSpace +import qualified Numeric.LinearAlgebra as LA +import Prelude hiding (zipWith) + +--- + +main :: IO () +main = do + img <- fileY "data/gray512.tif" + imgF <- fileY "data/gray512.tif" + rgba <- fileRGB "data/512x512.tif" + a512 <- fileRGB "data/512x512.tif" + -- a1024 <- fileRGB "data/1024x1024.tif" + -- a2048 <- fileRGB "data/2048x2048.tif" + -- a4096 <- fileRGB "data/4096x4096.tif" + -- let aHuge :: Raster S u 46500 46500 Word8 + -- aHuge = fromFunction S Par (\(r :. c) -> fromIntegral $ r * c) + defaultMain + [ + creation + , io + , colouring img + , massivOps img + , localOps rgba a512 img + , hmatrix + , conversions img + , focalOps img imgF + , compositeOps a512 -- a1024 a2048 a4096 + ] + +creation :: Benchmark +creation = bgroup "Raster Creation" + [ bench "constant 256x256" $ nf (_array . constantB) 5 + , bench "constant 512x512" $ nf (_array . constantB') 5 + , bench "fromFunction 256x256" $ nf (_array . functionB) (\(r :. c) -> r * c) + , bench "fromFunction 512x512" $ nf (_array . functionB') (\(r :. c) -> r * c) + , bench "fromVector Unboxed Int - 256x256" $ nf (_array . vectorB) uv + , bench "fromVector Boxed Int - 256x256" $ nf (_array . vectorB') bv ] + where uv = U.fromList ([1..65536] :: [Int]) + bv = V.fromList ([1..65536] :: [Int]) + +io :: Benchmark +io = bgroup "IO" + [ bench "fromRGBA 512x512" $ nfIO (_array . _red <$> fileRGB @Word8 @512 @512 "data/512x512.tif") + , bench "fromRGBA (Word8 -> Double) 512x512" $ nfIO (_array . _red <$> fileRGB @Double @512 @512 "data/512x512.tif") + , bench "fromGray Multiband 512x512" $ nfIO (_array <$> fileY @Word8 @512 @512 "data/512x512.tif") + , bench "fromGray Singleband 512x512" $ nfIO (_array <$> gray "data/gray512.tif") ] + +colouring :: Raster S p 512 512 Word8 -> Benchmark +colouring img = bgroup "Colouring" + [ bench "classify 512x512" $ nf (_array . strict S . classify invisible cr . lazy) img + , bench "grayscale 512x512" $ nf (_array . strict S . grayscale . lazy) img ] + where cr = greenRed [25, 50, 75, 100, 125, 150, 175, 200, 225, 255] + +massivOps :: Raster S p 512 512 Word8 -> Benchmark +massivOps img = bgroup "Massiv Operations" + [ bench "strict S . lazy" $ nf (_array . strict S . lazy) img + , bench "strict U . lazy" $ nf (_array . strict U . lazy) img + , bench "strict P . lazy" $ nf (_array . strict P . lazy) img ] + +localOps :: RGBARaster p 512 512 Word8 -> RGBARaster p 512 512 Double -> Raster S p 512 512 Word8 -> Benchmark +localOps (RGBARaster r g b _) (RGBARaster rF gF bF _) img = bgroup "Local Operations" + [ bench "fmap (+ 17) . lazy" $ nf (_array . strict S . fmap (+ 17) . lazy) img + , bench "zipWith (+)" $ nf (_array . strict S . zipWith (+) r) g + , bench "zipWith (/)" $ nf (_array . strict S . zipWith (/) doubles) doubles + , bench "(+)" $ nf (_array . strict S . (+ lazy r)) (lazy g) + , bench "(/)" $ nf (_array . strict S . (/ lazy doubles)) (lazy doubles) + , bench "lmax" $ nf (_array . strict S . lmax img) img + , bench "lmin" $ nf (_array . strict S . lmin img) img + , bench "lmean (Word8)" $ nf (_array . strict S . lmean @Word8 @Double) rs + , bench "lmean (Double)" $ nf (_array . strict S . lmean @Double @Double) rsF + , bench "lvariety" $ nf (_array . strict S . lvariety) rs + , bench "lmajority" $ nf (_array . strict S . lmajority) rs + , bench "lminority" $ nf (_array . strict S . lminority) rs + , bench "lvariance (Word8)" $ nf (fmap (_array . strict S) . lvariance) rs + , bench "lvariance (Double)" $ nf (fmap (_array . strict S) . lvariance) rsF ] + where rs = lazy r :| [lazy g, lazy b] + rsF = lazy rF :| [lazy gF, lazy bF] + +hmatrix :: Benchmark +hmatrix = bgroup "HMatrix" + [ bench "linearSolveLS" $ nf (LA.linearSolveLS zing) (LA.col [8,8,8,8,8,8,8,8,8]) + , bench "manual - MxM" $ nf (leftPseudo <>) (LA.col [8,8,8,8,8,8,8,8,8]) + , bench "manual - MxV" $ nf (leftPseudo LA.#>) (LA.vector [8,8,8,8,8,8,8,8,8]) ] + +conversions :: Raster S p 512 512 Word8 -> Benchmark +conversions img = bgroup "Numeric Conversion" + [ bench "Double -> Double via id" $ nf id tau + , bench "Double -> Double via realToFrac" $ nf (realToFrac @Double @Double) tau + , bench "Word -> Double via realToFrac" $ nf (realToFrac @Word8 @Double) 5 + , bench "realToFrac on Raster" $ nf (_array . strict S . fmap (realToFrac @Word8 @Double) . lazy) img ] + +focalOps :: Raster S p 512 512 Word8 -> Raster S p 512 512 Double -> Benchmark +focalOps img imgF = bgroup "Focal Operations" + [ bgroup "fsum" + [ bench "512" $ nf (_array . strict S . fsum) img + -- , bench "46500" $ nf (_array . strict S . fsum) huge + -- , bench "maybing" $ nf (_array . strict B . fmap Just . lazy) img + , bench "nodata" $ nf (_array . strict S . nodatafsum) img + ] + , bgroup "fmean" + [ bench "Word8" $ nf (_array . strict S . fmean @Word8 @Double) img + , bench "Double" $ nf (_array . strict S . fmean @Double @Double) imgF + ] + , bench "fmax" $ nf (_array . strict S . fmax) img + , bench "fmin" $ nf (_array . strict S . fmin) img + , bench "fmajority" $ nf (_array . strict S . fmajority) img + , bench "fminority" $ nf (_array . strict S . fminority) img + , bench "fvariety" $ nf (_array . strict S . fvariety) img + , bench "fpercentage" $ nf (_array . strict S . fpercentage) img + , bench "fpercentile" $ nf (_array . strict S . fpercentile) img + , bench "flinkage" $ nf (_array . strict S . flinkage) img + , bench "flength" $ nf (_array . strict S . flength . flinkage) img + , bench "fpartition" $ nf (_array . strict B . fpartition) img + , bench "fshape" $ nf (_array . strict B . fshape) img + , bench "ffrontage" $ nf (_array . strict S . ffrontage . fshape) img + , bench "farea" $ nf (_array . strict S . farea . fshape) img + , bgroup "fvolume" + [ bench "Word8 -> Double" $ nf (_array . strict S . fvolume . strict S . fmap (realToFrac @Word8 @Double) . lazy) img + , bench "Double" $ nf (_array . strict S . fvolume) imgF + ] + , bgroup "fgradient" + [ bench "Word8" $ nf (_array . strict S . fgradient . strict S . fmap wtod . lazy) img + , bench "Double" $ nf (_array . strict S . fgradient) imgF + ] + , bgroup "faspect" + [ bench "Unsafe (Word8)" $ nf (_array . strict S . faspect' . strict S . fmap wtod . lazy) img + , bench "Unsafe (Double)" $ nf (_array .strict S . faspect') imgF + , bench "Safe (Word8)" $ nf (_array . strict B . faspect . strict S . fmap wtod . lazy) img + , bench "Safe (Double)" $ nf (_array . strict B . faspect) imgF + ] + , bgroup "fdownstream" + [ bench "Word8" $ nf (_array . strict S . fdownstream . strict S . fmap wtod . lazy) img + , bench "Double" $ nf (_array. strict S . fdownstream) imgF + ] + , bgroup "fupstream" + [ bench "Word8" $ nf (_array . strict S . fupstream . strict S . fdownstream . strict S . fmap wtod . lazy) img + , bench "Double" $ nf (_array . strict S . fupstream . strict S . fdownstream) imgF + ] + ] + +wtod :: Word8 -> Double +wtod = realToFrac + +compositeOps :: RGBARaster p 512 512 Double + -- -> RGBARaster p 1024 1024 Double + -- -> RGBARaster p 2048 2048 Double + -- -> RGBARaster p 4096 4096 Double + -> Benchmark +compositeOps i@(RGBARaster r g _ _) = + bgroup "Composite Operations" + [ bench "NDVI" $ nf (_array . strict S . ndvi (lazy g)) (lazy r) + , bench "EVI (512)" $ nf (_array . strict S . evi) i + -- , bench "EVI (1024)" $ nf (_array . strict S . evi) a1024 + -- , bench "EVI (2048)" $ nf (_array . strict S . evi) a2048 + -- , bench "EVI (4096)" $ nf (_array . strict S . evi) a4096 + , bench "EVI + Colour" $ nf (_array . strict S . classify invisible cr . evi) i + , bench "EVI + Colour + PNG (D)" $ nf (png . classify invisible cr . evi) i + , bench "EVI + Colour + PNG (S)" $ nf (png . strict S . classify invisible cr . evi) i ] + where cr = greenRed $ fmap (10 ^) [1..10] + +fromRight :: Either a b -> b +fromRight (Right b) = b +fromRight _ = error "Was Left" + +constantB :: Int -> Raster S p 256 256 Int +constantB = constant S Par + +constantB' :: Int -> Raster S p 512 512 Int +constantB' = constant S Par + +functionB :: (Ix2 -> Int) -> Raster S p 256 256 Int +functionB = fromFunction S Par + +functionB' :: (Ix2 -> Int) -> Raster S p 512 512 Int +functionB' = fromFunction S Par + +vectorB :: U.Vector Int -> Raster U p 256 256 Int +vectorB = fromRight . fromVector Par + +vectorB' :: V.Vector Int -> Raster B p 256 256 Int +vectorB' = fromRight . fromVector Par + +doubles :: Raster U p 512 512 Double +doubles = fromRight . fromVector Par $ U.fromList ([1..262144] :: [Double]) + +zing :: LA.Matrix Double +zing = LA.matrix 3 [ -0.5, -0.5, 1 + , -0.5, 0, 1 + , -0.5, 0.5, 1 + , 0, -0.5, 1 + , 0, 0, 1 + , 0, 0.5, 1 + , 0.5, -0.5, 1 + , 0.5, 0, 1 + , 0.5, 0.5, 1 ] + +gray :: FilePath -> IO (Raster S p 512 512 Word8) +gray fp = do + i <- fromGray fp + case i of + Left err -> error err + Right img -> pure img + +fileY :: (Elevator a, KnownNat r, KnownNat c) => FilePath -> IO (Raster S p r c a) +fileY fp = do + i <- fromGray fp + case i of + Left err -> error err + Right img -> pure img +{-# INLINE fileY #-} + +fileRGB :: (Elevator a, KnownNat r, KnownNat c) => FilePath -> IO (RGBARaster p r c a) +fileRGB fp = do + i <- fromRGBA fp + case i of + Left err -> error err + Right img -> pure img +{-# INLINE fileRGB #-} + +-- | See: https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index +ndvi :: (KnownNat r, KnownNat c) => Raster D p r c Double -> Raster D p r c Double -> Raster D p r c Double +ndvi nir red = (nir - red) / (nir + red) +{-# INLINE ndvi #-} + +-- | See: https://en.wikipedia.org/wiki/Enhanced_vegetation_index +evi :: (Fractional a, Storable a, KnownNat r, KnownNat c) => RGBARaster p r c a -> Raster D p r c a +evi (RGBARaster r g b _) = 2.5 * (numer / denom) + where nir = lazy g -- fudging it. + numer = nir - lazy r + denom = nir + (6 * lazy r) - (7.5 * lazy b) + 1 +{-# INLINE evi #-} + +nodatafsum :: Raster S p 512 512 Word8 -> Raster DW p 512 512 Word8 +nodatafsum = fmap (maybe 0 getSum) . fmonoid . strict B . fmap check . lazy + where check 0 = Nothing + check n = Just $ Sum n +{-# INLINE nodatafsum #-} diff --git a/data/1024x1024.tif b/data/1024x1024.tif new file mode 100644 # file too large to diff: data/1024x1024.tif diff --git a/data/256x256.tif b/data/256x256.tif new file mode 100644 Binary files /dev/null and b/data/256x256.tif differ diff --git a/data/512x512.tif b/data/512x512.tif new file mode 100644 Binary files /dev/null and b/data/512x512.tif differ diff --git a/data/gray512.tif b/data/gray512.tif new file mode 100644 Binary files /dev/null and b/data/gray512.tif differ diff --git a/lib/Geography/MapAlgebra.hs b/lib/Geography/MapAlgebra.hs new file mode 100644 --- /dev/null +++ b/lib/Geography/MapAlgebra.hs @@ -0,0 +1,1247 @@ +{-# LANGUAGE Rank2Types, DataKinds, KindSignatures, ScopedTypeVariables #-} +{-# LANGUAGE FlexibleInstances, FlexibleContexts #-} +{-# LANGUAGE ApplicativeDo, BangPatterns, UnboxedTuples, TypeInType #-} +{-# LANGUAGE DerivingStrategies, GeneralizedNewtypeDeriving, DeriveAnyClass #-} + +-- | +-- Module : Geography.MapAlgebra +-- Copyright : (c) Colin Woodbury, 2018 +-- License : BSD3 +-- Maintainer: Colin Woodbury +-- +-- This library is an implementation of /Map Algebra/ as described in the +-- book /GIS and Cartographic Modeling/ (GaCM) by Dana Tomlin. The fundamental +-- primitive is the `Raster`, a rectangular grid of data that usually describes +-- some area on the earth. A `Raster` need not contain numerical data, however, +-- and need not just represent satellite imagery. It is essentially a matrix, +-- which of course forms a `Functor`, and thus is available for all the +-- operations we would expect to run on any Functor. /GIS and Cartographic Modeling/ +-- doesn't lean on this fact, and so describes many seemingly custom +-- operations which to Haskell are just applications of `fmap` or `zipWith` +-- with pure functions. +-- +-- Here are the main classes of operations ascribed to /Map Algebra/ and their +-- corresponding approach in Haskell: +-- +-- * Single-Raster /Local Operations/ -> `fmap` with pure functions +-- * Multi-Raster /Local Operations/ -> `foldl` with `zipWith` and pure functions +-- * /Focal Operations/ -> 'massiv'-based smart `Stencil` operations +-- * /Zonal Operations/ -> Not yet implemented +-- +-- Whether it is meaningful to perform operations between two given +-- `Raster`s (i.e. whether the Rasters properly overlap on the earth) is not +-- handled in this library and is left to the application. +-- +-- The "colour ramp" generation functions (like `greenRed`) gratefully borrow colour +-- sets from Gretchen N. Peterson's book /Cartographer's Toolkit/. +-- +-- === A Word on Massiv: Fused, Parallel Arrays +-- Thanks to the underlying `Array` library [massiv](https://hackage.haskell.org/package/massiv), +-- most operations over and between Rasters are /fused/, meaning that no extra memory +-- is allocated in between each step of a composite operation. +-- +-- Take the [Enhanced Vegetation Index](https://en.wikipedia.org/wiki/Enhanced_vegetation_index) +-- calculation: +-- +-- @ +-- evi :: Raster D p r c Double -> Raster D p r c Double -> Raster D p r c Double -> Raster D p r c Double +-- evi nir red blue = 2.5 * (numer / denom) +-- where numer = nir - red +-- denom = nir + (6 * red) - (7.5 * blue) + 1 +-- @ +-- +-- 8 binary operators are used here, but none allocate new memory. It's only when +-- some `lazy` Raster is made `strict` that calculations occur and memory is allocated. +-- +-- Provided your machine has more than 1 CPU, Rasters read by functions like +-- `fromRGBA` will automatically be in `Par`allel mode. This means that +-- forcing calculations with `strict` will cause evaluation to be done +-- with every CPU your machine has. The effect of this is quite potent for Focal +-- Operations, which yield special, cache-friendly windowed (`DW`) Rasters. +-- +-- Familiarity with Massiv will help in using this library. A guide +-- [can be found here](https://github.com/lehins/massiv). +-- +-- === Compilation Options for Best Performance +-- When using this library, always compile your project with @-threaded@ and @-with-rtsopts=-N@. +-- These will ensure your executables will automatically use all the available CPU cores. +-- +-- As always, @-O2@ is your friend. The @{-\# INLINE ... \#-}@ pragma is also very likely +-- to improve the performance of code that uses functions from this library. Make sure +-- to benchmark proactively. +-- +-- For particularly mathy operations like `fmean`, compiling with @-fllvm@ grants +-- about a 2x speedup. + +module Geography.MapAlgebra + ( + -- * Types + -- ** Rasters + Raster(..) + , lazy, strict + , RGBARaster(..) + -- *** Creation + , constant, fromFunction, fromVector, fromRGBA, fromGray + -- *** Colouring + -- | The `M.Map`s here can be used with `classify` to + -- transform /ranges/ of values into certain colours in \(\mathcal{O}(n\log(n))\). + -- Each Map-generating function (like `greenRed`) creates a "colour ramp" of 10 colours. So, it expects + -- to be given a list of 10 "break points" which become the Map's keys. Any less than 10 will result + -- in the later colours not being used. Any more than 10 will be ignored. The list of break points is + -- assumed to be sorted. + -- `invisible` can be used as the default value to `classify`, to make invisible any value that falls outside + -- the range of the Maps. + -- + -- If you aren't interested in colour but still want to render your `Raster`, + -- consider `grayscale`. Coloured `Raster`s can be unwrapped with `_array` and then + -- output with functions like `writeImage`. + , grayscale + , invisible + , greenRed, spectrum, blueGreen, purpleYellow, brownBlue + , grayBrown, greenPurple, brownYellow, purpleGreen, purpleRed + -- *** Output and Display + -- | For coloured output, first use `classify` over your `Raster` to produce a + -- @Raster u p r c (Pixel RGBA Word8)@. Then unwrap it with `_array` and output + -- with something like `writeImage`. + -- + -- For quick debugging, you can visualize a `Raster` with `display`. + , writeImage, writeImageAuto + , png, display + -- ** Projections + , Projection(..) + , reproject + , Sphere, LatLng, WebMercator + , Point(..) + -- * Map Algebra + -- ** Local Operations + -- | Operations between `Raster`s. All operations are element-wise: + -- + -- @ + -- 1 1 + 2 2 == 3 3 + -- 1 1 2 2 3 3 + -- + -- 2 2 * 3 3 == 6 6 + -- 2 2 3 3 6 6 + -- @ + -- + -- If an operation you need isn't available here, use our `zipWith`: + -- + -- @ + -- zipWith :: (a -> b -> d) -> Raster u p r c a -> Raster u p r c b -> Raster D p r c d + -- + -- -- Your operation. + -- foo :: Int -> Int -> Int + -- + -- bar :: Raster u p r c Int -> Raster u p r c Int -> Raster D p r c Int + -- bar a b = zipWith foo a b + -- @ + , zipWith + -- *** Unary + -- | If you want to do simple unary @Raster -> Raster@ operations (called + -- /LocalCalculation/ in GaCM), `Raster` is a `Functor` so you can use + -- `fmap` as normal: + -- + -- @ + -- myRaster :: Raster D p r c Int + -- + -- -- Add 17 to every value in the Raster. + -- fmap (+ 17) myRaster + -- @ + , classify + -- *** Binary + -- | You can safely use these with the `foldl` family on any `Foldable` of + -- `Raster`s. You would likely want @foldl1'@ which is provided by both List + -- and Vector. + -- + -- Keep in mind that `Raster` `D` has a `Num` instance, so you can use + -- all the normal math operators with them as well. + , lmax, lmin + -- *** Other + -- | There is no binary form of these functions that exists without + -- producing numerical error, so you can't use the `foldl` family with these. + -- Consider the average operation, where the following is /not/ true: + -- \[ + -- \forall abc \in \mathbb{R}. \frac{\frac{a + b}{2} + c}{2} = \frac{a + b + c}{3} + -- \] + , lmean, lvariety, lmajority, lminority, lvariance + -- ** Focal Operations + -- | Operations on one `Raster`, given some polygonal neighbourhood. + -- Your `Raster` must be of a `Manifest` type (i.e. backed by real memory) before + -- you attempt any focal operations. Without this constraint, wayward users + -- run the risk of setting up operations that would perform terribly. + -- Use `strict` to easily convert a lazy `Raster` to a memory-backed one. + -- + -- @ + -- myRaster :: Raster D p r c Float + -- + -- averaged :: Raster DW p r c Float + -- averaged = fmean $ strict P myRaster + -- @ + , fsum, fproduct, fmonoid, fmean + , fmax, fmin + , fmajority, fminority, fvariety + , fpercentage, fpercentile + -- *** Lineal + -- | Focal operations that assume that groups of data points represent line-like objects + -- in a `Raster`. GaCM calls these /lineal characteristics/ and describes them fully + -- on page 18 and 19. + , Line(..) + , flinkage, flength + -- *** Areal + -- | Focal operations that assume that groups of data points represent 2D areas + -- in a `Raster`. GaCM calls these /areal characteristics/ and describes them fully + -- on page 20 and 21. + , Corners(..), Surround(..) + , fpartition, fshape, ffrontage, farea + -- *** Surficial + -- | Focal operations that work over elevation `Raster`s. GaCM calls elevation + -- features /surficial characteristics/ and describes them fully on page 21 + -- and 22. + -- + -- Some of these operations require finding a "best-fit plane" that + -- approximates the surficial shape of each pixel. Each pixel has 9 "facet points" + -- calculated for it based on its surrounding pixels. We then use these facets to determine + -- a plane which adheres to this equation: + -- + -- \[ + -- ax + by + c = z + -- \] + -- This is a linear equation that we can solve for in the form \(Ax = B\). + -- For facet points \((x_i, y_i, z_i)\), we have: + -- + -- \[ + -- \begin{bmatrix} + -- x_0 & y_0 & 1 \\ + -- x_1 & y_1 & 1 \\ + -- \vdots & \vdots & \vdots \\ + -- x_n & y_n & 1 + -- \end{bmatrix} \begin{bmatrix} + -- a\\ + -- b\\ + -- c + -- \end{bmatrix} = \begin{bmatrix} + -- z_0\\ + -- z_1\\ + -- \vdots\\ + -- z_n + -- \end{bmatrix} + -- \] + -- + -- Since this system of equations is "over determined", we rework the above to + -- find the coefficients of the best-fitting plane via: + -- \[ + -- \begin{bmatrix} + -- a\\ + -- b\\ + -- c + -- \end{bmatrix} = \boxed{(A^{T}A)^{-1}A^{T}}B + -- \] + -- The boxed section is called the "left pseudo inverse" and is available as `leftPseudo`. + -- The actual values of \(A\) don't matter for our purposes, hence \(A\) can be fixed to + -- avoid redundant calculations. + , Drain(..), Direction(..) + , direction, directions, drainage + , fvolume, fgradient, faspect, faspect', fdownstream, fupstream + -- * Utilities + , leftPseudo, tau + ) where + +import Control.Concurrent (getNumCapabilities) +import Control.DeepSeq (NFData(..), deepseq) +import Data.Bits (testBit) +import Data.Bool (bool) +import qualified Data.ByteString.Lazy as BL +import Data.Default (Default, def) +import Data.Foldable +import qualified Data.List as L +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.List.NonEmpty as NE +import qualified Data.Map.Strict as M +import qualified Data.Massiv.Array as A +import Data.Massiv.Array hiding (zipWith) +import Data.Massiv.Array.IO +import qualified Data.Massiv.Array.Manifest.Vector as A +import Data.Massiv.Array.Unsafe as A +import Data.Proxy (Proxy(..)) +import Data.Semigroup +import qualified Data.Set as S +import Data.Typeable (Typeable) +import qualified Data.Vector.Generic as GV +import qualified Data.Vector.Storable as VS +import Data.Word +import GHC.TypeLits +import Graphics.ColorSpace (Elevator, RGBA, Y, Pixel(..), ColorSpace) +import qualified Numeric.LinearAlgebra as LA +import qualified Prelude as P +import Prelude hiding (zipWith) +import Text.Printf (printf) + +-- + +-- | A location on the Earth in some `Projection`. +data Point p = Point { x :: !Double, y :: !Double } deriving (Eq, Show) + +-- | The Earth is not a sphere. Various schemes have been invented +-- throughout history that provide `Point` coordinates for locations on the +-- earth, although all are approximations and come with trade-offs. We call +-- these "Projections", since they are a mapping of `Sphere` coordinates to +-- some other approximation. The Projection used most commonly for mapping on +-- the internet is called `WebMercator`. +-- +-- A Projection is also known as a Coordinate Reference System (CRS). +-- +-- Use `reproject` to convert `Point`s between various Projections. +-- +-- __Note:__ Full support for Projections is still pending. +class Projection p where + -- | Convert a `Point` in this Projection to one of radians on a perfect `Sphere`. + toSphere :: Point p -> Point Sphere + + -- | Convert a `Point` of radians on a perfect sphere to that of a specific Projection. + fromSphere :: Point Sphere -> Point p + +-- | Reproject a `Point` from one `Projection` to another. +reproject :: (Projection p, Projection r) => Point p -> Point r +reproject = fromSphere . toSphere +{-# INLINE reproject #-} + +-- | A perfect geometric sphere. The earth isn't actually shaped this way, +-- but it's a convenient middle-ground for converting between various +-- `Projection`s. +data Sphere + +instance Projection Sphere where + toSphere = id + fromSphere = id + +-- | Latitude (north-south position) and Longitude (east-west position). +data LatLng + +-- instance Projection LatLng where +-- toSphere = undefined +-- fromSphere = undefined + +-- | The most commonly used `Projection` for mapping in internet applications. +data WebMercator + +-- instance Projection WebMercator where +-- toSphere = undefined +-- fromSphere = undefined + +-- | A rectangular grid of data representing some area on the earth. +-- +-- * @u@: What is the /underlying representation/ of this Raster? (see 'massiv') +-- * @p@: What `Projection` is this Raster in? +-- * @r@: How many rows does this Raster have? +-- * @c@: How many columns does this Raster have? +-- * @a@: What data type is held in this Raster? +-- +-- By having explicit p, r, and c, we make impossible any operation between +-- two Rasters of differing size or projection. Conceptually, we consider +-- Rasters of different size and projection to be /entirely different types/. +-- Example: +-- +-- @ +-- -- | A lazy 256x256 Raster with the value 5 at every index. Uses DataKinds +-- -- and "type literals" to achieve the same-size guarantee. +-- myRaster :: Raster D WebMercator 256 256 Int +-- myRaster = constant D Par 5 +-- +-- >>> length myRaster +-- 65536 +-- @ +newtype Raster u p (r :: Nat) (c :: Nat) a = Raster { _array :: Array u Ix2 a } + +-- | Warning: This will evaluate (at most) the 10x10 top-left corner of your +-- `Raster` for display. This should only be used for debugging. +instance (Show a, Load (EltRepr u Ix2) Ix2 a, Size u Ix2 a) => Show (Raster u p r c a) where + show (Raster a) = show . computeAs B $ extract' zeroIndex (r :. c) a + where (r :. c) = liftIndex (P.min 10) $ size a + +instance (Eq a, Unbox a) => Eq (Raster U p r c a) where + Raster a == Raster b = a == b + {-# INLINE (==) #-} + +instance (Eq a, Storable a) => Eq (Raster S p r c a) where + Raster a == Raster b = a == b + {-# INLINE (==) #-} + +instance (Eq a, Prim a) => Eq (Raster P p r c a) where + Raster a == Raster b = a == b + {-# INLINE (==) #-} + +instance (Eq a, NFData a) => Eq (Raster N p r c a) where + Raster a == Raster b = a == b + {-# INLINE (==) #-} + +instance Eq a => Eq (Raster B p r c a) where + Raster a == Raster b = a == b + {-# INLINE (==) #-} + +instance Eq a => Eq (Raster D p r c a) where + Raster a == Raster b = a == b + {-# INLINE (==) #-} + +instance Functor (Raster DW p r c) where + fmap f (Raster a) = Raster $ fmap f a + {-# INLINE fmap #-} + +instance Functor (Raster D p r c) where + fmap f (Raster a) = Raster $ fmap f a + {-# INLINE fmap #-} + +instance Functor (Raster DI p r c) where + fmap f (Raster a) = Raster $ fmap f a + {-# INLINE fmap #-} + +instance (KnownNat r, KnownNat c) => Applicative (Raster D p r c) where + pure = constant D Par + {-# INLINE pure #-} + + -- TODO: Use strict ($)? + fs <*> as = zipWith ($) fs as + {-# INLINE (<*>) #-} + +instance Semigroup a => Semigroup (Raster D p r c a) where + a <> b = zipWith (<>) a b + {-# INLINE (<>) #-} + +instance (Monoid a, KnownNat r, KnownNat c) => Monoid (Raster D p r c a) where + mempty = constant D Par mempty + {-# INLINE mempty #-} + + a `mappend` b = zipWith mappend a b + {-# INLINE mappend #-} + +instance (Num a, KnownNat r, KnownNat c) => Num (Raster D p r c a) where + a + b = zipWith (+) a b + {-# INLINE (+) #-} + + a - b = zipWith (-) a b + {-# INLINE (-) #-} + + a * b = zipWith (*) a b + {-# INLINE (*) #-} + + abs = fmap abs + {-# INLINE abs #-} + + signum = fmap signum + {-# INLINE signum #-} + + fromInteger = constant D Par . fromInteger + {-# INLINE fromInteger #-} + +instance (Fractional a, KnownNat r, KnownNat c) => Fractional (Raster D p r c a) where + a / b = zipWith (/) a b + {-# INLINE (/) #-} + + fromRational = constant D Par . fromRational + {-# INLINE fromRational #-} + +-- TODO: more explicit implementations? +-- | `length` has a specialized \(\mathcal{O}(1)\) implementation. +instance Foldable (Raster D p r c) where + foldMap f (Raster a) = foldMap f a + {-# INLINE foldMap #-} + + sum (Raster a) = A.sum a + {-# INLINE sum #-} + + product (Raster a) = A.product a + {-# INLINE product #-} + + -- | \(\mathcal{O}(1)\). + length (Raster a) = (\(r :. c) -> r * c) $ A.size a + {-# INLINE length #-} + +-- | \(\mathcal{O}(1)\). Force a `Raster`'s representation to `D`, allowing it +-- to undergo various operations. All operations between `D` `Raster`s are fused +-- and allocate no extra memory. +lazy :: Source u Ix2 a => Raster u p r c a -> Raster D p r c a +lazy (Raster a) = Raster $ delay a +{-# INLINE lazy #-} + +-- | Evaluate some lazy (`D`, `DW`, or `DI`) `Raster` to some explicit `Manifest` type +-- (i.e. to a real memory-backed Array). Will follow the `Comp`utation strategy +-- of the underlying 'massiv' `Array`. +-- +-- __Note:__ If using the `Par` computation strategy, make sure you're compiling with +-- @-with-rtsopts=-N@ to automatically use all available CPU cores at runtime. Otherwise +-- your "parallel" operations will only execute on one core. +strict :: (Load v Ix2 a, Mutable u Ix2 a) => u -> Raster v p r c a -> Raster u p r c a +strict u (Raster a) = Raster $ computeAs u a +{-# INLINE strict #-} + +-- | Create a `Raster` of any size which has the same value everywhere. +constant :: (KnownNat r, KnownNat c, Construct u Ix2 a) => u -> Comp -> a -> Raster u p r c a +constant u c a = fromFunction u c (const a) +{-# INLINE constant #-} + +-- | \(\mathcal{O}(1)\). Create a `Raster` from a function of its row and column number respectively. +fromFunction :: forall u p r c a. (KnownNat r, KnownNat c, Construct u Ix2 a) => + u -> Comp -> (Ix2 -> a) -> Raster u p r c a +fromFunction u c f = Raster $ makeArrayR u c sh f + where sh = fromInteger (natVal (Proxy :: Proxy r)) :. fromInteger (natVal (Proxy :: Proxy c)) +{-# INLINE fromFunction #-} + +-- | \(\mathcal{O}(1)\). Create a `Raster` from the values of any `GV.Vector` type. +-- Will fail if the size of the Vector does not match the declared size of the `Raster`. +fromVector :: forall v p r c a. (KnownNat r, KnownNat c, GV.Vector v a, Mutable (A.ARepr v) Ix2 a, Typeable v) => + Comp -> v a -> Either String (Raster (A.ARepr v) p r c a) +fromVector comp v | (r * c) == GV.length v = Right . Raster $ A.fromVector comp (r :. c) v + | otherwise = Left $ printf "Expected Pixel Count: %d - Actual: %d" (r * c) (GV.length v) + where r = fromInteger $ natVal (Proxy :: Proxy r) + c = fromInteger $ natVal (Proxy :: Proxy c) +{-# INLINE fromVector #-} + +-- | An RGBA image whose colour bands are distinct. +data RGBARaster p r c a = RGBARaster { _red :: !(Raster S p r c a) + , _green :: !(Raster S p r c a) + , _blue :: !(Raster S p r c a) + , _alpha :: !(Raster S p r c a) } + +-- | Read any image type into a `Raster` of distinct colour bands +-- with the cell type you declare. If the source image stores its +-- values as `Int` but you declare `Double`, the conversion will happen +-- automatically. +-- +-- Will fail if the declared size of the `Raster` +-- does not match the actual size of the input image. +fromRGBA :: forall p r c a. (Elevator a, KnownNat r, KnownNat c) => FilePath -> IO (Either String (RGBARaster p r c a)) +fromRGBA fp = do + cap <- getNumCapabilities + img <- setComp (bool Par Seq $ cap == 1) <$> readImageAuto fp + let rows = fromInteger $ natVal (Proxy :: Proxy r) + cols = fromInteger $ natVal (Proxy :: Proxy c) + (r :. c) = size img + if r == rows && c == cols + then do + (ar, ag, ab, aa) <- spreadRGBA img + pure . Right $ RGBARaster (Raster ar) (Raster ag) (Raster ab) (Raster aa) + else pure . Left $ printf "Expected Size: %d x %d - Actual Size: %d x %d" rows cols r c +{-# INLINE fromRGBA #-} + +spreadRGBA :: (Index ix, Elevator e) + => A.Array S ix (Pixel RGBA e) + -> IO (A.Array S ix e, A.Array S ix e, A.Array S ix e, A.Array S ix e) +spreadRGBA arr = do + let sz = A.size arr + mr <- A.unsafeNew sz + mb <- A.unsafeNew sz + mg <- A.unsafeNew sz + ma <- A.unsafeNew sz + flip A.imapP_ arr $ \ix (PixelRGBA r g b a) -> do + A.unsafeWrite mr ix r + A.unsafeWrite mg ix g + A.unsafeWrite mb ix b + A.unsafeWrite ma ix a + ar <- A.unsafeFreeze (getComp arr) mr + ag <- A.unsafeFreeze (getComp arr) mg + ab <- A.unsafeFreeze (getComp arr) mb + aa <- A.unsafeFreeze (getComp arr) ma + return (ar, ag, ab, aa) +{-# INLINE spreadRGBA #-} + +-- | Read a grayscale image. If the source file has more than one colour band, +-- they'll be combined automatically. +fromGray :: forall p r c a. (Elevator a, KnownNat r, KnownNat c) => FilePath -> IO (Either String (Raster S p r c a)) +fromGray fp = do + cap <- getNumCapabilities + img <- setComp (bool Par Seq $ cap == 1) <$> readImageAuto fp + let rows = fromInteger $ natVal (Proxy :: Proxy r) + cols = fromInteger $ natVal (Proxy :: Proxy c) + (r :. c) = size img + pure . bool (Left $ printf "Expected Size: %d x %d - Actual Size: %d x %d" rows cols r c) (Right $ f img) $ r == rows && c == cols + where f :: Image S Y a -> Raster S p r c a + f img = Raster . A.fromVector (getComp img) (size img) . VS.unsafeCast $ A.toVector img +{-# INLINE fromGray #-} + +-- | An invisible pixel (alpha channel set to 0). +invisible :: Pixel RGBA Word8 +invisible = PixelRGBA 0 0 0 0 + +-- | Construct a colour ramp. +-- ramp :: Ord k => [(Word8, Word8, Word8)] -> [k] -> M.Map k PixelRGBA8 +ramp :: Ord k => [(Word8, Word8, Word8)] -> [k] -> M.Map k (Pixel RGBA Word8) +ramp colours breaks = M.fromList . P.zip breaks $ P.map (\(r,g,b) -> PixelRGBA r g b maxBound) colours +{-# INLINE ramp #-} + +-- | From page 32 of /Cartographer's Toolkit/. +greenRed :: Ord k => [k] -> M.Map k (Pixel RGBA Word8) +greenRed = ramp colours + where colours = [ (0, 48, 0), (31, 79, 20), (100, 135, 68), (148, 193, 28), (193, 242, 3) + , (241, 255, 159), (249, 228, 227), (202, 145, 150), (153, 101, 97), (142, 38 ,18) ] + +-- | From page 33 of /Cartographer's Toolkit/. +spectrum :: Ord k => [k] -> M.Map k (Pixel RGBA Word8) +spectrum = ramp colours + where colours = [ (0, 22, 51), (51, 18, 135), (150, 0, 204), (242, 13, 177), (255, 61, 61) + , (240, 152, 56), (248, 230, 99), (166, 249, 159), (184, 249, 212), (216, 230, 253) ] + +-- | From page 34 of /Cartographer's Toolkit/. +blueGreen :: Ord k => [k] -> M.Map k (Pixel RGBA Word8) +blueGreen = ramp colours + where colours = [ (29, 43, 53), (37, 44, 95), (63, 70, 134), (89, 112, 147), (87, 124, 143) + , (117, 160, 125), (188, 219, 173), (239, 253, 163), (222, 214, 67), (189, 138, 55) ] + +-- | From page 35 of /Cartographer's Toolkit/. +purpleYellow :: Ord k => [k] -> M.Map k (Pixel RGBA Word8) +purpleYellow = ramp colours + where colours = [ (90, 89, 78), (73, 65, 132), (107, 86, 225), (225, 67, 94), (247, 55, 55) + , (251, 105, 46), (248, 174, 66), (249, 219, 25), (255, 255, 0), (242, 242, 242) ] + +-- | From page 36 of /Cartographer's Toolkit/. +brownBlue :: Ord k => [k] -> M.Map k (Pixel RGBA Word8) +brownBlue = ramp colours + where colours = [ (27, 36, 43), (86, 52, 42), (152, 107, 65), (182, 176, 152), (215, 206, 191) + , (198, 247, 0), (53, 227, 0), (30, 158, 184), (22, 109, 138), (12, 47, 122) ] + +-- | From page 37 of /Cartographer's Toolkit/. +grayBrown :: Ord k => [k] -> M.Map k (Pixel RGBA Word8) +grayBrown = ramp colours + where colours = [ (64, 57, 88), (95, 96, 116), (158, 158, 166), (206, 208, 197), (215, 206, 191) + , (186, 164, 150), (160, 124, 98), (117, 85, 72), (90, 70, 63), (39, 21, 17) ] + +-- | From page 38 of /Cartographer's Toolkit/. +greenPurple :: Ord k => [k] -> M.Map k (Pixel RGBA Word8) +greenPurple = ramp colours + where colours = [ (89, 168, 15), (158, 213, 76), (196, 237, 104), (226, 255, 158), (240, 242, 221) + , (248, 202, 140), (233, 161, 137), (212, 115, 132), (172, 67, 123), (140, 40, 110) ] + +-- | From page 39 of /Cartographer's Toolkit/. +brownYellow :: Ord k => [k] -> M.Map k (Pixel RGBA Word8) +brownYellow = ramp colours + where colours = [ (96, 72, 96), (120, 72, 96), (168, 96, 96), (192, 120, 96), (240, 168, 72) + , (248, 202, 140), (254, 236, 174), (255, 244, 194), (255, 247, 219), (255, 252, 246) ] + +-- | From page 40 of /Cartographer's Toolkit/. +purpleGreen :: Ord k => [k] -> M.Map k (Pixel RGBA Word8) +purpleGreen = ramp colours + where colours = [ (80, 73, 113), (117, 64, 152), (148, 116, 180), (199, 178, 214), (223, 204, 228) + , (218, 234, 193), (171, 214, 155), (109, 192, 103), (13, 177, 75), (57, 99, 83) ] + +-- | From page 41 of /Cartographer's Toolkit/. +purpleRed :: Ord k => [k] -> M.Map k (Pixel RGBA Word8) +purpleRed = ramp colours + where colours = [ (51, 60, 255), (76, 60, 233), (99, 60, 211), (121, 60, 188), (155, 60, 155) + , (166, 60, 143), (188, 60, 121), (206, 60, 94), (217, 60, 83), (255, 60, 76) ] + +-- | Convert a `Raster` into grayscale pixels, suitable for easy output with functions +-- like `writeImage`. +grayscale :: Functor (Raster u p r c) => Raster u p r c a -> Raster u p r c (Pixel Y a) +grayscale = fmap PixelY +{-# INLINE grayscale #-} + +-- | View a `Raster` as grayscale with the default image viewer of your OS. +-- +-- For more direct control, consider `displayImage` from 'massiv-io'. +display :: (Functor (Raster u p r c), Load u Ix2 (Pixel Y a), Elevator a) => Raster u p r c a -> IO () +display = displayImage . computeAs S . _array . grayscale + +-- | Render a PNG-encoded `BL.ByteString` from a coloured `Raster`. +-- Useful for returning a `Raster` from a webserver endpoint. +png :: (Source u Ix2 (Pixel cs a), ColorSpace cs a) => Raster u p r c (Pixel cs a) -> BL.ByteString +png (Raster a) = encode PNG def a +{-# INLINE png #-} + +-- | Called /LocalClassification/ in GaCM. The first argument is the value +-- to give to any index whose value is less than the lowest break in the `M.Map`. +-- +-- This is a glorified `fmap` operation, but we expose it for convenience. +classify :: (Ord a, Functor f) => b -> M.Map a b -> f a -> f b +classify d m r = fmap f r + where f a = maybe d snd $ M.lookupLE a m +{-# INLINE classify #-} + +-- | Finds the minimum value at each index between two `Raster`s. +lmin :: (Ord a, Source u Ix2 a) => Raster u p r c a -> Raster u p r c a -> Raster D p r c a +lmin = zipWith P.min +{-# INLINE lmin #-} + +-- | Finds the maximum value at each index between two `Raster`s. +lmax :: (Ord a, Source u Ix2 a) => Raster u p r c a -> Raster u p r c a -> Raster D p r c a +lmax = zipWith P.max +{-# INLINE lmax #-} + +-- | Averages the values per-index of all `Raster`s in a collection. +lmean :: (Real a, Fractional b, KnownNat r, KnownNat c) => NonEmpty (Raster D p r c a) -> Raster D p r c b +lmean (a :| [b]) = Raster $ A.zipWith (\n m -> realToFrac (n + m) / 2) (_array a) (_array b) +lmean (a :| [b,c]) = Raster $ A.zipWith3 (\n m o -> realToFrac (n + m + o) / 3) (_array a) (_array b) (_array c) +lmean (a :| as) = (\n -> realToFrac n / len) <$> foldl' (+) a as + where len = 1 + fromIntegral (length as) +{-# INLINE lmean #-} + +-- | The count of unique values at each shared index. +lvariety :: (KnownNat r, KnownNat c, Eq a) => NonEmpty (Raster D p r c a) -> Raster D p r c Word +lvariety = fmap (fromIntegral . length . NE.nub) . sequenceA +{-# INLINE lvariety #-} + +-- | The most frequently appearing value at each shared index. +lmajority :: (KnownNat r, KnownNat c, Ord a) => NonEmpty (Raster D p r c a) -> Raster D p r c a +lmajority = fmap majo . sequenceA +{-# INLINE lmajority #-} + +-- | Find the most common value in some `Foldable`. +majo :: (Foldable t, Ord a) => t a -> a +majo = fst . g . f + where f = foldl' (\m a -> M.insertWith (+) a 1 m) M.empty + g = L.foldl1' (\(a,c) (k,v) -> if c < v then (k,v) else (a,c)) . M.toList +{-# INLINE majo #-} + +-- | The least frequently appearing value at each shared index. +lminority :: (KnownNat r, KnownNat c, Ord a) => NonEmpty (Raster D p r c a) -> Raster D p r c a +lminority = fmap mino . sequenceA +{-# INLINE lminority #-} + +-- | Find the least common value in some `Foldable`. +mino :: (Foldable t, Ord a) => t a -> a +mino = fst . g . f + where f = foldl' (\m a -> M.insertWith (+) a 1 m) M.empty + g = L.foldl1' (\(a,c) (k,v) -> if c > v then (k,v) else (a,c)) . M.toList +{-# INLINE mino #-} + +-- | A measure of how spread out a dataset is. This calculation will fail +-- with `Nothing` if a length 1 list is given. +lvariance :: (Real a, KnownNat r, KnownNat c) => NonEmpty (Raster D p r c a) -> Maybe (Raster D p r c Double) +lvariance (_ :| []) = Nothing +lvariance rs = Just (f <$> sequenceA rs) + where len = realToFrac $ length rs + avg ns = (\z -> realToFrac z / len) $ foldl' (+) 0 ns + f os@(n :| ns) = foldl' (\acc m -> acc + ((realToFrac m - av) ^ 2)) ((realToFrac n - av) ^ 2) ns / (len - 1) + where av = avg os +{-# INLINE lvariance #-} + +-- Old implementation that was replaced with `sequenceA` usage above. I wonder if this is faster? +-- Leaving it here in case we feel like comparing later. +--listEm :: (Projection p, KnownNat r, KnownNat c) => NonEmpty (Raster p r c a) -> Raster p r c (NonEmpty a) +--listEm = sequenceA +--listEm (r :| rs) = foldl' (\acc s -> zipWith cons s acc) z rs +-- where z = (\a -> a :| []) <$> r +--{-# INLINE [2] listEm #-} + +-- | Combine two `Raster`s, element-wise, with a binary operator. +zipWith :: (Source u Ix2 a, Source u Ix2 b) => + (a -> b -> d) -> Raster u p r c a -> Raster u p r c b -> Raster D p r c d +zipWith f (Raster a) (Raster b) = Raster $ A.zipWith f a b +{-# INLINE zipWith #-} + +-- | Focal Addition. +fsum :: (Num a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a +fsum (Raster a) = Raster $ mapStencil (neighbourhoodStencil f (Fill 0)) a + where f nw no ne we fo ea sw so se = nw + no + ne + we + fo + ea + sw + so + se +{-# INLINE fsum #-} + +-- | Focal Product. +fproduct :: (Num a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a +fproduct (Raster a) = Raster $ mapStencil (neighbourhoodStencil f (Fill 1)) a + where f nw no ne we fo ea sw so se = nw * no * ne * we * fo * ea * sw * so * se + +-- | Focal Monoid - combine all elements of a neighbourhood via their `Monoid` instance. +-- In terms of precedence, the neighbourhood focus is the "left-most", and all other +-- elements are "added" to it. +-- +-- This is not mentioned in GaCM, but seems useful nonetheless. +fmonoid :: (Monoid a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a +fmonoid (Raster a) = Raster $ mapStencil (neighbourhoodStencil f (Fill mempty)) a + where f nw no ne we fo ea sw so se = fo `mappend` nw `mappend` no `mappend` ne `mappend` we `mappend` ea `mappend` sw `mappend` so `mappend` se + +-- | Focal Mean. +fmean :: (Real a, Fractional b, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c b +fmean = fmap (\n -> realToFrac n / 9) . fsum +{-# INLINE fmean #-} + +-- | Focal Maximum. +fmax :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a +fmax (Raster a) = Raster $ mapStencil (neighbourhoodStencil f Edge) a + where f nw no ne we fo ea sw so se = P.max nw . P.max no . P.max ne . P.max we . P.max fo . P.max ea . P.max sw $ P.max so se +{-# INLINE fmax #-} + +-- | Focal Minimum. +fmin :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a +fmin (Raster a) = Raster $ mapStencil (neighbourhoodStencil f Edge) a + where f nw no ne we fo ea sw so se = P.min nw . P.min no . P.min ne . P.min we . P.min fo . P.min ea . P.min sw $ P.min so se +{-# INLINE fmin #-} + +-- | Focal Variety - the number of unique values in each neighbourhood. +fvariety :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Word +fvariety (Raster a) = Raster $ mapStencil (neighbourhoodStencil f Edge) a + where f nw no ne we fo ea sw so se = fromIntegral . length $ L.nub [ nw, no, ne, we, fo, ea, sw, so, se ] +{-# INLINE fvariety #-} + +-- | Focal Majority - the most frequently appearing value in each neighbourhood. +fmajority :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a +fmajority (Raster a) = Raster $ mapStencil (neighbourhoodStencil f Continue) a + where f nw no ne we fo ea sw so se = majo [ nw, no, ne, we, fo, ea, sw, so, se ] +{-# INLINE fmajority #-} + +-- | Focal Minority - the least frequently appearing value in each neighbourhood. +fminority :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a +fminority (Raster a) = Raster $ mapStencil (neighbourhoodStencil f Continue) a + where f nw no ne we fo ea sw so se = mino [ nw, no, ne, we, fo, ea, sw, so, se ] +{-# INLINE fminority #-} + +-- | Focal Percentage - the percentage of neighbourhood values that are equal +-- to the neighbourhood focus. Not to be confused with `fpercentile`. +fpercentage :: (Eq a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Double +fpercentage (Raster a) = Raster $ mapStencil (neighbourhoodStencil f Continue) a + where f nw no ne we fo ea sw so se = ( bool 0 1 (nw == fo) + + bool 0 1 (no == fo) + + bool 0 1 (ne == fo) + + bool 0 1 (we == fo) + + bool 0 1 (ea == fo) + + bool 0 1 (sw == fo) + + bool 0 1 (so == fo) + + bool 0 1 (se == fo) ) / 8 +{-# INLINE fpercentage #-} + +-- | Focal Percentile - the percentage of neighbourhood values that are /less/ +-- than the neighbourhood focus. Not to be confused with `fpercentage`. +fpercentile :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Double +fpercentile (Raster a) = Raster $ mapStencil (neighbourhoodStencil f Continue) a + where f nw no ne we fo ea sw so se = ( bool 0 1 (nw < fo) + + bool 0 1 (no < fo) + + bool 0 1 (ne < fo) + + bool 0 1 (we < fo) + + bool 0 1 (ea < fo) + + bool 0 1 (sw < fo) + + bool 0 1 (so < fo) + + bool 0 1 (se < fo) ) / 8 +{-# INLINE fpercentile #-} + +-- `Fill def` has the highest chance of the edge pixel and the off-the-edge pixel +-- having a different value. This is until the following is addressed: +-- https://github.com/fosskers/mapalgebra/pull/3#issuecomment-379943208 +-- | Focal Linkage - a description of how each neighbourhood focus is connected +-- to its neighbours. Foci of equal value to none of their neighbours will have +-- a value of 0. +flinkage :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Line +flinkage (Raster a) = Raster $ mapStencil (neighbourhoodStencil linkage (Fill def)) a +{-# INLINE flinkage #-} + +-- | A description of lineal directions with the same encoding as `Drain`. +-- See `flinkage` and `flength`. +newtype Line = Line { _line :: Word8 } + deriving stock (Eq, Ord, Show) + deriving newtype (Storable, Prim, Default) + +instance NFData Line where + rnf (Line a) = deepseq a () + +linkage :: Eq a => a -> a -> a -> a -> a -> a -> a -> a -> a -> Line +linkage nw no ne we fo ea sw so se = Line $ axes + diags + where axes = bool 0 2 (no == fo) + bool 0 8 (we == fo) + bool 0 16 (ea == fo) + bool 0 64 (so == fo) + diags = bool 0 1 (nw == fo && not (testBit axes 1 || testBit axes 3)) + + bool 0 4 (ne == fo && not (testBit axes 1 || testBit axes 4)) + + bool 0 32 (sw == fo && not (testBit axes 3 || testBit axes 6)) + + bool 0 128 (se == fo && not (testBit axes 4 || testBit axes 6)) +{-# INLINE linkage #-} + +-- | Directions that neighbourhood foci can be connected by. +data Direction = East | NorthEast | North | NorthWest | West | SouthWest | South | SouthEast + deriving (Eq, Ord, Enum, Show) + +-- | Focal Length - the length of the lineal structure at every location. The result is in +-- "pixel units", where 1 is the height/width of one pixel. +flength :: Functor (Raster u p r c) => Raster u p r c Line -> Raster u p r c Double +flength = fmap f + where half = 1 / 2 + root = 1 / sqrt 2 + f (Line a) = bool 0 half (testBit a 1) + + bool 0 half (testBit a 3) + + bool 0 half (testBit a 4) + + bool 0 half (testBit a 6) + + bool 0 root (testBit a 0) + + bool 0 root (testBit a 2) + + bool 0 root (testBit a 5) + + bool 0 root (testBit a 7) +{-# INLINE flength #-} + +-- | A layout of the areal conditions of a single `Raster` pixel. +-- It describes whether each pixel corner is occupied by the same +-- "areal zone" as the pixel centre. +data Corners = Corners { _topLeft :: !Surround + , _bottomLeft :: !Surround + , _bottomRight :: !Surround + , _topRight :: !Surround } deriving (Eq, Show) + +instance NFData Corners where + rnf (Corners tl bl br tr) = tl `deepseq` bl `deepseq` br `deepseq` tr `deepseq` () + +-- | A state of surroundedness of a pixel corner. +-- For the examples below, the bottom-left pixel is considered the focus and +-- we're wondering about the surroundedness of its top-right corner. +data Surround = Complete -- ^ A corner has three of the same opponent against it. + -- + -- The corner is considered "occupied" by the opponent value, + -- thus forming a diagonal areal edge. + -- + -- @ + -- [ 1 1 ] + -- [ 0 1 ] + -- @ + | OneSide -- ^ One edge of a corner is touching an opponent, but + -- the other edge touches a friend. + -- + -- @ + -- [ 1 1 ] or [ 0 1 ] + -- [ 0 0 ] [ 0 1 ] + -- @ + | Open -- ^ A corner is surrounded by friends. + -- + -- @ + -- [ 0 0 ] or [ 0 0 ] or [ 1 0 ] + -- [ 0 0 ] [ 0 1 ] [ 0 0 ] + -- @ + | RightAngle -- ^ Similar to `Complete`, except that the diagonal + -- opponent doesn't match the other two. The corner + -- is considered surrounded, but not "occupied". + -- + -- @ + -- [ 1 2 ] + -- [ 0 1 ] + -- @ + | OutFlow -- ^ Similar to `Complete`, except that the area of the + -- focus surrounds the diagonal neighbour. + -- + -- @ + -- [ 0 1 ] + -- [ 0 0 ] + -- @ + deriving (Eq, Ord, Show) + +instance NFData Surround where + rnf s = case s of + Complete -> () + OneSide -> () + Open -> () + RightAngle -> () + OutFlow -> () + +-- | Imagining a 2x2 neighbourhood with its focus in the bottom-left, +-- what `Surround` relationship does the focus have with the other pixels? +surround :: Eq a => a -> a -> a -> a -> Surround +surround fo tl tr br + | up && tl == tr && tr == br = Complete + | up && right = RightAngle + | (up && diag) || (diag && right) = OneSide + | diag && fo == tl && fo == br = OutFlow + | otherwise = Open + where up = fo /= tl + diag = fo /= tr + right = fo /= br +{-# INLINE surround #-} + +-- | What is the total length of the areal edges (if there are any) at a given pixel? +frontage :: Corners -> Double +frontage (Corners tl bl br tr) = f tl + f bl + f br + f tr + where f Complete = 1 / sqrt 2 + f OneSide = 1 / 2 + f Open = 0 + f RightAngle = 1 + f OutFlow = 1 / sqrt 2 + +-- | Focal Partition - the areal form of each location, only considering +-- the top-right edge. +fpartition :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Corners +fpartition (Raster a) = Raster $ mapStencil partStencil a +{-# INLINE fpartition #-} + +partStencil :: (Eq a, Default a) => Stencil Ix2 a Corners +partStencil = makeStencil Reflect (2 :. 2) (1 :. 0) $ \f -> do + tl <- f (-1 :. 0) + tr <- f (-1 :. 1) + br <- f (0 :. 1) + fo <- f (0 :. 0) + pure $ Corners (surround fo tl tl fo) Open (surround fo fo br br) (surround fo tl tr br) +{-# INLINE partStencil #-} + +-- | Like `fpartition`, but considers the `Surround` of all corners. Is alluded to +-- in GaCM but isn't given its own operation name. +-- +-- If preparing for `ffrontage` or `farea`, you almost certainly want this function and +-- not `fpartition`. +fshape :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Corners +fshape (Raster a) = Raster $ mapStencil (neighbourhoodStencil f Reflect) a + where f nw no ne we fo ea sw so se = Corners (surround fo no nw we) + (surround fo so sw we) + (surround fo so se ea) + (surround fo no ne ea) +{-# INLINE fshape #-} + +-- | Focal Frontage - the length of areal edges between each pixel and its neighbourhood. +-- +-- Usually, the output of `fshape` is the appropriate input for this function. +ffrontage :: Functor (Raster u p r c) => Raster u p r c Corners -> Raster u p r c Double +ffrontage = fmap frontage +{-# INLINE ffrontage #-} + +-- | The area of a 1x1 square is 1. It has 8 right-triangular sections, +-- each with area 1/8. +area :: Corners -> Double +area (Corners tl bl br tr) = 1 - f tl - f bl - f br - f tr + where f Complete = 1/8 + f OutFlow = -1/8 -- For areas that "invade" their neighbours. + f _ = 0 +{-# INLINE area #-} + +-- | Focal Area - the area of the shape made up by a neighbourhood focus and its +-- surrounding pixels. Each pixel is assumed to have length and width of 1. +-- +-- Usually, the output of `fshape` is the appropriate input for this function. +farea :: Functor (Raster u p r c) => Raster u p r c Corners -> Raster u p r c Double +farea = fmap area +{-# INLINE farea #-} + +-- | Focal Volume - the surficial volume under each pixel, assuming the `Raster` +-- represents elevation in some way. +fvolume :: (Fractional a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a +fvolume (Raster a) = Raster $ mapStencil (neighbourhoodStencil volume Reflect) a +{-# INLINE fvolume #-} + +volume :: Fractional a => a -> a -> a -> a -> a -> a -> a -> a -> a -> a +volume tl up tr le fo ri bl bo br = + ((fo * 8) -- Simple algebra to reorganize individual volume calculations for each subtriangle. + + nw + no + + no + ne + + ne + ea + + ea + se + + se + so + + so + sw + + sw + we + + we + nw) / 24 + where nw = (tl + up + le + fo) / 4 + no = (up + fo) / 2 + ne = (up + tr + fo + ri) / 4 + we = (le + fo) / 2 + ea = (fo + ri) / 2 + sw = (le + fo + bl + bo) / 4 + so = (fo + bo) / 2 + se = (fo + ri + bo + br) / 4 +{-# INLINE volume #-} + +-- | Direct access to the entire neighbourhood. +neighbourhood :: Applicative f => (a -> a -> a -> a -> a -> a -> a -> a -> a -> b) -> (Ix2 -> f a) -> f b +neighbourhood g f = g <$> f (-1 :. -1) <*> f (-1 :. 0) <*> f (-1 :. 1) + <*> f (0 :. -1) <*> f (0 :. 0) <*> f (0 :. 1) + <*> f (1 :. -1) <*> f (1 :. 0) <*> f (1 :. 1) +{-# INLINE neighbourhood #-} + +neighbourhoodStencil :: Default a => (a -> a -> a -> a -> a -> a -> a -> a -> a -> b) -> Border a -> Stencil Ix2 a b +neighbourhoodStencil f b = makeStencil b (3 :. 3) (1 :. 1) (neighbourhood f) +{-# INLINE neighbourhoodStencil #-} + +-- | Get the surficial facets for each pixel and apply some function to them. +facetStencil :: (Fractional a, Default a) => (a -> a -> a -> a -> a -> a -> a -> a -> a -> b) -> Stencil Ix2 a b +facetStencil f = makeStencil Reflect (3 :. 3) (1 :. 1) (neighbourhood g) + where g nw no ne we fo ea sw so se = f ((nw + no + we + fo) / 4) + ((no + fo) / 2) + ((no + ne + fo + ea) / 4) + ((we + fo) / 2) + fo + ((fo + ea) / 2) + ((we + fo + sw + so) / 4) + ((fo + so) / 2) + ((fo + ea + so + se) / 4) +{-# INLINE facetStencil #-} + +-- | The first part to the "left pseudo inverse" needed to calculate +-- a best-fitting plane of 9 points. +leftPseudo :: LA.Matrix Double +leftPseudo = LA.inv (aT <> a) <> aT + where aT = LA.tr' a + a = LA.matrix 3 [ -0.5, -0.5, 1 + , -0.5, 0, 1 + , -0.5, 0.5, 1 + , 0, -0.5, 1 + , 0, 0, 1 + , 0, 0.5, 1 + , 0.5, -0.5, 1 + , 0.5, 0, 1 + , 0.5, 0.5, 1 ] + +-- TODO: newtype wrapper for `Radians`? +-- | Focal Gradient - a measurement of surficial slope for each pixel relative to +-- the horizonal cartographic plane. Results are in radians, with a flat plane +-- having a slope angle of 0 and a near-vertical plane approaching \(\tau / 4\). +fgradient :: (Manifest u Ix2 Double) => Raster u p r c Double -> Raster DW p r c Double +fgradient (Raster a) = Raster $ mapStencil (facetStencil gradient) a +{-# INLINE fgradient #-} + +-- | \(\tau\). One full rotation of the unit circle. +tau :: Double +tau = 6.283185307179586 + +-- | Given a list of \(z\) values, find the slope of the best-fit +-- plane that matches those points. +-- +-- See: https://stackoverflow.com/a/16669463/643684 +gradient :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double +gradient nw no ne we fo ea sw so se = (tau / 2) - acos (normal vs LA.! 2) + where vs = [ nw, no, ne, we, fo, ea, sw, so, se ] + +-- | Given a list of \(z\) values, find a normal vector that /points down/ +-- from a best-fit plane toward the cartographic plane. +normal :: [Double] -> LA.Vector Double +normal = LA.normalize . zcoord (-1) . normal' + +-- | A non-normalized, non-Z-corrected normal. Handy for `faspect`, +-- which needs to drop the Z and renormalize. +normal' :: [Double] -> LA.Vector Double +normal' vs = leftPseudo LA.#> LA.vector vs + +-- | Replace the Z-coordinate of a Vector. +zcoord :: Double -> LA.Vector Double -> LA.Vector Double +zcoord n v = LA.vector [ v LA.! 0, v LA.! 1, n ] + +-- | Focal Aspect - the compass direction toward which the surface +-- descends most rapidly. Results are in radians, with 0 or \(\tau\) being North, +-- \(\tau / 4\) being East, and so on. For areas that are essentially flat, their +-- aspect will be `Nothing`. +faspect :: Manifest u Ix2 Double => Raster u p r c Double -> Raster DW p r c (Maybe Double) +faspect (Raster a) = Raster $ mapStencil (facetStencil f) a + where f nw no ne we fo ea sw so se = case normal' [ nw, no, ne, we, fo, ea, sw, so, se ] of + n | ((n LA.! 0) =~ 0) && ((n LA.! 1) =~ 0) -> Nothing + | otherwise -> Just $ angle (LA.normalize $ zcoord 0 n) axis + axis = LA.vector [1, 0, 0] +{-# INLINE faspect #-} + +-- | Like `faspect`, but slightly faster. Beware of nonsense results when the plane is flat. +faspect' :: Manifest u Ix2 Double => Raster u p r c Double -> Raster DW p r c Double +faspect' (Raster a) = Raster $ mapStencil (facetStencil f) a + where f nw no ne we fo ea sw so se = angle (LA.normalize $ zcoord 0 $ normal' [ nw, no, ne, we, fo, ea, sw, so , se ]) axis + axis = LA.vector [1, 0, 0] +{-# INLINE faspect' #-} + +-- | Approximate Equality. Considers two `Double` to be equal if they are +-- less than \(/tau / 1024\) apart. +(=~) :: Double -> Double -> Bool +a =~ b = abs (a - b) < 0.0061359 + +-- | Given two normalized (length 1) vectors in R3, find the angle between them. +angle :: LA.Vector Double -> LA.Vector Double -> Double +angle u v = acos $ LA.dot u v + +-- | The main type for `fdownstream` and `fupstream`, used to calculate +-- Focal Drainage. This scheme for encoding drainage patterns is described +-- on page 81 of GaCM. +-- +-- ==== __Full Explanation__ +-- +-- Fluid can flow in or out of a square pixel in one of 256 ways. Imagine a pit, +-- whose neighbours are all higher in elevation: liquid would flow in from all +-- eight compass directions, but no liquid would flow out. Consider then +-- a neighbourhood of random heights - fluid might flow in or out of the focus +-- in any permutation of the eight directions. +-- +-- The scheme for encoding these permutations in a single `Word8` as described +-- in GaCM is this: +-- +-- Flow in a particular direction is represented by a power of 2: +-- +-- @ +-- [ 1 2 4 ] +-- [ 8 16 ] +-- [ 32 64 128 ] +-- @ +-- +-- Direction values are summed to make the encoding. +-- If there were drainage to the North, East, and SouthEast, we'd see a sum +-- of \(2 + 16 + 128 = 146\) to uniquely represent this. +-- +-- Analysing a drainage pattern from a `Drain` is just as easy: check if the bit corresponding +-- to the desired direction is flipped. The `direction` function handles this. +newtype Drain = Drain { _drain :: Word8 } + deriving stock (Eq, Ord, Show) + deriving newtype (Storable, Prim, Default) + +instance NFData Drain where + rnf (Drain a) = deepseq a () + +-- | Focal Drainage - downstream portion. This indicates the one or more compass +-- directions of steepest descent from each location. Appropriate as the input +-- to `fupstream`. +-- +-- __Note:__ Peak-like surfaces will not flow equally in all 8 directions. Consider this +-- neighbourhood: +-- +-- @ +-- [ 1 1 1 ] +-- [ 1 3 1 ] +-- [ 1 1 1 ] +-- @ +-- +-- According to the rules in GaCM for calculating the intermediate surficial "facet" +-- points for the focus, 3, we arrive at the following facet height matrix: +-- +-- @ +-- [ 1.5 2 1.5 ] +-- [ 2 3 2 ] +-- [ 1.5 2 1.5 ] +-- @ +-- +-- With these numbers it's clear that the corners would yield a steeper angle, +-- so our resulting `Drain` would only contain the directions +-- of the diagonals. +fdownstream :: Manifest u Ix2 Double => Raster u p r c Double -> Raster DW p r c Drain +fdownstream (Raster a) = Raster $ mapStencil (facetStencil downstream) a +{-# INLINE fdownstream #-} + +downstream :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Drain +downstream nw no ne we fo ea sw so se = Drain . snd $ foldl' f (0, 0) angles + where f (!curr, !s) (!a, !d) | a =~ curr = (curr, s + d) + | a > curr = (a, d) + | otherwise = (curr, s) + angles = [ (fo - nw, 1) + , (fo - no, 2) + , (fo - ne, 4) + , (fo - we, 8) + , (fo - ea, 16) + , (fo - sw, 32) + , (fo - so, 64) + , (fo - se, 128) ] + +-- | Focal Drainage - upstream portion. This indicates the one of more compass +-- directions from which liquid would flow into each surface location. +-- See also `fdownstream`. +fupstream :: Manifest u Ix2 Drain => Raster u p r c Drain -> Raster DW p r c Drain +fupstream (Raster a) = Raster $ mapStencil (neighbourhoodStencil f $ Fill (Drain 0)) a + where f nw no ne we _ ea sw so se = Drain $ bool 0 1 (testBit (_drain nw) 7) + + bool 0 2 (testBit (_drain no) 6) + + bool 0 4 (testBit (_drain ne) 5) + + bool 0 8 (testBit (_drain we) 4) + + bool 0 16 (testBit (_drain ea) 3) + + bool 0 32 (testBit (_drain sw) 2) + + bool 0 64 (testBit (_drain so) 1) + + bool 0 128 (testBit (_drain se) 0) +{-# INLINE fupstream #-} + +-- | Does a given `Drain` indicate flow in a certain `Direction`? +direction :: Direction -> Drain -> Bool +direction dir (Drain d) = case dir of + NorthWest -> testBit d 0 + North -> testBit d 1 + NorthEast -> testBit d 2 + West -> testBit d 3 + East -> testBit d 4 + SouthWest -> testBit d 5 + South -> testBit d 6 + SouthEast -> testBit d 7 + +-- | All `Direction`s that a `Drain` indicates flow toward. +directions :: Drain -> S.Set Direction +directions d = S.fromList $ foldl' (\acc dir -> bool acc (dir : acc) $ direction dir d) [] dirs + where dirs = [NorthWest, North, NorthEast, West, East, SouthWest, South, SouthEast] + +-- | The opposite of `directions`. +drainage :: S.Set Direction -> Drain +drainage = Drain . S.foldl' f 0 + where f acc d = case d of + NorthWest -> acc + 1 + North -> acc + 2 + NorthEast -> acc + 4 + West -> acc + 8 + East -> acc + 16 + SouthWest -> acc + 32 + South -> acc + 64 + SouthEast -> acc + 128 diff --git a/mapalgebra.cabal b/mapalgebra.cabal new file mode 100644 --- /dev/null +++ b/mapalgebra.cabal @@ -0,0 +1,93 @@ +-- This file has been generated from package.yaml by hpack version 0.20.0. +-- +-- see: https://github.com/sol/hpack +-- +-- hash: 8c5d8cd118654ac1b12534eabd521ebaabddb34bc070fb3780ac606cf3d2be9a + +name: mapalgebra +version: 0.1.0 +synopsis: Efficient, polymorphic Map Algebra. +description: Efficient, polymorphic Map Algebra. + . + This library is an implementation of /Map Algebra/ as described in the book /GIS and Cartographic Modeling/ by Dana Tomlin. The fundamental type, the `Raster`, is typesafe. Rasters of different size and projection are considered different types, and so cannot be combined in any way. + . + Also featured are op fusion (i.e. "lazy Rasters"), extremely fast Focal Operations, and typesafe NoData handling. Please see the main module for a more detailed introduction. +category: Geography +homepage: https://github.com/fosskers/mapalgebra +author: Colin Woodbury +maintainer: colin@fosskers.ca +copyright: 2018 Colin Woodbury +license: BSD3 +license-file: LICENSE +build-type: Simple +cabal-version: >= 1.10 + +extra-source-files: + CHANGELOG.md + data/1024x1024.tif + data/256x256.tif + data/512x512.tif + data/gray512.tif + README.md + +library + hs-source-dirs: + lib + ghc-options: -fwarn-unused-imports -fwarn-unused-binds -fwarn-name-shadowing -fwarn-unused-matches -fwarn-incomplete-patterns -Wincomplete-uni-patterns + build-depends: + base >=4.9 && <4.12 + , bytestring + , containers + , data-default + , deepseq + , hmatrix >=0.18 && <0.19 + , massiv >=0.1 && <0.2 + , massiv-io >=0.1 && <0.2 + , vector >=0.11 && <0.13 + exposed-modules: + Geography.MapAlgebra + default-language: Haskell2010 + +test-suite mapalgebra-test + type: exitcode-stdio-1.0 + main-is: Test.hs + hs-source-dirs: + test + ghc-options: -fwarn-unused-imports -fwarn-unused-binds -fwarn-name-shadowing -fwarn-unused-matches -fwarn-incomplete-patterns -Wincomplete-uni-patterns -threaded -with-rtsopts=-N + build-depends: + HUnit-approx >=1.1 && <1.2 + , QuickCheck + , base >=4.9 && <4.12 + , bytestring + , containers + , data-default + , deepseq + , hmatrix >=0.18 && <0.19 + , mapalgebra + , massiv >=0.1 && <0.2 + , massiv-io >=0.1 && <0.2 + , tasty >=0.11 && <2.0 + , tasty-hunit >=0.9 && <0.11 + , tasty-quickcheck >=0.8 && <0.10 + , vector >=0.11 && <0.13 + default-language: Haskell2010 + +benchmark mapalgebra-bench + type: exitcode-stdio-1.0 + main-is: Bench.hs + hs-source-dirs: + bench + ghc-options: -fwarn-unused-imports -fwarn-unused-binds -fwarn-name-shadowing -fwarn-unused-matches -fwarn-incomplete-patterns -Wincomplete-uni-patterns -threaded -O2 -with-rtsopts=-N + build-depends: + base >=4.9 && <4.12 + , bytestring + , containers + , criterion >=1.1 && <1.4 + , data-default + , deepseq + , hmatrix >=0.18 && <0.19 + , mapalgebra + , massiv >=0.1 && <0.2 + , massiv-io >=0.1 && <0.2 + , vector >=0.11 && <0.13 + default-language: Haskell2010 diff --git a/test/Test.hs b/test/Test.hs new file mode 100644 --- /dev/null +++ b/test/Test.hs @@ -0,0 +1,338 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE TypeApplications, ImplicitParams #-} + +module Main ( main ) where + +import Data.Int +import Data.List.NonEmpty (NonEmpty(..)) +import Data.Massiv.Array as A +import qualified Data.Set as S +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as U +import qualified Data.Vector.Storable as VS +import Data.Word +import Geography.MapAlgebra +import qualified Numeric.LinearAlgebra as LA +import Prelude as P +import Test.HUnit.Approx +import qualified Test.QuickCheck.Arbitrary as QC +import Test.Tasty +import Test.Tasty.HUnit +import Test.Tasty.QuickCheck + +--- + +main :: IO () +main = defaultMain suite + +suite :: TestTree +suite = testGroup "Unit Tests" + [ testGroup "Raster Creation" + [ testCase "constant (256x256)" $ length (lazy small) @?= 65536 + , testCase "constant (2^16 x 2^16)" $ length lazybig @?= 4294967296 + , testCase "Image Reading (RGBA)" $ do + i <- fileRGBA + fmap (getComp . _array . _red) i @?= Right Par + ] + , testGroup "Typeclass Ops" + [ testCase "(==)" $ assertBool "(==) doesn't work" (small == small) + , testCase "(+)" $ strict P (lazy one + lazy one) @?= two + ] + , testGroup "Folds" + [ testCase "sum (small)" $ P.sum (lazy small) @?= 327680 + -- , testCase "sum (large)" $ P.sum lazybig @?= 21474836480 + ] + , testGroup "Local Ops" + [ testCase "(+)" $ P.sum (lazy small + lazy small) @?= (327680 * 2) + , testCase "lmin" $ strict P (lmin one two) @?= one + , testCase "lvariety" $ (strict P . lvariety . fmap lazy $ one :| [two]) @?= two + , testCase "lmajority" $ (strict P . lmajority . fmap lazy $ one :| [one, two]) @?= one + , testCase "lminority" $ (strict P . lminority . fmap lazy $ one :| [one, two]) @?= two + -- , testCase "(+) big" $ strict P (lazy big + lazy big) @?= bog + ] + , testGroup "Focal Ops" + [ testCase "fvariety" $ strict P (fvariety one) @?= one + , testCase "fmax" $ strict P (fmax one) @?= one + , testCase "fmin" $ strict P (fmin one) @?= one + , testGroup "flinkage" + [ testCase "single point" singlePoint + , testCase "2x2 same" twoByTwoSame + , testCase "2x2 diff" twoByTwoDiff + , testCase "3x3" threeByThree + ] + , testCase "flength" flengthTest + , testCase "fpartition" fpartitionTest + , testCase "fshape" fshapeTest + , testCase "ffrontage" ffrontageTest + , testGroup "farea" + [ testCase "3x3 Open" fareaOpen + , testCase "3x3 Centre" fareaCentre + , testCase "4x4 Complex" fareaComplex + ] + , testGroup "fvolume" + [ testCase "3x3 Flat" fvolumeFlat + , testCase "3x3 Hill" fvolumeHill + ] + , testProperty "Least Squares" leastSquares + , testGroup "fgradient" + [ testCase "3x3 Flat" fgradientFlat + , testCase "3x3 (tau/8)" fgradient45 + ] + , testGroup "faspect" + [ testCase "3x3 Flat" faspectFlat + , testCase "3x3 East" faspectEast + , testCase "3x3 South" faspect45 + ] + , testGroup "fdownstream" + [ testCase "3x3 Spikey" fdownstream4 + , testCase "3x3 Flat" fdownstreamFlat + , testCase "3x3 Peak" fdownstreamPeak + , testCase "3x3 Pit" fdownstreamPit + ] + , testGroup "fupstream" + [ testCase "3x3 Peak" fupstreamPeak + , testCase "3x3 Flat" fupstreamFlat + ] + ] + ] + +one :: Raster P p 7 7 Word +one = constant P Seq 1 + +two :: Raster P p 7 7 Word +two = constant P Seq 2 + +small :: Raster P p 256 256 Int +small = constant P Seq 5 + +lazybig :: Raster D p 65536 65536 Int +lazybig = constant D Par 5 + +-- big :: Raster P p 65536 65536 Word8 +-- big = constant P Par 5 + +-- bog :: Raster P p 65536 65536 Word8 +-- bog = constant P Par 10 + +fileRGBA :: IO (Either String (RGBARaster p 512 512 Word8)) +fileRGBA = fromRGBA "data/512x512.tif" + +singlePoint :: Assertion +singlePoint = actual @?= expected + where expected :: Raster B p 1 1 Line + expected = constant B Seq (Line 0) + actual :: Raster B p 1 1 Line + actual = strict B . flinkage $ constant P Seq (1 :: Int) + +twoByTwoSame :: Assertion +twoByTwoSame = actual @?= expected + where expected :: Raster S p 2 2 Line + expected = fromRight . fromVector Seq . VS.fromList + $ P.map (Line . _drain . drainage . S.fromList) [ [ East, South ] + , [ West, South ] + , [ North,East ] + , [ West, North ] ] + actual :: Raster S p 2 2 Line + actual = fromRight . fmap (strict S . flinkage) . fromVector Seq $ U.fromList ([1,1,1,1] :: [Int]) + +twoByTwoDiff :: Assertion +twoByTwoDiff = actual @?= expected + where expected :: Raster S p 2 2 Line + expected = fromRight . fromVector Seq . VS.fromList + $ P.map (Line . _drain . drainage . S.fromList) [ [ SouthEast ] + , [ SouthWest ] + , [ NorthEast ] + , [ NorthWest ] ] + actual :: Raster S p 2 2 Line + actual = fromRight . fmap (strict S . flinkage) . fromVector Seq $ U.fromList ([1,2,2,1] :: [Int]) + +threeByThree :: Assertion +threeByThree = actual @?= expected + where expected :: Raster S p 3 3 Line + expected = fromRight . fromVector Seq . VS.fromList + $ P.map (Line . _drain . drainage . S.fromList) [ [ ] + , [ South ] + , [ ] + , [ East ] + , [ North, West, South, East ] + , [ West ] + , [ ] + , [ North ] + , [ ] ] + actual :: Raster S p 3 3 Line + actual = fromRight . fmap (strict S . flinkage) . fromVector Seq $ U.fromList ([1,2,1,2,2,2,1,2,1] :: [Int]) + +flengthTest :: Assertion +flengthTest = actual @?= expected + where actual :: Raster U p 3 3 Double + actual = strict U . flength . flinkage . fromRight . fromVector Seq $ VS.fromList ([1,2,1,2,2,2,1,2,1] :: [Int]) + expected :: Raster U p 3 3 Double + expected = fromRight . fromVector Seq $ U.fromList [ 0, 0.5, 0, 0.5, 2, 0.5, 0, 0.5, 0 ] + +fromRight :: Either a b -> b +fromRight (Right b) = b +fromRight _ = error "Was Left" + +fpartitionTest :: Assertion +fpartitionTest = actual @?= expected + where expected :: Raster B p 2 2 Corners + expected = fromRight . fromVector Seq $ V.fromList [ Corners Open Open Open Open + , Corners Open Open Open Open + , Corners OneSide Open OneSide Complete + , Corners Open Open Open Open ] + actual :: Raster B p 2 2 Corners + actual = strict B . fpartition . fromRight . fromVector Seq $ U.fromList ([1,1,2,1] :: [Int]) + +fshapeTest :: Assertion +fshapeTest = actual @?= expected + where expected :: Raster B p 3 3 Corners + expected = fromRight . fromVector Seq $ V.fromList [ Corners Open Open OutFlow Open + , Corners Open Open Open Open + , Corners Open OutFlow Open Open + , Corners Open Open Open Open + , Corners Complete Complete Complete Complete + , Corners Open Open Open Open + , Corners Open Open Open OutFlow + , Corners Open Open Open Open + , Corners OutFlow Open Open Open ] + actual :: Raster B p 3 3 Corners + actual = strict B . fshape . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,0,1,1,1,1] :: [Int]) + +ffrontageTest :: Assertion +ffrontageTest = let ?epsilon = 0.001 in actual @?~ expected + where expected :: Double + expected = 1 + (1 / sqrt 2) + actual :: Double + actual = flip index' (1 :. 1) . _array . strict S $ ffrontage rast + rast :: Raster DW p 4 4 Corners + rast = fshape . fromRight . fromVector Seq $ U.fromList ( [1,1,1,0 + ,1,0,0,0 + ,1,0,0,1 + ,1,0,1,1] :: [Int] ) + +fareaOpen :: Assertion +fareaOpen = actual @?= expected + where expected :: Raster U p 3 3 Double + expected = fromRight . fromVector Seq $ U.fromList [1,1,1,1,1,1,1,1,1] + actual :: Raster U p 3 3 Double + actual = strict U . farea . fshape . fromRight . fromVector Seq $ U.fromList ([0,0,0,0,0,0,0,0,0] :: [Int]) + +fareaCentre :: Assertion +fareaCentre = actual @?= expected + where expected :: Raster U p 3 3 Double + expected = fromRight . fromVector Seq $ U.fromList [ 1 + 1/8, 1, 1 + 1/8 + , 1, 1/2, 1 + , 1 + 1/8, 1, 1 + 1/8 ] + actual :: Raster U p 3 3 Double + actual = strict U . farea . fshape . fromRight . fromVector Seq $ U.fromList ([0,0,0,0,1,0,0,0,0] :: [Int]) + +fareaComplex :: Assertion +fareaComplex = let ?epsilon = 0.001 in actual @?~ (7 / 8) + where actual :: Double + actual = flip index' (1 :. 1) . _array . strict P $ farea rast + rast :: Raster DW p 4 4 Corners + rast = fshape . fromRight . fromVector Seq $ U.fromList ( [1,1,1,0 + ,1,0,0,0 + ,1,0,0,1 + ,1,0,1,1] :: [Int] ) + +fvolumeFlat :: Assertion +fvolumeFlat = (strict U $ fvolume expected) @?= expected + where expected :: Raster U p 3 3 Double + expected = fromRight . fromVector Seq $ U.fromList [8,8,8,8,8,8,8,8,8] + +fvolumeHill :: Assertion +fvolumeHill = (flip index' (1 :. 1) $ _array actual) @?= expected + where expected :: Double + expected = P.sum [20,20,16,20,16,16,16,16,12,16,12,12] / 12 + actual :: Raster U p 3 3 Double + actual = strict U . fvolume @Double . fromRight . fromVector Seq $ U.fromList [24,24,24 + ,16,16,16 + ,8,8,8] + +newtype Vec = Vec [Double] deriving (Show) + +instance Arbitrary Vec where + arbitrary = Vec <$> QC.vector 9 + +-- | A QuickCheck property to test whether my custom Least Squares is as +-- accurate as the one provided by HMatrix. +leastSquares :: Vec -> Bool +leastSquares (Vec vs) = f 0 && f 1 && f 2 + where m = head . LA.toColumns $ LA.linearSolveLS zing (LA.col vs) + v = leftPseudo LA.#> LA.vector vs + f i = (m LA.! i) =~ (v LA.! i) + +-- | Approximate Equality. +(=~) :: Double -> Double -> Bool +a =~ b = abs (a - b) < 0.0001 + +zing :: LA.Matrix Double +zing = LA.matrix 3 [ -0.5, -0.5, 1 + , -0.5, 0, 1 + , -0.5, 0.5, 1 + , 0, -0.5, 1 + , 0, 0, 1 + , 0, 0.5, 1 + , 0.5, -0.5, 1 + , 0.5, 0, 1 + , 0.5, 0.5, 1 ] + +fgradientFlat :: Assertion +fgradientFlat = actual @?= expected + where expected :: Raster U p 3 3 Double + expected = fromRight . fromVector Seq $ U.fromList [0,0,0,0,0,0,0,0,0] + actual :: Raster U p 3 3 Double + actual = strict U . fgradient . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,1,1,1,1,1] :: [Double]) + +fgradient45 :: Assertion +fgradient45 = let ?epsilon = 0.0001 in (flip index' (1 :. 1) $ _array actual) @?~ (tau / 8) + where actual :: Raster U p 3 3 Double + actual = strict U . fgradient . fromRight . fromVector Seq $ U.fromList ([3,3,3,2,2,2,1,1,1] :: [Double]) + +faspectFlat :: Assertion +faspectFlat = (flip index' (1 :. 1) $ _array actual) @?= Nothing + where actual :: Raster B p 3 3 (Maybe Double) + actual = strict B . faspect . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,1,1,1,1,1] :: [Double]) + +faspect45 :: Assertion +faspect45 = (flip index' (1 :. 1) $ _array actual) @?= Just (tau / 2) + where actual :: Raster B p 3 3 (Maybe Double) + actual = strict B . faspect . fromRight . fromVector Seq $ U.fromList ([3,3,3,2,2,2,1,1,1] :: [Double]) + +faspectEast :: Assertion +faspectEast = let ?epsilon = 0.0001 in (flip index' (1 :. 1) $ _array actual) @?~ (tau / 4) + where actual :: Raster B p 3 3 Double + actual = strict B . faspect' . fromRight . fromVector Seq $ U.fromList ([3,2,1,3,2,1,3,2,1] :: [Double]) + +fdownstream4 :: Assertion +fdownstream4 = (flip index' (1 :. 1) $ _array actual) @?= drainage (S.fromList [North,South,East,West]) + where actual :: Raster S p 3 3 Drain + actual = strict S . fdownstream . fromRight . fromVector Seq $ U.fromList ([3,1,3,1,2,1,3,1,3] :: [Double]) + +fdownstreamFlat :: Assertion +fdownstreamFlat = (flip index' (1 :. 1) $ _array actual) @?= drainage (S.fromList [East ..]) + where actual :: Raster S p 3 3 Drain + actual = strict S . fdownstream . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,1,1,1,1,1] :: [Double]) + +fdownstreamPeak :: Assertion +fdownstreamPeak = (flip index' (1 :. 1) $ _array actual) @?= drainage (S.fromList [NorthEast, NorthWest, SouthWest, SouthEast]) + where actual :: Raster S p 3 3 Drain + actual = strict S . fdownstream . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,3,1,1,1,1] :: [Double]) + +fdownstreamPit :: Assertion +fdownstreamPit = (flip index' (1 :. 1) $ _array actual) @?= Drain 0 + where actual :: Raster S p 3 3 Drain + actual = strict S . fdownstream . fromRight . fromVector Seq $ U.fromList ([2,2,2,2,1,2,2,2,2] :: [Double]) + +fupstreamFlat :: Assertion +fupstreamFlat = (flip index' (1 :. 1) $ _array actual) @?= drainage (S.fromList [East ..]) + where actual :: Raster S p 3 3 Drain + actual = strict S . fupstream . strict S . fdownstream . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,1,1,1,1,1] :: [Double]) + +fupstreamPeak :: Assertion +fupstreamPeak = (flip index' (1 :. 1) $ _array actual) @?= Drain 0 + where actual :: Raster S p 3 3 Drain + actual = strict S . fupstream . strict S . fdownstream . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,3,1,1,1,1] :: [Double])