Geography.MapAlgebra

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 Rasters (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, 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 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 Parallel 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.

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.

Types

Rasters

newtype Raster u p (r :: Nat) (c :: Nat) a

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

Constructors

Raster

Fields

• _array :: Array u Ix2 a

Instances

Functor (Raster k DI p r c)

Methods

fmap :: (a -> b) -> Raster k DI p r c a -> Raster k DI p r c b

(<$) :: a -> Raster k DI p r c b -> Raster k DI p r c a

Functor (Raster k D p r c)

Methods

fmap :: (a -> b) -> Raster k D p r c a -> Raster k D p r c b

(<$) :: a -> Raster k D p r c b -> Raster k D p r c a

Functor (Raster k DW p r c)

Methods

fmap :: (a -> b) -> Raster k DW p r c a -> Raster k DW p r c b

(<$) :: a -> Raster k DW p r c b -> Raster k DW p r c a

(KnownNat r, KnownNat c) => Applicative (Raster k D p r c)

Methods

pure :: a -> Raster k D p r c a

(<*>) :: Raster k D p r c (a -> b) -> Raster k D p r c a -> Raster k D p r c b

liftA2 :: (a -> b -> c) -> Raster k D p r c a -> Raster k D p r c b -> Raster k D p r c c

(*>) :: Raster k D p r c a -> Raster k D p r c b -> Raster k D p r c b

(<*) :: Raster k D p r c a -> Raster k D p r c b -> Raster k D p r c a

Foldable (Raster k D p r c)

length has a specialized \(\mathcal{O}(1)\) implementation.

Methods

fold :: Monoid m => Raster k D p r c m -> m

foldMap :: Monoid m => (a -> m) -> Raster k D p r c a -> m

foldr :: (a -> b -> b) -> b -> Raster k D p r c a -> b

foldr' :: (a -> b -> b) -> b -> Raster k D p r c a -> b

foldl :: (b -> a -> b) -> b -> Raster k D p r c a -> b

foldl' :: (b -> a -> b) -> b -> Raster k D p r c a -> b

foldr1 :: (a -> a -> a) -> Raster k D p r c a -> a

foldl1 :: (a -> a -> a) -> Raster k D p r c a -> a

toList :: Raster k D p r c a -> [a]

null :: Raster k D p r c a -> Bool

length :: Raster k D p r c a -> Int

elem :: Eq a => a -> Raster k D p r c a -> Bool

maximum :: Ord a => Raster k D p r c a -> a

minimum :: Ord a => Raster k D p r c a -> a

sum :: Num a => Raster k D p r c a -> a

product :: Num a => Raster k D p r c a -> a

Eq a => Eq (Raster k D p r c a)

Methods

(==) :: Raster k D p r c a -> Raster k D p r c a -> Bool

(/=) :: Raster k D p r c a -> Raster k D p r c a -> Bool

Eq a => Eq (Raster k B p r c a)

Methods

(==) :: Raster k B p r c a -> Raster k B p r c a -> Bool

(/=) :: Raster k B p r c a -> Raster k B p r c a -> Bool

(Eq a, NFData a) => Eq (Raster k N p r c a)

Methods

(==) :: Raster k N p r c a -> Raster k N p r c a -> Bool

(/=) :: Raster k N p r c a -> Raster k N p r c a -> Bool

(Eq a, Prim a) => Eq (Raster k P p r c a)

Methods

(==) :: Raster k P p r c a -> Raster k P p r c a -> Bool

(/=) :: Raster k P p r c a -> Raster k P p r c a -> Bool

(Eq a, Storable a) => Eq (Raster k S p r c a)

Methods

(==) :: Raster k S p r c a -> Raster k S p r c a -> Bool

(/=) :: Raster k S p r c a -> Raster k S p r c a -> Bool

(Eq a, Unbox a) => Eq (Raster k U p r c a)

Methods

(==) :: Raster k U p r c a -> Raster k U p r c a -> Bool

(/=) :: Raster k U p r c a -> Raster k U p r c a -> Bool

(Fractional a, KnownNat r, KnownNat c) => Fractional (Raster k D p r c a)

Methods

(/) :: Raster k D p r c a -> Raster k D p r c a -> Raster k D p r c a

recip :: Raster k D p r c a -> Raster k D p r c a

fromRational :: Rational -> Raster k D p r c a

(Num a, KnownNat r, KnownNat c) => Num (Raster k D p r c a)

Methods

(+) :: Raster k D p r c a -> Raster k D p r c a -> Raster k D p r c a

(-) :: Raster k D p r c a -> Raster k D p r c a -> Raster k D p r c a

(*) :: Raster k D p r c a -> Raster k D p r c a -> Raster k D p r c a

negate :: Raster k D p r c a -> Raster k D p r c a

abs :: Raster k D p r c a -> Raster k D p r c a

signum :: Raster k D p r c a -> Raster k D p r c a

fromInteger :: Integer -> Raster k D p r c a

(Show a, Load (EltRepr u Ix2) Ix2 a, Size u Ix2 a) => Show (Raster k u p r c a)

Warning: This will evaluate (at most) the 10x10 top-left corner of your Raster for display. This should only be used for debugging.

Methods

showsPrec :: Int -> Raster k u p r c a -> ShowS

show :: Raster k u p r c a -> String

showList :: [Raster k u p r c a] -> ShowS

Semigroup a => Semigroup (Raster k D p r c a)

Methods

(<>) :: Raster k D p r c a -> Raster k D p r c a -> Raster k D p r c a

sconcat :: NonEmpty (Raster k D p r c a) -> Raster k D p r c a

stimes :: Integral b => b -> Raster k D p r c a -> Raster k D p r c a

(Monoid a, KnownNat r, KnownNat c) => Monoid (Raster k D p r c a)

Methods

mempty :: Raster k D p r c a

mappend :: Raster k D p r c a -> Raster k D p r c a -> Raster k D p r c a

mconcat :: [Raster k D p r c a] -> Raster k D p r c a

lazy :: Source u Ix2 a => Raster u p r c a -> Raster D p r c a

\(\mathcal{O}(1)\). Force a Raster's representation to D, allowing it to undergo various operations. All operations between D Rasters are fused and allocate no extra memory.

strict :: (Load v Ix2 a, Mutable u Ix2 a) => u -> Raster v p r c a -> Raster u p r c a

Evaluate some lazy (D, DW, or DI) Raster to some explicit Manifest type (i.e. to a real memory-backed Array). Will follow the Computation 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.

data RGBARaster p r c a

An RGBA image whose colour bands are distinct.

Constructors

RGBARaster

Fields

• _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)

