Copyright  (c) Colin Woodbury 2018 

License  BSD3 
Maintainer  Colin Woodbury <colin@fosskers.ca> 
Safe Haskell  None 
Language  Haskell2010 
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:
 SingleRaster Local Operations >
fmap
with pure functions  MultiRaster Local Operations >
foldl
withzipWith
and pure functions  Focal Operations >
massiv
based smartStencil
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,
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 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, cachefriendly 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 withrtsopts=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.
 newtype Raster u p (r :: Nat) (c :: Nat) a = Raster {}
 lazy :: Source u Ix2 a => Raster u p r c a > Raster D p r c a
 strict :: (Load v Ix2 a, Mutable u Ix2 a) => u > Raster v p r c a > Raster u p r c a
 data RGBARaster p r c a = RGBARaster {}
 constant :: (KnownNat r, KnownNat c, Construct u Ix2 a) => u > Comp > a > Raster u p r c a
 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
 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)
 fromRGBA :: forall p r c a. (Elevator a, KnownNat r, KnownNat c) => FilePath > IO (Either String (RGBARaster p r c a))
 fromGray :: forall p r c a. (Elevator a, KnownNat r, KnownNat c) => FilePath > IO (Either String (Raster S p r c a))
 grayscale :: Functor (Raster u p r c) => Raster u p r c a > Raster u p r c (Pixel Y a)
 invisible :: Pixel RGBA Word8
 greenRed :: Ord k => [k] > Map k (Pixel RGBA Word8)
 spectrum :: Ord k => [k] > Map k (Pixel RGBA Word8)
 blueGreen :: Ord k => [k] > Map k (Pixel RGBA Word8)
 purpleYellow :: Ord k => [k] > Map k (Pixel RGBA Word8)
 brownBlue :: Ord k => [k] > Map k (Pixel RGBA Word8)
 grayBrown :: Ord k => [k] > Map k (Pixel RGBA Word8)
 greenPurple :: Ord k => [k] > Map k (Pixel RGBA Word8)
 brownYellow :: Ord k => [k] > Map k (Pixel RGBA Word8)
 purpleGreen :: Ord k => [k] > Map k (Pixel RGBA Word8)
 purpleRed :: Ord k => [k] > Map k (Pixel RGBA Word8)
 writeImage :: (Source r Ix2 (Pixel cs e), ColorSpace cs e) => FilePath > Image r cs e > IO ()
 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 ()
 png :: (Source u Ix2 (Pixel cs a), ColorSpace cs a) => Raster u p r c (Pixel cs a) > ByteString
 display :: (Functor (Raster u p r c), Load u Ix2 (Pixel Y a), Elevator a) => Raster u p r c a > IO ()
 class Projection p where
 reproject :: (Projection p, Projection r) => Point p > Point r
 data Sphere
 data LatLng
 data WebMercator
 data Point p = Point {}
 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
 classify :: (Ord a, Functor f) => b > Map a b > f a > f b
 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
 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
 lmean :: (Real a, Fractional b, KnownNat r, KnownNat c) => NonEmpty (Raster D p r c a) > Raster D p r c b
 lvariety :: (KnownNat r, KnownNat c, Eq a) => NonEmpty (Raster D p r c a) > Raster D p r c Word
 lmajority :: (KnownNat r, KnownNat c, Ord a) => NonEmpty (Raster D p r c a) > Raster D p r c a
 lminority :: (KnownNat r, KnownNat c, Ord a) => NonEmpty (Raster D p r c a) > Raster D p r c a
 lvariance :: (Real a, KnownNat r, KnownNat c) => NonEmpty (Raster D p r c a) > Maybe (Raster D p r c Double)
 fsum :: (Num a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c a
 fproduct :: (Num a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c a
 fmonoid :: (Monoid a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c a
 fmean :: (Real a, Fractional b, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c b
 fmax :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c a
 fmin :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c a
 fmajority :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c a
 fminority :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c a
 fvariety :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c Word
 fpercentage :: (Eq a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c Double
 fpercentile :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c Double
 newtype Line = Line {}
 flinkage :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c Line
 flength :: Functor (Raster u p r c) => Raster u p r c Line > Raster u p r c Double
 data Corners = Corners {
 _topLeft :: !Surround
 _bottomLeft :: !Surround
 _bottomRight :: !Surround
 _topRight :: !Surround
 data Surround
 = Complete
  OneSide
  Open
  RightAngle
  OutFlow
 fpartition :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c Corners
 fshape :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c Corners
 ffrontage :: Functor (Raster u p r c) => Raster u p r c Corners > Raster u p r c Double
 farea :: Functor (Raster u p r c) => Raster u p r c Corners > Raster u p r c Double
 newtype Drain = Drain {}
 data Direction
 direction :: Direction > Drain > Bool
 directions :: Drain > Set Direction
 drainage :: Set Direction > Drain
 fvolume :: (Fractional a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c a
 fgradient :: Manifest u Ix2 Double => Raster u p r c Double > Raster DW p r c Double
 faspect :: Manifest u Ix2 Double => Raster u p r c Double > Raster DW p r c (Maybe Double)
 faspect' :: Manifest u Ix2 Double => Raster u p r c Double > Raster DW p r c Double
 fdownstream :: Manifest u Ix2 Double => Raster u p r c Double > Raster DW p r c Drain
 fupstream :: Manifest u Ix2 Drain => Raster u p r c Drain > Raster DW p r c Drain
 leftPseudo :: Matrix Double
 tau :: Double
Types
Rasters
newtype Raster u p (r :: Nat) (c :: Nat) a Source #
A rectangular grid of data representing some area on the earth.
u
: What is the underlying representation of this Raster? (seemassiv
)p
: WhatProjection
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 samesize guarantee. myRaster :: Raster D WebMercator 256 256 Int myRaster = constant D Par 5 >>> length myRaster 65536
Functor (Raster k DI p r c) Source #  
Functor (Raster k D p r c) Source #  
Functor (Raster k DW p r c) Source #  
(KnownNat r, KnownNat c) => Applicative (Raster k D p r c) Source #  
Foldable (Raster k D p r c) Source # 

Eq a => Eq (Raster k D p r c a) Source #  
Eq a => Eq (Raster k B p r c a) Source #  
(Eq a, NFData a) => Eq (Raster k N p r c a) Source #  
(Eq a, Prim a) => Eq (Raster k P p r c a) Source #  
(Eq a, Storable a) => Eq (Raster k S p r c a) Source #  
(Eq a, Unbox a) => Eq (Raster k U p r c a) Source #  
(Fractional a, KnownNat r, KnownNat c) => Fractional (Raster k D p r c a) Source #  
(Num a, KnownNat r, KnownNat c) => Num (Raster k D p r c a) Source #  
(Show a, Load (EltRepr u Ix2) Ix2 a, Size u Ix2 a) => Show (Raster k u p r c a) Source #  Warning: This will evaluate (at most) the 10x10 topleft corner of your

Semigroup a => Semigroup (Raster k D p r c a) Source #  
(Monoid a, KnownNat r, KnownNat c) => Monoid (Raster k D p r c a) Source #  
strict :: (Load v Ix2 a, Mutable u Ix2 a) => u > Raster v p r c a > Raster u p r c a Source #
Evaluate some lazy (D
, DW
, or DI
) Raster
to some explicit Manifest
type
(i.e. to a real memorybacked 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
withrtsopts=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 Source #
An RGBA image whose colour bands are distinct.
Creation
constant :: (KnownNat r, KnownNat c, Construct u Ix2 a) => u > Comp > a > Raster u p r c a Source #
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 Source #
\(\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) Source #
fromRGBA :: forall p r c a. (Elevator a, KnownNat r, KnownNat c) => FilePath > IO (Either String (RGBARaster p r c a)) Source #
fromGray :: forall p r c a. (Elevator a, KnownNat r, KnownNat c) => FilePath > IO (Either String (Raster S p r c a)) Source #
Read a grayscale image. If the source file has more than one colour band, they'll be combined automatically.
Colouring
The Map
s here can be used with classify
to
transform ranges of values into certain colours in \(\mathcal{O}(n\log(n))\).
Each Mapgenerating 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 :: Functor (Raster u p r c) => Raster u p r c a > Raster u p r c (Pixel Y a) Source #
Convert a Raster
into grayscale pixels, suitable for easy output with functions
like writeImage
.
greenRed :: Ord k => [k] > Map k (Pixel RGBA Word8) Source #
From page 32 of Cartographer's Toolkit.
spectrum :: Ord k => [k] > Map k (Pixel RGBA Word8) Source #
From page 33 of Cartographer's Toolkit.
blueGreen :: Ord k => [k] > Map k (Pixel RGBA Word8) Source #
From page 34 of Cartographer's Toolkit.
purpleYellow :: Ord k => [k] > Map k (Pixel RGBA Word8) Source #
From page 35 of Cartographer's Toolkit.
brownBlue :: Ord k => [k] > Map k (Pixel RGBA Word8) Source #
From page 36 of Cartographer's Toolkit.
grayBrown :: Ord k => [k] > Map k (Pixel RGBA Word8) Source #
From page 37 of Cartographer's Toolkit.
greenPurple :: Ord k => [k] > Map k (Pixel RGBA Word8) Source #
From page 38 of Cartographer's Toolkit.
brownYellow :: Ord k => [k] > Map k (Pixel RGBA Word8) Source #
From page 39 of Cartographer's Toolkit.
purpleGreen :: Ord k => [k] > Map k (Pixel RGBA Word8) Source #
From page 40 of Cartographer's Toolkit.
purpleRed :: Ord k => [k] > Map k (Pixel RGBA Word8) Source #
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 (
being saved as Image
r RGBA
Double
)PNG
file would be written as
(
, thus using highest supported precision Image
r RGBA
Word16
)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 Source #
Render a PNGencoded 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 () Source #
View a Raster
as grayscale with the default image viewer of your OS.
For more direct control, consider displayImage
from 'massivio'.
Projections
class Projection p where Source #
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 tradeoffs. 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.
reproject :: (Projection p, Projection r) => Point p > Point r Source #
Reproject a Point
from one Projection
to another.
A perfect geometric sphere. The earth isn't actually shaped this way,
but it's a convenient middleground for converting between various
Projection
s.
data WebMercator Source #
The most commonly used Projection
for mapping in internet applications.
A location on the Earth in some Projection
.
Map Algebra
Local Operations
Operations between Raster
s. All operations are elementwise:
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 Source #
Combine two Raster
s, elementwise, 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
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 :: (Ord a, Source u Ix2 a) => Raster u p r c a > Raster u p r c a > Raster D p r c a Source #
Finds the maximum 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 Source #
Finds the minimum value at each index between two Raster
s.
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 Source #
Averages the values perindex of all Raster
s in a collection.
lvariety :: (KnownNat r, KnownNat c, Eq a) => NonEmpty (Raster D p r c a) > Raster D p r c Word Source #
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 Source #
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 Source #
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) Source #
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 memorybacked 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 Source #
Focal Addition.
fproduct :: (Num a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c a Source #
Focal Product.
fmonoid :: (Monoid a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c a Source #
Focal Monoid  combine all elements of a neighbourhood via their Monoid
instance.
In terms of precedence, the neighbourhood focus is the "leftmost", 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 Source #
Focal Mean.
fmax :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c a Source #
Focal Maximum.
fmin :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c a Source #
Focal Minimum.
fmajority :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c a Source #
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 Source #
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 Source #
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 Source #
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 Source #
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 linelike objects
in a Raster
. GaCM calls these lineal characteristics and describes them fully
on page 18 and 19.
flinkage :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c Line Source #
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 Source #
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.
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.
Corners  

A state of surroundedness of a pixel corner. For the examples below, the bottomleft pixel is considered the focus and we're wondering about the surroundedness of its topright corner.
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 [ 1 2 ] [ 0 1 ] 
OutFlow  Similar to [ 0 1 ] [ 0 0 ] 
fpartition :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c Corners Source #
Focal Partition  the areal form of each location, only considering the topright edge.
fshape :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c Corners Source #
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 Source #
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 Source #
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 Raster
s. GaCM calls elevation
features surficial characteristics and describes them fully on page 21
and 22.
Some of these operations require finding a "bestfit 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 bestfitting 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.
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.
Directions that neighbourhood foci can be connected by.
fvolume :: (Fractional a, Default a, Manifest u Ix2 a) => Raster u p r c a > Raster DW p r c a Source #
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 Source #
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 nearvertical plane approaching \(\tau / 4\).
faspect :: Manifest u Ix2 Double => Raster u p r c Double > Raster DW p r c (Maybe Double) Source #
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 Source #
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 Source #
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: Peaklike 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 Source #
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 Source #
The first part to the "left pseudo inverse" needed to calculate a bestfitting plane of 9 points.