Creation

constant :: (KnownNat r, KnownNat c, Construct u Ix2 a) => u -> Comp -> a -> Raster u p r c a

Create a Raster of any size which has the same value everywhere.

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

\(\mathcal{O}(1)\). Create a Raster from a function of its row and column number respectively.

fromVector :: forall v p r c a. (KnownNat r, KnownNat c, Vector v a, Mutable (ARepr v) Ix2 a, Typeable v) => Comp -> v a -> Either String (Raster (ARepr v) p r c a)

\(\mathcal{O}(1)\). Create a Raster from the values of any Vector type. Will fail if the size of the Vector does not match the declared size of the Raster.

fromRGBA :: forall p r c a. (Elevator a, KnownNat r, KnownNat c) => FilePath -> IO (Either String (RGBARaster 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.

fromGray :: forall p r c a. (Elevator a, KnownNat r, KnownNat c) => FilePath -> IO (Either String (Raster S p r c a))

Read a grayscale image. If the source file has more than one colour band, they'll be combined automatically.

Colouring

The Maps 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 Rasters can be unwrapped with _array and then 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)

Convert a Raster into grayscale pixels, suitable for easy output with functions like writeImage.

invisible :: Pixel RGBA Word8

An invisible pixel (alpha channel set to 0).

greenRed :: Ord k => [k] -> Map k (Pixel RGBA Word8)

From page 32 of Cartographer's Toolkit.

spectrum :: Ord k => [k] -> Map k (Pixel RGBA Word8)

From page 33 of Cartographer's Toolkit.

blueGreen :: Ord k => [k] -> Map k (Pixel RGBA Word8)

From page 34 of Cartographer's Toolkit.

purpleYellow :: Ord k => [k] -> Map k (Pixel RGBA Word8)

From page 35 of Cartographer's Toolkit.

brownBlue :: Ord k => [k] -> Map k (Pixel RGBA Word8)

From page 36 of Cartographer's Toolkit.

grayBrown :: Ord k => [k] -> Map k (Pixel RGBA Word8)

From page 37 of Cartographer's Toolkit.

greenPurple :: Ord k => [k] -> Map k (Pixel RGBA Word8)

From page 38 of Cartographer's Toolkit.

brownYellow :: Ord k => [k] -> Map k (Pixel RGBA Word8)

From page 39 of Cartographer's Toolkit.

purpleGreen :: Ord k => [k] -> Map k (Pixel RGBA Word8)

From page 40 of Cartographer's Toolkit.

purpleRed :: Ord k => [k] -> Map k (Pixel RGBA Word8)

From page 41 of Cartographer's Toolkit.

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 :: (Source r Ix2 (Pixel cs e), ColorSpace cs e) => FilePath -> Image r cs e -> IO ()

This function will guess an output file format from the file extension and will write to file any image with the colorspace that is supported by that format. Precision of the image might be adjusted using Elevator if precision of the source array is not supported by the image file format. For instance an (Image r RGBA Double) being saved as PNG file would be written as (Image r RGBA Word16), thus using highest supported precision Word16 for that format. If automatic colors space is also desired, writeImageAuto can be used instead.

Can throw ConvertError, EncodeError and other usual IO errors.

writeImageAuto :: (Source r Ix2 (Pixel cs e), ColorSpace cs e, ToYA cs e, ToRGBA cs e, ToYCbCr cs e, ToCMYK cs e) => FilePath -> Image r cs e -> IO ()

Write an image to file while performing all necessary precisiona and color space conversions.

png :: (Source u Ix2 (Pixel cs a), ColorSpace cs a) => Raster u p r c (Pixel cs a) -> ByteString

Render a PNG-encoded ByteString from a coloured Raster. Useful for returning a Raster from a webserver endpoint.

display :: (Functor (Raster u p r c), Load u Ix2 (Pixel Y a), Elevator a) => Raster u p r c a -> IO ()

View a Raster as grayscale with the default image viewer of your OS.

For more direct control, consider displayImage from 'massiv-io'.

Projections

class Projection p where

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 Points between various Projections.

Note: Full support for Projections is still pending.

Minimal complete definition

toSphere, fromSphere

Methods

toSphere :: Point p -> Point Sphere

Convert a Point in this Projection to one of radians on a perfect Sphere.

fromSphere :: Point Sphere -> Point p

Convert a Point of radians on a perfect sphere to that of a specific Projection.

Instances

Projection * Sphere

Methods

toSphere :: Point Sphere p -> Point * Sphere

fromSphere :: Point * Sphere -> Point Sphere p

reproject :: (Projection p, Projection r) => Point p -> Point r

Reproject a Point from one Projection to another.

data Sphere

A perfect geometric sphere. The earth isn't actually shaped this way, but it's a convenient middle-ground for converting between various Projections.

Instances

Projection * Sphere

Methods

toSphere :: Point Sphere p -> Point * Sphere

fromSphere :: Point * Sphere -> Point Sphere p

data LatLng

Latitude (north-south position) and Longitude (east-west position).

data WebMercator

The most commonly used Projection for mapping in internet applications.

data Point p

A location on the Earth in some Projection.

Constructors

Point

Fields

• x :: !Double
• y :: !Double

Instances

Eq (Point k p)

Methods

(==) :: Point k p -> Point k p -> Bool

(/=) :: Point k p -> Point k p -> Bool

Show (Point k p)

Methods

showsPrec :: Int -> Point k p -> ShowS

show :: Point k p -> String

showList :: [Point k p] -> ShowS

Map Algebra

Local Operations

Operations between Rasters. 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 :: (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

Combine two Rasters, element-wise, with a binary operator.

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 :: (Ord a, Functor f) => b -> Map a b -> f a -> f b

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 Map.

This is a glorified fmap operation, but we expose it for convenience.

Binary

You can safely use these with the foldl family on any Foldable of Rasters. 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 :: (Ord a, Source u Ix2 a) => Raster u p r c a -> Raster u p r c a -> Raster D p r c a

Finds the maximum value at each index between two Rasters.

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

Finds the minimum value at each index between two Rasters.

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 :: (Real a, Fractional b, KnownNat r, KnownNat c) => NonEmpty (Raster D p r c a) -> Raster D p r c b

Averages the values per-index of all Rasters in a collection.

lvariety :: (KnownNat r, KnownNat c, Eq a) => NonEmpty (Raster D p r c a) -> Raster D p r c Word

The count of unique values at each shared index.

lmajority :: (KnownNat r, KnownNat c, Ord a) => NonEmpty (Raster D p r c a) -> Raster D p r c a

The most 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

The least frequently appearing value at each shared index.

lvariance :: (Real a, KnownNat r, KnownNat c) => NonEmpty (Raster D p r c a) -> Maybe (Raster D p r c Double)

A measure of how spread out a dataset is. This calculation will fail with Nothing if a length 1 list is given.

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 :: (Num a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a

Focal Addition.

fproduct :: (Num a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a

Focal Product.

fmonoid :: (Monoid a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a

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.

fmean :: (Real a, Fractional b, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c b

Focal Mean.

fmax :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a

Focal Maximum.

fmin :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a

Focal Minimum.

fmajority :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a

Focal Majority - the most 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

Focal Minority - the least frequently appearing value in each neighbourhood.

fvariety :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Word

Focal Variety - the number of unique values in each neighbourhood.

fpercentage :: (Eq a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Double

Focal Percentage - the percentage of neighbourhood values that are equal to the neighbourhood focus. Not to be confused with fpercentile.

fpercentile :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Double

Focal Percentile - the percentage of neighbourhood values that are less than the neighbourhood focus. Not to be confused with fpercentage.

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.

newtype Line

A description of lineal directions with the same encoding as Drain. See flinkage and flength.

Constructors

Line

Fields

• _line :: Word8

Instances

Eq Line

Methods

(==) :: Line -> Line -> Bool

(/=) :: Line -> Line -> Bool

Ord Line

Methods

compare :: Line -> Line -> Ordering

(<) :: Line -> Line -> Bool

(<=) :: Line -> Line -> Bool

(>) :: Line -> Line -> Bool

(>=) :: Line -> Line -> Bool

max :: Line -> Line -> Line

min :: Line -> Line -> Line

Show Line

Methods

showsPrec :: Int -> Line -> ShowS

show :: Line -> String

showList :: [Line] -> ShowS

Storable Line

Methods

sizeOf :: Line -> Int

alignment :: Line -> Int

peekElemOff :: Ptr Line -> Int -> IO Line

pokeElemOff :: Ptr Line -> Int -> Line -> IO ()

peekByteOff :: Ptr b -> Int -> IO Line

pokeByteOff :: Ptr b -> Int -> Line -> IO ()

peek :: Ptr Line -> IO Line

poke :: Ptr Line -> Line -> IO ()

Default Line

Methods

def :: Line

NFData Line

Methods

rnf :: Line -> ()

Prim Line

Methods

sizeOf# :: Line -> Int#

alignment# :: Line -> Int#

indexByteArray# :: ByteArray# -> Int# -> Line

readByteArray# :: MutableByteArray# s -> Int# -> State# s -> (#TupleRep [RuntimeRep], LiftedRep, State# s, Line#)

writeByteArray# :: MutableByteArray# s -> Int# -> Line -> State# s -> State# s

setByteArray# :: MutableByteArray# s -> Int# -> Int# -> Line -> State# s -> State# s

indexOffAddr# :: Addr# -> Int# -> Line

readOffAddr# :: Addr# -> Int# -> State# s -> (#TupleRep [RuntimeRep], LiftedRep, State# s, Line#)

writeOffAddr# :: Addr# -> Int# -> Line -> State# s -> State# s

setOffAddr# :: Addr# -> Int# -> Int# -> Line -> State# s -> State# s

flinkage :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Line

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.

flength :: Functor (Raster u p r c) => Raster u p r c Line -> Raster u p r c Double

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.

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.

data Corners

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.

Constructors

Corners

Fields

• _topLeft :: !Surround
• _bottomLeft :: !Surround
• _bottomRight :: !Surround
• _topRight :: !Surround

Instances

Eq Corners

Methods

(==) :: Corners -> Corners -> Bool

(/=) :: Corners -> Corners -> Bool

Show Corners

Methods

showsPrec :: Int -> Corners -> ShowS

show :: Corners -> String

showList :: [Corners] -> ShowS

NFData Corners

Methods

rnf :: Corners -> ()

data Surround

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.

Constructors

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 ]

Instances

Eq Surround

Methods

(==) :: Surround -> Surround -> Bool

(/=) :: Surround -> Surround -> Bool

Ord Surround

Methods

compare :: Surround -> Surround -> Ordering

(<) :: Surround -> Surround -> Bool

(<=) :: Surround -> Surround -> Bool

(>) :: Surround -> Surround -> Bool

(>=) :: Surround -> Surround -> Bool

max :: Surround -> Surround -> Surround

min :: Surround -> Surround -> Surround

Show Surround

Methods

showsPrec :: Int -> Surround -> ShowS

show :: Surround -> String

showList :: [Surround] -> ShowS

NFData Surround

Methods

rnf :: Surround -> ()

fpartition :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Corners

Focal Partition - the areal form of each location, only considering the top-right edge.

fshape :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Corners

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.

ffrontage :: Functor (Raster u p r c) => Raster u p r c Corners -> Raster u p r c Double

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.

farea :: Functor (Raster u p r c) => Raster u p r c Corners -> Raster u p r c Double

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.

Surficial

Focal operations that work over elevation Rasters. 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.

newtype Drain

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.

Constructors

Drain

Fields

• _drain :: Word8

Instances

Eq Drain

Methods

(==) :: Drain -> Drain -> Bool

(/=) :: Drain -> Drain -> Bool

Ord Drain

Methods

compare :: Drain -> Drain -> Ordering

(<) :: Drain -> Drain -> Bool

(<=) :: Drain -> Drain -> Bool

(>) :: Drain -> Drain -> Bool

(>=) :: Drain -> Drain -> Bool

max :: Drain -> Drain -> Drain

min :: Drain -> Drain -> Drain

Show Drain

Methods

showsPrec :: Int -> Drain -> ShowS

show :: Drain -> String

showList :: [Drain] -> ShowS

Storable Drain

Methods

sizeOf :: Drain -> Int

alignment :: Drain -> Int

peekElemOff :: Ptr Drain -> Int -> IO Drain

pokeElemOff :: Ptr Drain -> Int -> Drain -> IO ()

peekByteOff :: Ptr b -> Int -> IO Drain

pokeByteOff :: Ptr b -> Int -> Drain -> IO ()

peek :: Ptr Drain -> IO Drain

poke :: Ptr Drain -> Drain -> IO ()

Default Drain

Methods

def :: Drain

NFData Drain

Methods

rnf :: Drain -> ()

Prim Drain

Methods

sizeOf# :: Drain -> Int#

alignment# :: Drain -> Int#

indexByteArray# :: ByteArray# -> Int# -> Drain

readByteArray# :: MutableByteArray# s -> Int# -> State# s -> (#TupleRep [RuntimeRep], LiftedRep, State# s, Drain#)

writeByteArray# :: MutableByteArray# s -> Int# -> Drain -> State# s -> State# s

setByteArray# :: MutableByteArray# s -> Int# -> Int# -> Drain -> State# s -> State# s

indexOffAddr# :: Addr# -> Int# -> Drain

readOffAddr# :: Addr# -> Int# -> State# s -> (#TupleRep [RuntimeRep], LiftedRep, State# s, Drain#)

writeOffAddr# :: Addr# -> Int# -> Drain -> State# s -> State# s

setOffAddr# :: Addr# -> Int# -> Int# -> Drain -> State# s -> State# s

data Direction

Directions that neighbourhood foci can be connected by.

Constructors

East

NorthEast

North

NorthWest

West

SouthWest

South

SouthEast

Instances

Enum Direction

Methods

succ :: Direction -> Direction

pred :: Direction -> Direction

toEnum :: Int -> Direction

fromEnum :: Direction -> Int

enumFrom :: Direction -> [Direction]

enumFromThen :: Direction -> Direction -> [Direction]

enumFromTo :: Direction -> Direction -> [Direction]

enumFromThenTo :: Direction -> Direction -> Direction -> [Direction]

Eq Direction

Methods

(==) :: Direction -> Direction -> Bool

(/=) :: Direction -> Direction -> Bool

Ord Direction

Methods

compare :: Direction -> Direction -> Ordering

(<) :: Direction -> Direction -> Bool

(<=) :: Direction -> Direction -> Bool

(>) :: Direction -> Direction -> Bool

(>=) :: Direction -> Direction -> Bool

max :: Direction -> Direction -> Direction

min :: Direction -> Direction -> Direction

Show Direction

Methods

showsPrec :: Int -> Direction -> ShowS

show :: Direction -> String

showList :: [Direction] -> ShowS

direction :: Direction -> Drain -> Bool

Does a given Drain indicate flow in a certain Direction?

directions :: Drain -> Set Direction

All Directions that a Drain indicates flow toward.

drainage :: Set Direction -> Drain

The opposite of directions.

fvolume :: (Fractional a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a

Focal Volume - the surficial volume under each pixel, assuming the Raster represents elevation in some way.

fgradient :: Manifest u Ix2 Double => Raster u p r c Double -> Raster DW p r c Double

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\).

faspect :: Manifest u Ix2 Double => Raster u p r c Double -> Raster DW p r c (Maybe Double)

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 Double

Like faspect, but slightly faster. Beware of nonsense results when the plane is flat.

fdownstream :: Manifest u Ix2 Double => Raster u p r c Double -> Raster DW p r c Drain

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.

fupstream :: Manifest u Ix2 Drain => Raster u p r c Drain -> Raster DW p r c Drain

Focal Drainage - upstream portion. This indicates the one of more compass directions from which liquid would flow into each surface location. See also fdownstream.

Utilities

leftPseudo :: Matrix Double

The first part to the "left pseudo inverse" needed to calculate a best-fitting plane of 9 points.

tau :: Double

\(\tau\). One full rotation of the unit circle.