RichDEM C++ Reference¶
 template <class T>

class
A2Array2D
¶ Public Functions

A2Array2D
(std::string prefix, int per_tile_width, int per_tile_height, int width, int height, int cachesize)¶
 template <class U>

A2Array2D
(std::string filename_template, const A2Array2D<U> &other, int cachesize)¶

T &
operator()
(int32_t tx, int32_t ty, int32_t x, int32_t y)¶

T &
operator()
(int32_t x, int32_t y)¶

void
makeQuadIndex
(int32_t x, int32_t y, int32_t &tx, int32_t &ty, int32_t &px, int32_t &py) const¶

int32_t
width
() const¶

int32_t
height
() const¶

int32_t
widthInTiles
() const¶

int32_t
heightInTiles
() const¶

int32_t
notNullTiles
() const¶

int32_t
tileWidth
(int32_t tx, int32_t ty) const¶

int32_t
tileHeight
(int32_t tx, int32_t ty) const¶

int32_t
stdTileHeight
() const¶

int32_t
stdTileWidth
() const¶

void
setAll
(const T &val)¶

bool
isNoData
(int32_t x, int32_t y)¶

bool
isNoData
(int32_t tx, int32_t ty, int32_t px, int32_t py)¶

bool
isReadonly
() const¶

void
setNoData
(const T &ndval)¶

int32_t
getEvictions
() const¶

bool
isNullTile
(int32_t tx, int32_t ty) const¶

bool
isEdgeCell
(int32_t x, int32_t y) const¶

bool
in_grid
(int32_t x, int32_t y) const¶

bool
isInteriorCell
(int32_t x, int32_t y) const¶

void
printStamp
(int size)¶

void
loadTile
(int tx, int ty)¶
Private Functions

void
_LoadTile
(int tile_x, int tile_y)¶
Private Members

int
quick_width_in_tiles
¶

int
quick_height_in_tiles
¶

std::vector<std::vector<WrappedArray2D>>
data
¶

LRU<WrappedArray2D *>
lru
¶

int32_t
not_null_tiles
= 0¶

int64_t
total_width_in_cells
= 0¶

int64_t
total_height_in_cells
= 0¶

int32_t
per_tile_width
= 0¶

int32_t
per_tile_height
= 0¶

int32_t
evictions
= 0¶

int64_t
cells_in_not_null_tiles
= 0¶

T
no_data_to_set
¶

bool
readonly
= true¶
Friends

friend
richdem::A2Array2D::A2Array2D

 template <class T>

class
Array2D
¶  #include <Array2D.hpp>
Class to hold and manipulate GDAL and native rasters.
Array2D manages a twodimensional raster dataset. Passed a request to load such data, it peeks at the file header and can either load data on construction or wait until a later point. It can also offload data to disk.
 Author
 Richard Barnes (rbarnes@umn.edu)
Array2D permits simple copy construction as well as templated copies, which transfer projections and geotransforms, but not the actual data. This is useful for say, create a flow directions raster which is homologous to a DEM.
Array2D implements two addressing schemes: “xy” and “i”. All methods are available in each scheme; users may use whichever is convenient. The xyscheme accesses raster cells by their xycoordinates. The ischeme accesses cells by their address in a flat array. Internally, xyaddresses are converted to iaddresses. iaddressing is frequently faster because it reduces the space needed to store coordinates and requires no addressing mathematics; however, xyaddressing may be more intuitive. It is suggested to develop algorithms using xyaddressing and then convert them to iaddressing if additional speed is desired. The results of the two versions can then be compared against each other to verify that using iaddressing has not introduced any errors.
Unnamed Group
Public Types

typedef int32_t
xy_t
¶ xyaddressing data type

typedef uint32_t
i_t
¶ iaddressing data type
Public Functions

Array2D
()¶

Array2D
(xy_t width, xy_t height, const T &val = T())¶ Creates a raster of the specified dimensions.
 Parameters
width
: Width of the rasterheight
: Height of the rasterval
: Initial value of all the raster’s cells. Defaults to the Array2D template type’s default value

Array2D
(T *data0, const xy_t width, const xy_t height)¶ Wraps a flat array in an Array2D object.
Wraps a flat array in an Array2D object. The Array2D does not take ownership of the data.
 Parameters
data0
: Pointer to data to wrapwidth
: Width of the dataheight
: Height of the data
 template <class U>

Array2D
(const Array2D<U> &other, const T &val = T())¶ Create a raster with the same properties and dimensions as another raster. No data is copied between the two.
 Parameters
other
: Raster whose properties and dimensions should be copiedval
: Initial value of all the raster’s cells.

Array2D
(const std::string &filename, bool native, xy_t xOffset = 0, xy_t yOffset = 0, xy_t part_width = 0, xy_t part_height = 0, bool exact = false, bool load_data = true)¶ TODO.

void
dumpData
()¶ Caches the raster data and all its properties to disk. Data is then purged from RAM.
 Post
 Calls to loadData() after this will result in data being loaded from the cache.

void
loadData
()¶ Loads data from disk into RAM.
If dumpData() has been previously called, data is loaded from the cache; otherwise, it is loaded from a GDAL file. No data is loaded if data is already present in RAM.

T *
getData
()¶ Returns a pointer to the internal data array.

bool
empty
() const¶ Returns TRUE if no data is present in RAM.

T
noData
() const¶ Returns the NoData value of the raster. Cells equal to this value sould generally not be used in calculations. But note that the isNoData() method is a much better choice for testing whether a cell is NoData or not.

T
min
() const¶ Finds the minimum value of the raster, ignoring NoData cells.

T
max
() const¶ Finds the maximum value of the raster, ignoring NoData cells.

void
replace
(const T oldval, const T newval)¶ Replace one cell value with another throughout the raster. Can operate on NoData cells.
 Parameters
oldval
: Value to be replacednewval
: Value to replace ‘oldval’ with

i_t
countval
(const T val) const¶ Counts the number of occurrences of a particular value in the raster. Can operate on NoData cells.
 Return
 The number of times ‘val’ appears in the raster. Will be 0 if raster is not loaded in RAM.
 Parameters
val
: Value to be be counted

void
iToxy
(const i_t i, xy_t &x, xy_t &y) const¶ Convert from index coordinates to x,y coordinates.
 Parameters
i
: Index coordinatex
: Xcoordinate of iy
: Ycoordinate of i

i_t
xyToI
(xy_t x, xy_t y) const¶ Convert from x,y coordinates to index coordinates.
 Return
 Returns the index coordinate i of (x,y)
 Parameters
x
: Xcoordinate to converty
: Ycoordinate to convert

i_t
nToI
(i_t i, xy_t dx, xy_t dy) const¶ Given a cell identified by an icoordinate, return the icoordinate of the neighbour identified by dx,dy.
 Return
 icoordinate of the neighbour. Usually referred to as ‘ni’
 Parameters
i
: icoordinate of cell whose neighbour needs to be identifieddx
: xdisplacement of the neighbour from idy
: ydisplacement of the neighbour from i

i_t
getN
(i_t i, uint8_t n) const¶ Given a cell identified by an icoordinate, return the icoordinate of the neighbour identified by n.
 Return
 icoordinate of the neighbour. Usually referred to as ‘ni’
 Parameters
i
: icoordinate of cell whose neighbour needs to be identifiedn
: Neighbour to be identified

int
nshift
(const uint8_t n) const¶ Return the offset of the neighbour cell identified by n.
 Return
 Offset of the neighbour n
 Parameters
n
: Neighbour for which offset should be retrieved

bool
operator==
(const Array2D<T> &o) const¶ Determine if two rasters are equivalent based on dimensions, NoData value, and their data.

bool
isNoData
(xy_t x, xy_t y) const¶ Whether or not a cell is NoData using x,y coordinates.
 Return
 Returns TRUE if the cell is NoData
 Parameters
x
: Xcoordinate of cell to testy
: Ycoordinate of cell to test

bool
isNoData
(i_t i) const¶ Whether or not a cell is NoData using i coordinates.
 Return
 Returns TRUE if the cell is NoData
 Parameters
i
: icoordinate of cell to test

void
flipVert
()¶ Flips the raster from top to bottom.

void
flipHorz
()¶ Flips the raster from sidetoside.

void
transpose
()¶ Flips the raster about its diagonal axis, like a matrix tranpose.

bool
inGrid
(xy_t x, xy_t y) const¶ Test whether a cell lies within the boundaries of the raster.
 Return
 TRUE if cell lies within the raster
 Parameters
x
: Xcoordinate of cell to testy
: Ycoordinate of cell to test

bool
isEdgeCell
(xy_t x, xy_t y) const¶ Test whether a cell lies on the boundary of the raster.
 Return
 TRUE if cell lies on the raster’s boundary
 Parameters
x
: Xcoordinate of cell to testy
: Xcoordinate of cell to test

bool
isTopLeft
(xy_t x, xy_t y) const¶ Determines whether an (x,y) pair is the top left of the DEM.
 Return
 True, if the (x,y) pair is the top left of the DEM; otherwise, false

bool
isTopRight
(xy_t x, xy_t y) const¶ Determines whether an (x,y) pair is the top right of the DEM.
 Return
 True, if the (x,y) pair is the top right of the DEM; otherwise, false

bool
isBottomLeft
(xy_t x, xy_t y) const¶ Determines whether an (x,y) pair is the bottom left of the DEM.
 Return
 True, if the (x,y) pair is the bottom left of the DEM; otherwise, false

bool
isBottomRight
(xy_t x, xy_t y) const¶ Determines whether an (x,y) pair is the bottom right of the DEM.
 Return
 True, if the (x,y) pair is the bottom right of the DEM; otherwise, false

bool
isTopRow
(xy_t x, xy_t y) const¶ Determines whether an (x,y) pair is in the top row of the DEM.
 Return
 True, if the (x,y) pair is in the top row of the DEM; otherwise, false

bool
isBottomRow
(xy_t x, xy_t y) const¶ Determines whether an (x,y) pair is in the bottom row of the DEM.
 Return
 True, if the (x,y) pair is in the bottom row of the DEM; otherwise, false

bool
isLeftCol
(xy_t x, xy_t y) const¶ Determines whether an (x,y) pair is in the left column of the DEM.
 Return
 True, if the (x,y) pair is in the left column of the DEM; otherwise, false

bool
isRightCol
(xy_t x, xy_t y) const¶ Determines whether an (x,y) pair is in the right column of the DEM.
 Return
 True, if the (x,y) pair is in the right column of the DEM; otherwise, false

bool
isEdgeCell
(i_t i) const¶ Test whether a cell lies on the boundary of the raster.
 Return
 TRUE if cell lies on the raster’s boundary
 Parameters
i
: icoordinate of cell to test

void
setNoData
(const T &ndval)¶ Sets the NoData value of the raster.
 Parameters
ndval
: Value to change NoData to

void
setAll
(const T val)¶ Sets all of the raster’s cells to ‘val’.
 Parameters
val
: Value to change the cells to

void
resize
(const xy_t width0, const xy_t height0, const T &val0 = T())¶ Resize the raster. Note: this clears all the raster’s data.
 Parameters
width0
: New width of the rasterheight0
: New height of the rasterval0
: Value to set all the cells to. Defaults to the raster’s template type default value

void
expand
(xy_t new_width, xy_t new_height, const T val)¶ Makes a raster larger and retains the raster’s old data, similar to resize.
Note: Using this command requires RAM equal to the sum of the old raster and the new raster. The old raster is placed in the upperleft of the new raster.
 Parameters
new_width
: New width of the raster. Must be >= the old width.new_height
: New height of the raster. Must be >= the old height.val
: Value to set the new cells to

void
countDataCells
() const¶ Counts the number of cells which are not NoData.

i_t
numDataCells
() const¶ Returns the number of cells which are not NoData. May count them.
 Return
 Returns the number of cells which are not NoData.

T &
operator()
(i_t i)¶ Return cell value based on icoordinate.
 Return
 The value of the cell identified by ‘i’
 Parameters
i
: icoordinate of cell whose data should be fetched.

T
operator()
(i_t i) const¶ Return cell value based on icoordinate.
 Return
 The value of the cell identified by ‘i’
 Parameters
i
: icoordinate of cell whose data should be fetched.

T &
operator()
(xy_t x, xy_t y)¶ Return cell value based on x,y coordinates.
 Return
 The value of the cell identified by x,y
 Parameters
x
: Xcoordinate of cell whose data should be fetched.y
: Ycoordinate of cell whose data should be fetched.

T
operator()
(xy_t x, xy_t y) const¶ Return cell value based on x,y coordinates.
 Return
 The value of the cell identified by x,y
 Parameters
x
: Xcoordinate of cell whose data should be fetched.y
: Ycoordinate of cell whose data should be fetched.

std::vector<T>
topRow
() const¶ Returns a copy of the top row of the raster.
 Return
 A vector containing a copy of the top row of the raster

std::vector<T>
bottomRow
() const¶ Returns a copy of the bottom row of the raster.
 Return
 A vector containing a copy of the bottom row of the raster

std::vector<T>
leftColumn
() const¶ Returns a copy of the left column of the raster.
Top to bottom is reoriented as left to right.
 Return
 A vector containing a copy of the left column of the raster

std::vector<T>
rightColumn
() const¶ Returns a copy of the right column of the raster.
Top to bottom is reoriented as left to right.
 Return
 A vector containing a copy of the right column of the raster

void
setRow
(xy_t y, const T &val)¶ Sets an entire row of a raster to a given value.
 Parameters
y
: The row to be setval
: The value to set the row to

void
setCol
(xy_t x, const T &val)¶ Sets an entire column of a raster to a given value.
 Parameters
x
: The column to be setval
: The value to set the column to

std::vector<T>
getRowData
(xy_t y) const¶ Returns a copy of an arbitrary row of the raster.
 Return
 A vector containing a copy of the selected row
 Parameters
y
: The row to retrieve

std::vector<T>
getColData
(xy_t x) const¶ Returns a copy of an arbitrary column of the raster.
 Return
 A vector containing a copy of the selected column
 Parameters
x
: The column to retrieve

void
clear
()¶ Clears all raster data from RAM.
 template <class U>

void
templateCopy
(const Array2D<U> &other)¶ Copies the geotransform, projection, and basename of another raster.
 Parameters
other
: Raster to copy from

void
printStamp
(int size, std::string msg = "") const¶ Output a square of cells useful for determining raster orientation.
This method prints out a square block of cells whose upperleft corner is the (integerdivision) center of the raster.
Stamps are only shown if the SHOW_STAMPS preprocessor variable is set.
Since algorithms may have to flip rasters horizontally or vertically before manipulating them, it is important that all algorithms work on data in the same orientation. This method, used in testing, helps a user ensure that their algorithm is orientating data correctly.
 Parameters
size
: Output stamp will be size x sizemsg
: Message to print prior to the stamp

void
printBlock
(const int radius, const xy_t x0, const xy_t y0, bool color = false, const std::string msg = "") const¶ Prints a square of cells centered at x,y. Useful for debugging.
 Parameters
radius
: Output stamp will be 2*radius x 2*radiusx0
: Xcoordinate of block centery0
: Ycoordinate of block centercolor
: Print the (x,y) cell in colour?msg
: Optional message to print above the block

void
printAll
(const std::string msg = "") const¶ Prints the entire array.
 Parameters
msg
: Optional message to print above the block

double
getCellArea
() const¶ Get the area of an individual cell in square projection units.
 Return
 The area of the cell in square projection units

double
getCellLengthX
() const¶ Get the length of a cell along the raster’s horizontal axis.
 Return
 The length of the cell along the raster’s horizontal axis

double
getCellLengthY
() const¶ Get the length of a cell along the raster’s horizontal axis.
 Return
 The length of the cell along the raster’s horizontal axis

void
scale
(const double x)¶ Multiplies the entire array by a scalar.
 Parameters
x
: Value to multiply array by

bool
owned
() const¶
Public Members
Private Functions

void
saveToCache
(const std::string &filename)¶ Saves raster to a simplystructure file on disk, possibly using compression.
 Post
 Using loadData() after running this function will result in data being loaded from the cache, rather than the original file (if any).
Private Members

ManagedVector<T>
data
¶ Holds the raster data in a 1D array this improves caching versus a 2D array

T
no_data
¶ NoData value of the raster.

bool
from_cache
¶ If TRUE, loadData() loads data from the cache assuming the Native format. Otherwise, it assumes it is loading from a GDAL file.
Friends

friend
richdem::Array2D::Array2D

class
GridCell
¶  #include <grid_cell.hpp>
Stores the (x,y) coordinates of a grid cell.
Subclassed by richdem::GridCellZ< elev_t >, richdem::GridCellZ< double >, richdem::GridCellZ< float >
 template <class elev_t>

class
GridCellZ
¶  #include <grid_cell.hpp>
Stores the (x,y,z) coordinates of a grid cell; useful for priority sorting with GridCellZ_pq.
Inherits from richdem::GridCell
Subclassed by richdem::GridCellZk< elev_t >
Public Members

elev_t
z
¶ Grid cell’s zcoordinate.

elev_t
 template <>

template<>
classGridCellZ
<double>¶  #include <grid_cell.hpp>
An (x,y,z) cell with NaNs taken as infinitely small numbers.
Inherits from richdem::GridCell
Public Members

double
z
¶ Grid cell’s zcoordinate.

double
 template <>

template<>
classGridCellZ
<float>¶  #include <grid_cell.hpp>
An (x,y,z) cell with NaNs taken as infinitely small numbers.
Inherits from richdem::GridCell
Public Functions

GridCellZ
(int x, int y, float z)¶

GridCellZ
()
Public Members

float
z
¶ Grid cell’s zcoordinate.

 template <class elev_t>

class
GridCellZk
¶  #include <grid_cell.hpp>
Stores the (x,y,z) coordinates of a grid cell and a priority indicator k; used by GridCellZk_pq.
Inherits from richdem::GridCellZ< elev_t >
Public Functions

GridCellZk
(int x, int y, elev_t z, int k)¶

GridCellZk
()¶

bool
operator<
(const GridCellZk<elev_t> &a) const¶

bool
operator>
(const GridCellZk<elev_t> &a) const¶
Public Members

int
k
¶ Used to store an integer to make sorting stable.

 template <typename T>

class
GridCellZk_pq
¶  #include <grid_cell.hpp>
A priority queue of GridCellZk, sorted by ascending height or, if heights are equal, by the order of insertion.
Inherits from std::priority_queue< GridCellZk< T >, std::vector< GridCellZk< T > >, std::greater< GridCellZk< T > > >
Private Members

uint64_t
count
= 0¶

uint64_t

class
LayoutfileReader
¶  #include <Layoutfile.hpp>
Used for reading a layoutfile describing a tiled dataset.
The class acts as a generator. The layoutfile is read on construction and its contents retrieved with next(). The Layoutfile specification can be found in Layoutfile.hpp.
Public Functions

LayoutfileReader
(std::string layout_filename)¶ Construct a new LayoutfileReader object reading from a given file.
 Author
 Richard Barnes
 Parameters
layout_filename
: Layoutfile to read from.

bool
next
()¶ Advance the reader to the next layoutfile entry.
 Return
 True if reader advanced successfully; false if not.

bool
newRow
() const¶  Return
 True if the current entry is the beginning of a new row.

const std::string &
getFilename
() const¶  Return
 The current entry’s filename without the path (e.g. “file.ext”).

const std::string &
getBasename
() const¶  Return
 The current entry’s filename without the path or extension (e.g. “file”).

const std::string
getFullPath
() const¶  Return
 The current entry’s path + filename (e.g. “path/to/file.ext”).

const std::string
getGridLocName
() const¶ Return a string representation of the current entry’s coordinates.
A layoutfile is a 2D grid of file names. This method returns the current entry’s position in that grid as
<X>_<Y>
 Return
 Current entry’s position as a string of the form <X>_<Y>

const std::string &
getPath
() const¶  Return
 Path of layoutfile: of “path/to/layoutfile.layout” returns “path/to/”.

bool
isNullTile
() const¶  Return
 True if the current entry was a blank

int
getX
() const¶  Return
 Xcoordinate of the current entry.

int
getY
() const¶  Return
 Ycoordinate of the current entry.


class
LayoutfileWriter
¶  #include <Layoutfile.hpp>
Used for creating a layoutfile describing a tiled dataset.
The class acts as an inverse generator. The layoutfile is created on construction and its contents appended to with addEntry(). The Layoutfile specification can be found in Layoutfile.hpp.
 template <class T>

class
LRU
¶  #include <lru.hpp>
A LeastRecently Used (LRU) cache.
Public Functions

void
insert
(const T &entry)¶ Insert an item into the LRU.
The item is either added to the queue or its entry is moved to the top of the queue. If the item is new and the length of the queue is greater than maxlen, then the least recently seen item is evicted from the queue.
 Parameters
entry
: The item to add to the queue.

int
size
() const¶ Returns the number of itmes in the LRU cache.
 Return
 Number of items in the LRU cache

int
getCapacity
() const¶ Returns the capacity of the LRU cache.
 Return
 The capacity of the LRU cache

void
 template <class T>

class
ManagedVector
¶  #include <ManagedVector.hpp>
ManagedVector works like a regular vector, but can wrap external memory.
Public Functions

ManagedVector
()¶ Creates an empty ManagedVector.

ManagedVector
(size_t size, T default_val = T())¶ Creates a ManagedVector with
size
members each set todefault_val
 Parameters
size
: Number of elements to be created in the vectordefault_val
: Initial value of the elements

ManagedVector
(T *data, size_t size)¶ Creates a ManagedVector which wraps
data0
of lengthsize0
 Parameters
data
: Memory to wrapsize
: Number of elements to wrap
 template <class U>

ManagedVector
(const ManagedVector<U> &other)¶

ManagedVector
(const ManagedVector<T> &other)¶
 template <class U>

ManagedVector
(ManagedVector<U> &&other)¶

~ManagedVector
()¶
 template <class U>

ManagedVector<T> &
operator=
(const ManagedVector<U> &other)¶
 template <class U>

ManagedVector<T> &
operator=
(ManagedVector<U> &&other)¶ Move assignment operator

T *
data
()¶ Get a raw pointer to the managed data
 Return
 A raw pointer to the managed data

const T *
data
() const¶ Get a raw constant pointer to the managed data
 Return
 A raw constant pointer to the managed data

bool
empty
() const¶ Are there more than zero elements being managed?
 Return
 True, if zero elements are managed; otherwise, false

size_t
size
() const¶ Get the number of elements being managed
 Return
 The number of elements being managed

bool
owned
() const¶ Determine whether the ManagedVector owns the memory it is managing
 Return
 True, if this ManagedVector owns its memory; otherwise, false

void
resize
(size_t new_size)¶

T &
operator[]
(size_t i)¶

const T &
operator[]
(size_t i) const¶
Private Members

bool
_owned
= true¶ If this is true, we are responsible for cleanup of the data.

size_t
_size
= 0¶ Number of elements being managed.
Friends

friend
richdem::ManagedVector::ManagedVector


class
ProgressBar
¶  #include <ProgressBar.hpp>
Manages a consolebased progress bar to keep the user entertained.
Defining the global
NOPROGRESS
will disable all progress operations, potentially speeding up a program. The look of the progress bar is shown in ProgressBar.hpp.Public Functions

void
start
(uint32_t total_work)¶ Start/reset the progress bar.
 Parameters
total_work
: The amount of work to be completed, usually specified in cells.

void
update
(uint32_t work_done0)¶ Update the visible progress bar, but only if enough work has been done.
Define the global
NOPROGRESS
flag to prevent this from having an effect. Doing so may speed up the program’s execution.

ProgressBar &
operator++
()¶ Increment by one the work done and update the progress bar.

double
stop
()¶ Stop the progress bar. Throws an exception if it wasn’t started.
 Return
 The number of seconds the progress bar was running.

double
time_it_took
()¶  Return
 Return the time the progress bar ran for.

uint32_t
cellsProcessed
() const¶
Private Functions

void
clearConsoleLine
() const¶ Clear current line on console so a new progress bar can be written.
Private Members

uint32_t
total_work
¶ Total work to be accomplished.

uint32_t
next_update
¶ Next point to update the visible progress bar.

uint32_t
call_diff
¶ Interval between updates in work units.

uint32_t
work_done
¶

uint16_t
old_percent
¶ Old percentage value (aka: should we update the progress bar) TODO: Maybe that we do not need this.

void

class
StreamLogger
¶ Public Functions

~StreamLogger
()¶
 template <typename T>

StreamLogger &
operator<<
(const T &t)¶

StreamLogger &
operator<<
(std::ostream &(*f)(std::ostream&))¶


class
TA_Setup_Curves_Vars
¶

class
TA_Setup_Vars
¶  #include <terrain_attributes.hpp>
Calculate a variety of terrain attributes.
This calculates a variety of terrain attributes according to the work of Burrough 1998’s “Principles of Geographical Information Systems” (p. 190). Algorithms and helpful ArcGIS pages are noted in comments in the code.
 Author
 Richard Barnes (rbarnes@umn.edu), Burrough (1998)
 Pre
 This function should never be called on a NoData cell
 Parameters
&elevations
: An elevation gridx0
: X coordinate of cell to perform calculation ony0
: Y coordinate of cell to perform calculation on&rise_over_run
: Returns riseoverrun slope as per Horn 1981&aspect
: Returns aspect as per Horn 1981 in degrees [0,360). Degrees increase in a clockwise fashion. 0 is north, 1 indicates a flat surface.&curvature
: Returns the difference of profile and planform curvatures (TODO: Clarify, this is from ArcGIS and was poorly described)&profile_curvature
: Returns the profile curvature as per Zevenbergen and Thorne 1987. 0 indicates a flat surface.&planform_curvature
: Returns the planform curvature as per Zevenbergen and Thorne 1987. 0 indicates a flat surface.

class
Timer
¶  #include <timer.hpp>
Used to time how intervals in code.
Such as how long it takes a given function to run, or how long I/O has taken.
Public Functions

void
start
()¶ Start the timers. Throws an exception if timer was already running.

double
stop
()¶ Stop the timer. Throws an exception if timer was already stopped. Calling this adds to the timer’s accumulated time.
 Return
 The accumulated time in seconds.

double
accumulated
()¶ Returns the timer’s accumulated time. Throws an exception if the timer is running.
 Return
 The timer’s accumulated time, in seconds.

double
lap
()¶ Returns the time between when the timer was started and the current moment. Throws an exception if the timer is not running.
 Return
 Time since the timer was started and current moment, in seconds.

void
reset
()¶ Stops the timer and resets its accumulated time. No exceptions are thrown ever.
Private Functions

double
timediff
(timeval beginning, timeval end)¶ Number of (fractional) seconds between two time objects.

void

class
WrappedArray2D
¶ Inherits from richdem::Array2D< T >
Public Functions

template<>
voidlazySetAll
()¶

template<>

namespace
richdem
¶ Typedefs

typedef uint8_t
d8_flowdir_t
¶

using
richdem::GridCellZ_pq = typedef std::priority_queue<GridCellZ<elev_t>, std::vector<GridCellZ<elev_t> >, std::greater<GridCellZ<elev_t> > >
A priority queue of GridCellZ, sorted by ascending height.

typedef char
label_t
¶
Enums
Functions

static std::string
trimStr
(std::string const &str)¶ Eliminate spaces from the beginning and end of str.

static std::string
GetBaseName
(std::string filename)¶ Get only the filename without its extension. That is, convert “path/to/file.ext” to “file”

void
ProcessMemUsage
(long &vmpeak, long &vmhwm)¶ Return memory statistics of the process.
This code is drawn from “http://stackoverflow.com/a/671389/752843”
 Parameters
vmpeak
: Peak virtual memory size (kB)vmhwm
: Peak resident set size (kB)

our_random_engine &
rand_engine
()¶

void
seed_rand
(unsigned long seed)¶

int
uniform_rand_int
(int from, int thru)¶

double
uniform_rand_real
(double from, double thru)¶

double
normal_rand
(double mean, double stddev)¶

RandomEngineState
SaveRandomState
()¶

void
SetRandomState
(const RandomEngineState &res)¶
 template <class T>

T
uniform_bits
()¶

std::string
PrintRichdemHeader
(int argc, char **argv)¶ Takes the program’s command line arguments and prints to stdout a header with a variety of useful information for identifying the particulars of what was run.
 template <class elev_t>

bool
HasDepressions
(const Array2D<elev_t> &elevations)¶ Determine if a DEM has depressions.
PriorityFlood starts on the edges of the DEM and then works its way inwards using a priority queue to determine the lowest cell which has a path to the edge. The neighbours of this cell are added to the priority queue. If the neighbours are lower than the cell which is adding them, then they are part of a depression and the question is answered.
 Author
 Richard Barnes (rbarnes@umn.edu)
 Pre
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM. Note that the NoData value is assumed to be a negative number less than any actual data value.
 Return
 True if the DEM contains depressions; otherwise, false.
 Correctness:
 The correctness of this command is determined by inspection. (TODO)
 Parameters
&elevations
: A grid of cell elevations
 template <class elev_t>

void
PriorityFlood_Original
(Array2D<elev_t> &elevations)¶ Fills all pits and removes all digital dams from a DEM.
PriorityFlood starts on the edges of the DEM and then works its way inwards using a priority queue to determine the lowest cell which has a path to the edge. The neighbours of this cell are added to the priority queue. If the neighbours are lower than the cell which is adding them, then they are raised to match its elevation; this fills depressions.
 Author
 Richard Barnes (rbarnes@umn.edu)
 Pre
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM. Note that the NoData value is assumed to be a negative number less than any actual data value.
 Post
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM.
 elevations contains no landscape depressions or digital dams.
 Correctness:
 The correctness of this command is determined by inspection. (TODO)
 Parameters
&elevations
: A grid of cell elevations
 template <class elev_t>

void
PriorityFlood_Barnes2014
(Array2D<elev_t> &elevations)¶ Fills all pits and removes all digital dams from a DEM, but faster.
PriorityFlood starts on the edges of the DEM and then works its way inwards using a priority queue to determine the lowest cell which has a path to the edge. The neighbours of this cell are added to the priority queue if they are higher. If they are lower, they are raised to the elevation of the cell adding them, thereby filling in pits. The neighbors are then added to a “pit” queue which is used to flood pits. Cells which are higher than a pit being filled are added to the priority queue. In this way, pits are filled without incurring the expense of the priority queue.
 Author
 Richard Barnes (rbarnes@umn.edu)
 Pre
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM. Note that the NoData value is assumed to be a negative number less than any actual data value.
 Post
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM.
 elevations contains no landscape depressions or digital dams.
 Correctness:
 The correctness of this command is determined by inspection. (TODO)
 Parameters
&elevations
: A grid of cell elevations
 template <class elev_t>

void
PriorityFloodEpsilon_Barnes2014
(Array2D<elev_t> &elevations)¶ Modifies floatingpoint cell elevations to guarantee drainage.
This version of PriorityFlood starts on the edges of the DEM and then works its way inwards using a priority queue to determine the lowest cell which has a path to the edge. The neighbours of this cell are added to the priority queue if they are higher. If they are lower, then their elevation is increased by a small amount to ensure that they have a drainage path and they are added to a “pit” queue which is used to flood pits. Cells which are higher than a pit being filled are added to the priority queue. In this way, pits are filled without incurring the expense of the priority queue.
 Author
 Richard Barnes (rbarnes@umn.edu)
 Pre
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM. Note that the NoData value is assumed to be a negative number less than any actual data value.
 Post
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM.
 elevations has no landscape depressions, digital dams, or flats.
 Correctness:
 The correctness of this command is determined by inspection. (TODO)
 Parameters
&elevations
: A grid of cell elevations
 template <>

void
PriorityFloodEpsilon_Barnes2014
(Array2D<uint8_t> &elevations)¶ PriorityFlood+Epsilon is not available for integer data types.
 template <>

void
PriorityFloodEpsilon_Barnes2014
(Array2D<uint16_t> &elevations)¶ PriorityFlood+Epsilon is not available for integer data types.
 template <>

void
PriorityFloodEpsilon_Barnes2014
(Array2D<int16_t> &elevations)¶ PriorityFlood+Epsilon is not available for integer data types.
 template <>

void
PriorityFloodEpsilon_Barnes2014
(Array2D<uint32_t> &elevations)¶ PriorityFlood+Epsilon is not available for integer data types.
 template <>

void
PriorityFloodEpsilon_Barnes2014
(Array2D<int32_t> &elevations)¶ PriorityFlood+Epsilon is not available for integer data types.
 template <class elev_t>

void
PriorityFloodFlowdirs_Barnes2014
(const Array2D<elev_t> &elevations, Array2D<d8_flowdir_t> &flowdirs)¶ Determines D8 flow directions and implicitly fills pits.
This version of PriorityFlood starts on the edges of the DEM and then works its way inwards using a priority queue to determine the lowest cell which has a path to the edge. The neighbours of this cell are given D8 flow directions which point to it. All depressions are implicitly filled and digital dams removed.
 Author
 Richard Barnes (rbarnes@umn.edu)
Based on Metz 2011.
 Pre
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM. Note that the NoData value is assumed to be a negative number less than any actual data value.
 Post
 flowdirs contains a D8 flow direction of each cell or a value NO_FLOW for those cells which are not part of the DEM.
 flowdirs has no cells which are not part of a continuous flow path leading to the edge of the DEM.
 Correctness:
 The correctness of this command is determined by inspection. (TODO)
 Parameters
&elevations
: A grid of cell elevations&flowdirs
: A grid of D8 flow directions
 template <class elev_t>

void
pit_mask
(const Array2D<elev_t> &elevations, Array2D<uint8_t> &pit_mask)¶ Indicates which cells are in pits.
This version of PriorityFlood starts on the edges of the DEM and then works its way inwards using a priority queue to determine the lowest cell which has a path to the edge. If a cell is lower than this cell then it is part of a pit and is given a value 1 to indicate this. The result is a grid where every cell which is in a pit is labeled.
 Author
 Richard Barnes (rbarnes@umn.edu)
 Pre
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM. Note that the NoData value is assumed to be a negative number less than any actual data value.
 Post
 pit_mask contains a 1 for each cell which is in a pit and a 0 for each cell which is not. The value 3 indicates NoData
 Correctness:
 The correctness of this command is determined by inspection. (TODO)
 Parameters
&elevations
: A grid of cell elevations&pit_mask
: A grid of indicating which cells are in pits
 template <class elev_t>

void
PriorityFloodWatersheds_Barnes2014
(Array2D<elev_t> &elevations, Array2D<int32_t> &labels, bool alter_elevations)¶ Gives a common label to all cells which drain to a common point.
All the edge cells of the DEM are given unique labels. This version of PriorityFlood starts on the edges of the DEM and then works its way inwards using a priority queue to determine the lowest cell which has a path to the edge. The neighbours of this cell are then given its label. All depressions are implicitly filled and digital dams removed. The result is a grid of cells where all cells with a common label drain to a common point.
 Author
 Richard Barnes (rbarnes@umn.edu)
 Pre
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM. Note that the NoData value is assumed to be a negative number less than any actual data value.
 Post
 elevations contains no depressions or digital dams, if alter_elevations** was set.
 labels contains a label for each cell indicating its membership in a given watershed. Cells bearing common labels drain to common points.
 Correctness:
 The correctness of this command is determined by inspection. (TODO)
 Parameters
elevations
: A grid of cell elevationslabels
: A grid to hold the watershed labelsalter_elevations
: If true, then elevations is altered as though PriorityFlood_Barnes2014() had been applied. The result is that all cells drain to the edges of the DEM. Otherwise, elevations is not altered.
 template <class elev_t>

void
PriorityFlood_Barnes2014_max_dep
(Array2D<elev_t> &elevations, uint64_t max_dep_size)¶ Fill depressions, but only if they’re small.
PriorityFlood starts on the edges of the DEM and then works its way inwards using a priority queue to determine the lowest cell which has a path to the edge. The neighbours of this cell are added to the priority queue if they are higher. If they are lower, they are raised to the elevation of the cell adding them, thereby filling in pits. The neighbors are then added to a “pit” queue which is used to flood pits. Cells which are higher than a pit being filled are added to the priority queue. In this way, pits are filled without incurring the expense of the priority queue.
 Author
 Richard Barnes (rbarnes@umn.edu)
When a depression is encountered this command measures its size before filling it. Only small depressions are filled.
 Pre
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM. Note that the NoData value is assumed to be a negative number less than any actual data value.
 Post
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM.
 elevations all landscape depressions <=max_dep_size are filled.
 Correctness:
 The correctness of this command is determined by inspection. (TODO)
 Parameters
&elevations
: A grid of cell elevationsmax_dep_size
: Depression must have <=max_dep_size cells to be filled
 template <class T>

void
Lindsay2016
(Array2D<T> &dem, int mode, bool fill_depressions, uint32_t maxpathlen, T maxdepth)¶
 template <>

void
Lindsay2016
(Array2D<uint8_t> &dem, int mode, bool fill_depressions, uint32_t maxpathlen, uint8_t maxdepth)¶
 template <>

void
Lindsay2016
(Array2D<int16_t> &dem, int mode, bool fill_depressions, uint32_t maxpathlen, int16_t maxdepth)¶
 template <>

void
Lindsay2016
(Array2D<uint16_t> &dem, int mode, bool fill_depressions, uint32_t maxpathlen, uint16_t maxdepth)¶
 template <class elev_t>

void
ProcessTraceQue_onepass
(Array2D<elev_t> &dem, Array2D<label_t> &labels, std::queue<int> &traceQueue, std::priority_queue<std::pair<elev_t, int>, std::vector<std::pair<elev_t, int>>, std::greater<std::pair<elev_t, int>>> &priorityQueue)¶
 template <class elev_t>

void
ProcessPit_onepass
(elev_t c_elev, Array2D<elev_t> &dem, Array2D<label_t> &labels, std::queue<int> &depressionQue, std::queue<int> &traceQueue, std::priority_queue<std::pair<elev_t, int>, std::vector<std::pair<elev_t, int>>, std::greater<std::pair<elev_t, int>>> &priorityQueue)¶
 template <class elev_t>

void
PriorityFlood_Zhou2016
(Array2D<elev_t> &dem)¶ Fills all pits and removes all digital dams from a DEM, quickly.
Works similarly to the PriorityFlood described by Barnes et al. (2014), but reduces the number of items which must pass through the priority queue, thus achieving greater efficiencies.
 Author
 G. Zhou, Z. Sun, S. Fu, Richard Barnes (this implementation)
 Pre
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM. Note that the NoData value is assumed to be a negative number less than any actual data value.
 Post
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM.
 elevations contains no landscape depressions or digital dams.
 Parameters
&dem
: A grid of cell elevations

static int
d8_masked_FlowDir
(const Array2D<int32_t> &flat_mask, const Array2D<int32_t> &labels, const int x, const int y)¶ Helper function to d8_flow_flats()
This determines a cell’s flow direction, taking into account flat membership. It is a helper function to d8_flow_flats()
 Author
 Richard Barnes (rbarnes@umn.edu)
 Return
 The flow direction of the cell
 Parameters
&flat_mask
: A mask from resolve_flats_barnes()&labels
: The labels output from resolve_flats_barnes()x
: x coordinate of celly
: y coordinate of cell
 template <class U>

void
d8_flow_flats
(const Array2D<int32_t> &flat_mask, const Array2D<int32_t> &labels, Array2D<U> &flowdirs)¶ Calculates flow directions in flats.
This determines flow directions within flats which have been resolved using resolve_flats_barnes().
 Author
 Richard Barnes (rbarnes@umn.edu)
Uses the helper function d8_masked_FlowDir()
 Pre
 flat_mask contains the number of increments to be applied to each cell to form a gradient which will drain the flat it is a part of.
 Any cell without a local gradient has a value of NO_FLOW in flowdirs**; all other cells have defined flow directions.
 If a cell is part of a flat, it has a value greater than zero in labels** indicating which flat it is a member of; otherwise, it has a value of 0.
 Post
 Every cell whose flow direction could be resolved by this algorithm (all drainable flats) will have a defined flow direction in flowdirs**. Any cells which could not be resolved (nondrainable flats) will still be marked NO_FLOW.
 Parameters
&flat_mask
: A mask from resolve_flats_barnes()&labels
: The labels output from resolve_flats_barnes()&flowdirs
: Returns flatresolved flow directions
 template <class U>

static void
BuildAwayGradient
(const Array2D<U> &flowdirs, Array2D<int32_t> &flat_mask, std::deque<GridCell> edges, std::vector<int> &flat_height, const Array2D<int32_t> &labels)¶ Build a gradient away from the high edges of the flats.
The queue of highedge cells developed in find_flat_edges() is copied into the procedure. A breadthfirst expansion labels cells by their distance away from terrain of higher elevation. The maximal distance encountered is noted.
 Author
 Richard Barnes (rbarnes@umn.edu)
 Pre
 Every cell in labels is marked either 0, indicating that the cell is not part of a flat, or a number greater than zero which identifies the flat to which the cell belongs.
 Any cell without a local gradient is marked NO_FLOW in flowdirs.
 Every cell in flat_mask is initialized to 0.
 edges contains, in no particular order, all the high edge cells of the DEM (those flat cells adjacent to higher terrain) which are part of drainable flats.
 Post
 flat_height will have an entry for each label value of labels indicating the maximal number of increments to be applied to the flat identified by that label.
 flat_mask will contain the number of increments to be applied to each cell to form a gradient away from higher terrain; cells not in a flat will have a value of 0.
 Parameters
&flowdirs
: A 2D array indicating each cell’s flow direction&flat_mask
: A 2D array for storing flat_mask&edges
: The highedge FIFO queue from find_flat_edges()&flat_height
: Vector with length equal to max number of labels&labels
: 2D array storing labels developed in label_this()
 template <class U>

static void
BuildTowardsCombinedGradient
(const Array2D<U> &flowdirs, Array2D<int32_t> &flat_mask, std::deque<GridCell> edges, std::vector<int> &flat_height, const Array2D<int32_t> &labels)¶ Builds gradient away from the low edges of flats, combines gradients.
The queue of lowedge cells developed in find_flat_edges() is copied into the procedure. A breadthfirst expansion labels cells by their distance away from terrain of lower elevation. This is combined with the gradient from BuildAwayGradient() to give the final increments of each cell in forming the flat mask.
 Author
 Richard Barnes (rbarnes@umn.edu)
 Pre
 Every cell in labels is marked either 0, indicating that the cell is not part of a flat, or a number greater than zero which identifies the flat to which the cell belongs.
 Any cell without a local gradient is marked NO_FLOW in flowdirs.
 Every cell in flat_mask has either a value of 0, indicating that the cell is not part of a flat, or a value greater than zero indicating the number of increments which must be added to it to form a gradient away from higher terrain.
 flat_height has an entry for each label value of labels indicating the maximal number of increments to be applied to the flat identified by that label in order to form the gradient away from higher terrain.
 edges contains, in no particular order, all the low edge cells of the DEM (those flat cells adjacent to lower terrain).
 Post
 flat_mask will contain the number of increments to be applied to each cell to form a superposition of the gradient away from higher terrain with the gradient towards lower terrain; cells not in a flat have a value of 0.
 Parameters
&flowdirs
: A 2D array indicating each cell’s flow direction&flat_mask
: A 2D array for storing flat_mask&edges
: The lowedge FIFO queue from find_flat_edges()&flat_height
: Vector with length equal to max number of labels&labels
: 2D array storing labels developed in label_this()
 template <class T>

static void
label_this
(int x0, int y0, const int label, Array2D<int32_t> &labels, const Array2D<T> &elevations)¶ Labels all the cells of a flat with a common label.
Performs a flood fill operation which labels all the cells of a flat with a common label. Each flat will have a unique label
 Author
 Richard Barnes (rbarnes@umn.edu)
 Pre
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM.
 labels has the same dimensions as elevations.
 **(x0,y0)** belongs to the flat which is to be labeled.
 label is a unique label which has not been previously applied to a flat.
 labels is initialized to zero prior to the first call to this function.
 Post
 **(x0,y0)** and every cell reachable from it by passing over only cells of the same elevation as it (all the cells in the flat to which c belongs) will be marked as label in labels.
 Parameters
x0
: xcoordinate of flood fill seedy0
: ycoordinate of floodfill seedlabel
: Label to apply to the cells&labels
: 2D array which will contain the labels&elevations
: 2D array of cell elevations
 template <class T, class U>

static void
find_flat_edges
(std::deque<GridCell> &low_edges, std::deque<GridCell> &high_edges, const Array2D<U> &flowdirs, const Array2D<T> &elevations)¶ Identifies cells adjacent to higher and lower terrain.
Cells adjacent to lower and higher terrain are identified and added to the appropriate queue
 Author
 Richard Barnes (rbarnes@umn.edu)
 Pre
 elevations contains the elevations of every cell or a value NoData for cells not part of the DEM.
 Any cell without a local gradient is marked NO_FLOW in flowdirs.
 Post
 high_edges will contain, in no particular order, all the high edge cells of the DEM: those flat cells adjacent to higher terrain.
 low_edges will contain, in no particular order, all the low edge cells of the DEM: those flat cells adjacent to lower terrain.
 Parameters
&low_edges
: Queue for storing cells adjacent to lower terrain&high_edges
: Queue for storing cells adjacent to higher terrain&flowdirs
: 2D array indicating flow direction for each cell&elevations
: 2D array of cell elevations
 template <class T, class U>

void
resolve_flats_barnes
(const Array2D<T> &elevations, const Array2D<U> &flowdirs, Array2D<int32_t> &flat_mask, Array2D<int32_t> &labels)¶ Performs the flat resolution by Barnes, Lehman, and Mulla.
TODO
 Author
 Richard Barnes (rbarnes@umn.edu)
 Pre
 elevations contains the elevations of every cell or the NoData value for cells not part of the DEM.
 Any cell without a local gradient is marked NO_FLOW in flowdirs.
 Post
 flat_mask will have a value greater than or equal to zero for every cell, indicating its number of increments. These can be used be used in conjunction with labels to determine flow directions without altering the DEM, or to alter the DEM in subtle ways to direct flow.
 labels will have values greater than or equal to 1 for every cell which is in a flat. Each flat’s cells will bear a label unique to that flat.
 Parameters
&elevations
: 2D array of cell elevations&flowdirs
: 2D array indicating flow direction of each cell&flat_mask
: 2D array which will hold incremental elevation mask&labels
: 2D array indicating flat membership
 template <class U>

void
d8_flats_alter_dem
(const Array2D<int32_t> &flat_mask, const Array2D<int32_t> &labels, Array2D<U> &elevations)¶ Alters the elevations of the DEM so that all flats drain.
This alters elevations within the DEM so that flats which have been resolved using resolve_flats_barnes() will drain.
 Author
 Richard Barnes (rbarnes@umn.edu)
 Pre
 flat_mask contains the number of increments to be applied to each cell to form a gradient which will drain the flat it is a part of.
 If a cell is part of a flat, it has a value greater than zero in labels** indicating which flat it is a member of; otherwise, it has a value of 0.
 Post
 Every cell whose part of a flat which could be drained will have its elevation altered in such a way as to guarantee that its flat does drain.
 Parameters
&flat_mask
: A mask from resolve_flats_barnes()&labels
: A grouping from resolve_flats_barnes()&elevations
: 2D array of elevations
 template <class T, class U>

void
barnes_flat_resolution_d8
(Array2D<T> &elevations, Array2D<U> &flowdirs, bool alter)¶

static float
dinf_masked_FlowDir
(const Array2D<int32_t> &flat_resolution_mask, const Array2D<int32_t> &groups, const int x, const int y)¶

void
dinf_flow_flats
(const Array2D<int32_t> &flat_resolution_mask, const Array2D<int32_t> &groups, Array2D<float> &flowdirs)¶
 template <class T>

void
resolve_flats_barnes_dinf
(const Array2D<T> &elevations, Array2D<float> &flowdirs)¶
 template <class T>

void
GradientTowardsLower
(const Array2D<T> &elevations, const Array2D<uint8_t> &flowdirs, flat_type &flats, Array2D<int32_t> &inc1)¶
 template <class T>

void
GradientAwayFromHigher
(const Array2D<T> &elevations, const Array2D<uint8_t> &flowdirs, flat_type &flats, Array2D<int32_t> &inc2)¶
 template <class T>

void
CombineGradients
(Array2D<T> &elevations, const Array2D<int32_t> &inc1, const Array2D<int32_t> &inc2, float epsilon)¶
 template <class T>

static int
d8_FlowDir
(const Array2D<T> &elevations, const int x, const int y)¶ Calculates the D8 flow direction of a cell.
This calculates the D8 flow direction of a cell using the D8 neighbour system, as defined in utility.h. Cells on the edge of the grid flow off the nearest edge.
 Author
 Richard Barnes (rbarnes@umn.edu)
Helper function for d8_flow_directions().
 Return
 The D8 flow direction of the cell
 Parameters
&elevations
: A DEMx
: x coordinate of celly
: y coordinate of cell
 template <class T, class U>

void
d8_flow_directions
(const Array2D<T> &elevations, Array2D<U> &flowdirs)¶ Calculates the D8 flow directions of a DEM.
This calculates the D8 flow directions of a DEM. Its argument ‘flowdirs’ will return a grid with flow directions using the D8 neighbour system, as defined in utility.h. The choice of data type for array2d must be able to hold exact values for all neighbour identifiers (usually [1,7]).
 Author
 Richard Barnes (rbarnes@umn.edu)
Uses d8_FlowDir() as a helper function.
 Parameters
&elevations
: A DEM&flowdirs
: Returns the flow direction of each cell
 template <class T>

static float
dinf_FlowDir
(const Array2D<T> &elevations, const int x, const int y)¶ Determine the Dinfinite flow direction of a cell.
This function determines the Dinfinite flow direction of a cell, as described by Tarboton (1997) and Barnes (2013, TODO). TODO
 Author
 Implementation by Richard Barnes (rbarnes@umn.edu)
 Return
 A floatingpoint value between [0,2*Pi) indicating flow direction
 Parameters
elevations
: A 2D grid of elevation datax
: xcoordinate of cell to determine flow direction fory
: ycoordinate of cell to determine flow direction for
 template <class T>

void
dinf_flow_directions
(const Array2D<T> &elevations, Array2D<float> &flowdirs)¶ Determine the Dinfinite flow direction of every cell in a grid.
This function runs dinf_FlowDir() on every cell in a grid which has a data value.
 Author
 Richard Barnes (rbarnes@umn.edu)
 Parameters
&elevations
: A 2D grid of elevation data&flowdirs
: A 2D grid which will contain the flow directions
 template <class E>

std::vector<float>
FM_Freeman
(const Array2D<E> &elevations, const double xparam)¶
 template <class E>

std::vector<float>
FM_Holmgren
(const Array2D<E> &elevations, const double xparam)¶
 template <class T>

static T
sgn
(T val)¶ Returns the sign (+1, 1, 0) of a number. Branchless.
 Author
 Richard Barnes (rbarnes@umn.edu)
 Return
 1 for a negative input, +1 for a positive input, and 0 for a zero input
 Parameters
val
: Input value
 template <class T, class U>

void
d8_flow_accum
(const Array2D<T> &flowdirs, Array2D<U> &area)¶ Calculates the D8 flow accumulation, given the D8 flow directions.
This calculates the D8 flow accumulation of a grid of D8 flow directions by calculating each cell’s dependency on its neighbours and then using a priorityqueue to process cells in a topofthewatersheddown fashion
 Author
 Richard Barnes (rbarnes@umn.edu)
 Parameters
&flowdirs
: A D8 flowdir grid from d8_flow_directions()&area
: Returns the upslope area of each cell
 template <class T, class U>

void
d8_upslope_cells
(int x0, int y0, int x1, int y1, const Array2D<T> &flowdirs, Array2D<U> &upslope_cells)¶ Calculates which cells ultimately D8flow through a given cell.
Given the coordinates x0,y0 of a cell and x1,y1 of another, possibly distinct, cell this draws a line between the two using the Bresenham LineDrawing Algorithm and returns a grid showing all the cells whose flow ultimately passes through the indicated cells.
 Author
 Richard Barnes (rbarnes@umn.edu)
The grid has the values:
1=Upslope cell 2=Member of initializing line All other cells have a noData() value
 Parameters
x0
: xcoordinate of start of liney0
: ycoordinate of start of linex1
: xcoordinate of end of liney1
: ycoordinate of end of line&flowdirs
: A D8 flowdir grid from d8_flow_directions()&upslope_cells
: A grid of 1/2/NoData, as in the description

static void
where_do_i_flow
(float flowdir, int &nhigh, int &nlow)¶

static void
area_proportion
(float flowdir, int nhigh, int nlow, float &phigh, float &plow)¶
 template <class T, class U>

void
dinf_upslope_area
(const Array2D<T> &flowdirs, Array2D<U> &area)¶ Calculate each cell’s Dinfinity flow accumulation value.
TODO
 Author
 Tarboton (1997), Richard Barnes (rbarnes@umn.edu)
 Parameters
flowdirs
: A grid of Dinfinite flow directions&area
: A grid of flow accumulation values
 template <class elev_t, class accum_t>

void
FA_Tarboton
(const Array2D<elev_t> &elevations, Array2D<accum_t> &accum)¶
 template <class elev_t, class accum_t>

void
FA_Holmgren
(const Array2D<elev_t> &elevations, Array2D<accum_t> &accum, double xparam)¶
 template <class elev_t, class accum_t>

void
FA_Quinn
(const Array2D<elev_t> &elevations, Array2D<accum_t> &accum)¶
 template <class elev_t, class accum_t>

void
FA_Freeman
(const Array2D<elev_t> &elevations, Array2D<accum_t> &accum, double xparam)¶
 template <class elev_t, class accum_t>

void
FA_FairfieldLeymarie
(const Array2D<elev_t> &elevations, Array2D<accum_t> &accum)¶
 template <class elev_t, class accum_t>

void
FA_Rho8
(const Array2D<elev_t> &elevations, Array2D<accum_t> &accum)¶
 template <class elev_t, class accum_t>

void
FA_D8
(const Array2D<elev_t> &elevations, Array2D<accum_t> &accum)¶
 template <class elev_t, class accum_t>

void
FA_OCallaghan
(const Array2D<elev_t> &elevations, Array2D<accum_t> &accum)¶
 template <class F, class E, class A, typename… Args>

void
FlowAccumulation
(F func, const Array2D<E> &elevations, Array2D<A> &accum, Args... args)¶
 template <class T, class U, class V>

void
TA_SPI
(const Array2D<T> &flow_accumulation, const Array2D<U> &riserun_slope, Array2D<V> &result)¶ Calculates the SPI terrain attribute.
\((\textit{CellSize}\cdot\textit{FlowAccumulation}+0.001)\cdot(\frac{1}{100}\textit{PercentSlope}+0.001)\)
 Author
 Richard Barnes (rbarnes@umn.edu)
 Pre
flow_accumulation
andpercent_slope
must be the same size Post
result
takes the properties and dimensions offlow_accumulation
 Parameters
&flow_accumulation
: A flow accumulation grid (dinf_upslope_area())&riserun_slope
: A percent_slope grid (d8_slope())&result
: Altered to return the calculated SPI
 template <class T, class U, class V>

void
TA_CTI
(const Array2D<T> &flow_accumulation, const Array2D<U> &riserun_slope, Array2D<V> &result)¶ Calculates the CTI terrain attribute.
\(\log{\frac{\textit{CellSize}\cdot\textit{FlowAccumulation}+0.001}{\frac{1}{100}\textit{PercentSlope}+0.001}}\)
 Author
 Richard Barnes (rbarnes@umn.edu)
 Pre
flow_accumulation
andpercent_slope
must be the same size Post
result
takes the properties and dimensions offlow_accumulation
 Parameters
&flow_accumulation
: A flow accumulation grid (dinf_upslope_area())&riserun_slope
: A percent_slope grid (d8_slope())&result
: Altered to return the calculated SPI
 template <class T>

static TA_Setup_Vars
TerrainSetup
(const Array2D<T> &elevations, const int x, const int y, const float zscale)¶
 template <class T>

static TA_Setup_Curves_Vars
TerrainCurvatureSetup
(const Array2D<T> &elevations, const int x, const int y, const float zscale)¶
 template <class T>

static double
Terrain_Aspect
(const Array2D<T> &elevations, const int x, const int y, const float zscale)¶ Calculates aspect in degrees in the manner of Horn 1981.
 Return
 Aspect in degrees in the manner of Horn 1981
 template <class T>

static double
Terrain_Slope_RiseRun
(const Array2D<T> &elevations, const int x, const int y, const float zscale)¶ Calculates the rise/run slope along the maximum gradient on a fitted surface over a 3x3 be neighbourhood in the manner of Horn 1981.
 Return
 Rise/run slope
 template <class T>

static double
Terrain_Curvature
(const Array2D<T> &elevations, const int x, const int y, const float zscale)¶
 template <class T>

static double
Terrain_Planform_Curvature
(const Array2D<T> &elevations, const int x, const int y, const float zscale)¶
 template <class T>

static double
Terrain_Profile_Curvature
(const Array2D<T> &elevations, const int x, const int y, const float zscale)¶
 template <class T>

static double
Terrain_Slope_Percent
(const Array2D<T> &elevations, const int x, const int y, const float zscale)¶
 template <class T>

static double
Terrain_Slope_Radian
(const Array2D<T> &elevations, const int x, const int y, const float zscale)¶
 template <class T>

static double
Terrain_Slope_Degree
(const Array2D<T> &elevations, const int x, const int y, const float zscale)¶
 template <class F, class T>

static void
TerrainProcessor
(F func, const Array2D<T> &elevations, const float zscale, Array2D<float> &output)¶ Calculate a variety of terrain attributes.
This calculates a variety of terrain attributes according to the work of Burrough 1998’s “Principles of Geographical Information Systems” (p. 190). This function calls d8_terrain_attrib_helper to calculate the actual attributes. This function may perform some postprocessing (such as on slope), but it’s purpose is essentially to just scan the grid and pass off the work to d8_terrain_attrib_helper().
 Author
 Richard Barnes (rbarnes@umn.edu), Burrough (1998)
Possible attribute values are
 TATTRIB_CURVATURE
 TATTRIB_PLANFORM_CURVATURE
 TATTRIB_PROFILE_CURVATURE
 TATTRIB_ASPECT
 TATTRIB_SLOPE_RISERUN
 TATTRIB_SLOPE_PERCENT
 TATTRIB_SLOPE_RADIAN
 TATTRIB_SLOPE_DEGREE
 Post
output
takes the properties and dimensions ofelevations
 Parameters
func
: The attribute function to be used&elevations
: An elevation gridzscale
: Value by which to scale elevation&output
: A grid to hold the results
 template <class T>

void
TA_slope_riserun
(const Array2D<T> &elevations, Array2D<float> &slopes, float zscale = 1.0f)¶ Calculates the slope as rise/run.
Calculates the slope using Horn 1981, as per Burrough 1998’s “Principles of Geographical Information Systems” (p. 190)
 Author
 Richard Barnes (rbarnes@umn.edu), Horn (1981)
 Parameters
&elevations
: An elevation grid&slopes
: A slope gridzscale
: DEM is scaled by this factor prior to calculation
 template <class T>

void
TA_slope_percentage
(const Array2D<T> &elevations, Array2D<float> &slopes, float zscale = 1.0f)¶ Calculates the slope as percentage.
Calculates the slope using Horn 1981, as per Burrough 1998’s “Principles of Geographical Information Systems” (p. 190)
 Author
 Richard Barnes (rbarnes@umn.edu), Horn (1981)
 Parameters
&elevations
: An elevation grid&slopes
: A slope gridzscale
: DEM is scaled by this factor prior to calculation
 template <class T>

void
TA_slope_degrees
(const Array2D<T> &elevations, Array2D<float> &slopes, float zscale = 1.0f)¶ Calculates the slope as degrees.
Calculates the slope using Horn 1981, as per Burrough 1998’s “Principles of Geographical Information Systems” (p. 190)
 Author
 Richard Barnes (rbarnes@umn.edu), Horn (1981)
 Parameters
&elevations
: An elevation grid&slopes
: A slope gridzscale
: DEM is scaled by this factor prior to calculation
 template <class T>

void
TA_slope_radians
(const Array2D<T> &elevations, Array2D<float> &slopes, float zscale = 1.0f)¶ Calculates the slope as radians.
Calculates the slope using Horn 1981, as per Burrough 1998’s “Principles of Geographical Information Systems” (p. 190)
 Author
 Richard Barnes (rbarnes@umn.edu), Horn (1981)
 Parameters
&elevations
: An elevation grid&slopes
: A slope gridzscale
: DEM is scaled by this factor prior to calculation
 template <class T>

void
TA_aspect
(const Array2D<T> &elevations, Array2D<float> &aspects, float zscale = 1.0f)¶ Calculates the terrain aspect.
Calculates the aspect per Horn 1981, as described by Burrough 1998’s “Principles of Geographical Information Systems” (p. 190)
 Author
 Richard Barnes (rbarnes@umn.edu), Horn (1981)
 Parameters
&elevations
: An elevation grid&aspects
: An aspect gridzscale
: DEM is scaled by this factor prior to calculation
 template <class T>

void
TA_curvature
(const Array2D<T> &elevations, Array2D<float> &curvatures, float zscale = 1.0f)¶ Calculates the terrain curvature per Zevenbergen and Thorne 1987.
Calculates the curvature per Zevenbergen and Thorne 1987, as described by Burrough 1998’s “Principles of Geographical Information Systems” (p. 190)
 Author
 Richard Barnes (rbarnes@umn.edu), Horn (1981)
 Parameters
&elevations
: An elevation grid&curvatures
: A curvature gridzscale
: DEM is scaled by this factor prior to calculation
 template <class T>

void
TA_planform_curvature
(const Array2D<T> &elevations, Array2D<float> &planform_curvatures, float zscale = 1.0f)¶ Calculates the terrain planform curvature per Zevenbergen and Thorne 1987.
Calculates the curvature per Zevenbergen and Thorne 1987, as described by Burrough 1998’s “Principles of Geographical Information Systems” (p. 190)
 Author
 Richard Barnes (rbarnes@umn.edu), Horn (1981)
 Parameters
&elevations
: An elevation grid&planform_curvatures
: A planform curvature gridzscale
: DEM is scaled by this factor prior to calculation
 template <class T>

void
TA_profile_curvature
(const Array2D<T> &elevations, Array2D<float> &profile_curvatures, float zscale = 1.0f)¶ Calculates the terrain profile curvature per Zevenbergen and Thorne 1987.
Calculates the curvature per Zevenbergen and Thorne 1987, as described by Burrough 1998’s “Principles of Geographical Information Systems” (p. 190)
 Author
 Richard Barnes (rbarnes@umn.edu), Horn (1981)
 Parameters
&elevations
: An elevation grid&profile_curvatures
: A profile curvature gridzscale
: DEM is scaled by this factor prior to calculation
Variables

const int
dx
[9] = {0, 1, 1, 0, 1, 1, 1, 0, 1}¶ x offsets of D8 neighbours, from a central cell

const int
dy
[9] = {0, 0, 1, 1, 1, 0, 1, 1, 1}¶ y offsets of D8 neighbours, from a central cell

const bool
n_diag
[9] = {0, 0, 1, 0, 1, 0, 1, 0, 1}¶ True along diagonal directions, false along north, south, east, west.

const int
d8_inverse
[9] = {0,5,6,7,8,1,2,3,4}¶ Directions from neighbours to the central cell. Neighbours are labeled 08. This is the inverse direction leading from a neighbour to the central cell.

const wchar_t
fdarrows
[9] = {L'·',L'←',L'↖',L'↑',L'↗',L'→',L'↘',L'↓',L'↙'}¶ Arrows indicating flow directions.

const double
SQRT2
= 1.414213562373095048801688724209698078569671875376948¶ sqrt(2), used to generate distances from a central cell to its neighbours

const double
dr
[9] = {0,1,SQRT2,1,SQRT2,1,SQRT2,1,SQRT2}¶ Distances from a central cell to each of its 8 neighbours.

const uint8_t
d8_arcgis
[9] = {0,16,32,64,128,1,2,4,8}¶ Convert from RichDEM flowdirs to ArcGIS flowdirs.

const uint8_t
FLOWDIR_NO_DATA
= 255¶ Used to indicate that a flowdir cell is NoData.

const d8_flowdir_t
NO_FLOW
= 0¶ Value used to indicate that a cell does not have a defined flow direction.

const int32_t
ACCUM_NO_DATA
= 1¶ Value used to indicate that a flow accumulation cell is NoData.

const uint8_t
GRID_LEFT
= 1¶ Indicates a tile is on the LHS of a DEM.

const uint8_t
GRID_TOP
= 2¶ Indicates a tile is on the top of a DEM.

const uint8_t
GRID_RIGHT
= 4¶ Indicates a tile is on the RHS of a DEM.

const uint8_t
GRID_BOTTOM
= 8¶ Indicates a tile is on the bottom of a DEM.

const std::string
git_hash
= "NO HASH SPECIFIED!"¶ Git hash of program’s source (used if RICHDEM_GIT_HASH is undefined)

const std::string
compilation_datetime
= __DATE__ " " __TIME__¶ Date and time of when the program was compiled (used if RICHDEM_COMPILE_TIME is undefined)
Richard Barnes.

const std::string
program_identifier
= program_name + " (hash=" + git_hash + ", compiled="+compilation_datetime + ")"¶ Richdem vX.X.X (hash=GIT HASH, compiled=COMPILATION DATE TIME)

const float
d8_to_dinf
[9] ={1, 4*M_PI/4, 3*M_PI/4, 2*M_PI/4, 1*M_PI/4, 0, 7*M_PI/4, 6*M_PI/4, 5*M_PI/4}¶

const int
dy_e1
[8] = { 0 , 1 , 1 , 0 , 0 , 1 , 1 , 0 }¶

const int
dx_e1
[8] = { 1 , 0 , 0 , 1 , 1 , 0 , 0 , 1 }¶

const int
dy_e2
[8] = {1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 }¶

const int
dx_e2
[8] = { 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 }¶

const double
ac
[8] = { 0., 1., 1., 2., 2., 3., 3., 4.}¶

const double
af
[8] = { 1., 1., 1., 1., 1., 1., 1., 1.}¶

const int
dinf_dx
[9] = {1, 1, 0, 1, 1, 1, 0, 1, 1}¶

const int
dinf_dy
[9] = {0, 1, 1, 1, 0, 1, 1, 1, 0}¶

typedef uint8_t

namespace
std
¶

file
Array2D.hpp
 #include “gdal.hpp”#include <vector>#include <iostream>#include <fstream>#include <iomanip>#include <cassert>#include <algorithm>#include <typeinfo>#include <stdexcept>#include <limits>#include <ctime>#include <unordered_set>#include <map>#include “richdem/common/logger.hpp”#include “richdem/common/version.hpp”#include “richdem/common/constants.hpp”#include “richdem/common/ManagedVector.hpp”
Defines a 2D array object with many convenient methods for working with raster data, along with several functions for checking file data types.
Richard Barnes (rbarnes@umn.edu), 2015

file
communicationthreads.hpp
 #include <cereal/types/string.hpp>#include <cereal/types/vector.hpp>#include <cereal/types/map.hpp>#include <cereal/archives/binary.hpp>#include <sstream>#include <vector>#include <queue>#include <map>#include <atomic>#include <iterator>#include <cassert>#include <iostream>#include <list>#include <chrono>
Defines

_unused
(x)¶
Functions
 template <class Fn>

void
CommInit
(int n, Fn &&fn, int *argc, char ***argv)¶
 template <class T, class U>

void
CommSend
(const T *a, const U *b, int dest, int tag)¶
 template <class T>

void
CommSend
(const T *a, nullptr_t, int dest, int tag)¶

int
CommGetTag
(int from)¶

int
CommGetSource
()¶

int
CommRank
()¶

int
CommSize
()¶

void
CommAbort
(int errorcode)¶
 template <class T, class U>

void
CommRecv
(T *a, U *b, int from)¶
 template <class T>

void
CommRecv
(T *a, nullptr_t, int from)¶
 template <class T>

void
CommBroadcast
(T *datum, int root)¶

void
CommFinalize
()¶

int
CommBytesSent
()¶

int
CommBytesRecv
()¶

void
CommBytesReset
()¶

void
CommBarrier
()¶


file
communication.hpp
 #include <mpi.h>#include <cereal/types/string.hpp>#include <cereal/types/vector.hpp>#include <cereal/types/map.hpp>#include <cereal/archives/binary.hpp>#include <sstream>#include <vector>#include <iterator>#include <cassert>#include <iostream>#include <thread>#include <chrono>
Abstract calls to MPI, allowing for transparent serialization and communication stats.
Richard Barnes (rbarnes@umn.edu), 2015
Defines

_unused
(x) Used to hide the fact that some variables are used only for assertions.
Typedefs

typedef uint64_t
comm_count_type
¶ Data type used for storing Tx/Rx byte counts.
Functions

void
CommInit
(int *argc, char ***argv)¶ Initiate communication (wrapper for MPI_Init)
 template <class T, class U>

msg_type
CommPrepare
(const T *a, const U *b)¶ Convert up to two objects into a combined serialized representation.
 template <class T>

msg_type
CommPrepare
(const T *a, std::nullptr_t)¶ Convert one object into a serialized representation.
 template <class T, class U>

void
CommSend
(const T *a, const U *b, int dest, int tag)¶ Serialize and send up to two objects.
 template <class T>

void
CommSend
(const T *a, std::nullptr_t, int dest, int tag)¶ Serialize and send a single object.

void
CommISend
(msg_type &msg, int dest, int tag)¶ Send a preserialized object using nonblocking communication.
The object must be preserialized because the buffer containing the serialization must persist until the communication is complete. It makes more sense to manage this buffer outside of this library.

int
CommGetTag
(int from)¶ Check tag of incoming message. Blocksing.

int
CommRank
()¶ Get my unique process identifier (i.e. rank)

int
CommSize
()¶ How many processes are active?

void
CommAbort
(int errorcode)¶ Abort; If any process calls this it will kill all the processes.
 template <class T, class U>

void
CommRecv
(T *a, U *b, int from)¶ Receive up to two objects and deserialize them.
 template <class T>

void
CommRecv
(T *a, std::nullptr_t, int from)¶ Receive one object and deserialize it.
 template <class T>

void
CommBroadcast
(T *datum, int root)¶ Broadcast a message to all of the processes. (TODO: An integer message?)

void
CommFinalize
()¶ Wrap things up politely; call this when all communication is done.

comm_count_type
CommBytesSent
()¶ Get the number of bytes sent by this process.
 Return
 Number of bytes sent by this process

comm_count_type
CommBytesRecv
()¶ Get the number of bytes received by this process.
 Return
 Number of bytes received by this process

void
CommBytesReset
()¶ Reset message size statistics to zero.
Variables

comm_count_type
bytes_sent
= 0¶ Number of bytes sent.

comm_count_type
bytes_recv
= 0¶ Number of bytes received.


file
constants.hpp
Defines a number of constants used by many of the algorithms.
RichDEM uses the following D8 neighbourhood. This is used by the dx[] and dy[] variables, among many others.
234 105 876
ArcGIS uses the following bits to indicate flow toward a given neighbour:
32 64 128 16 0 1 8 4 2

file
gdal.hpp

file
grid_cell.hpp
 #include <vector>#include <queue>#include <cmath>
Defines structures for addressing grid cells and associated queues.
Richard Barnes (rbarnes@umn.edu), 2015

file
Layoutfile.hpp
 #include “richdem/common/logger.hpp”#include <string>#include <vector>#include <fstream>#include <sstream>#include <iostream>#include <cassert>#include <stdexcept>
Defines classes used for reading and writing tiled datasets.
A layout file is a text file with the format:
file1.tif, file2.tif, file3.tif, file4.tif, file5.tif, file6.tif, file7.tif , file8.tif, ,
where each of fileX.tif is a tile of the larger DEM collectively described by all of the files. All of fileX.tif must have the same shape; the layout file specifies how fileX.tif are arranged in relation to each other in space. Blanks between commas indicate that there is no tile there: the algorithm will treat such gaps as places to route flow towards (as if they are oceans). Note that the files need not have TIF format: they can be of any type which GDAL can read. Paths to fileX.tif are taken to be relative to the layout file.
Richard Barnes (rbarnes@umn.edu), 2015

file
logger.hpp
 #include <iostream>#include <sstream>#include <string>

file
ManagedVector.hpp
 #include <memory>

file
memory.hpp
 #include <fstream>#include <string>
Defines functions for calculating memory usage.
Richard Barnes (rbarnes@umn.edu), 2015

file
ProgressBar.hpp
 #include <string>#include <iostream>#include <iomanip>#include <sys/time.h>#include <stdexcept>#include “richdem/common/timer.hpp”
Defines a handy progress bar object so users don’t get impatient.
The progress bar indicates to the user how much work has been completed, how much is left, and how long it is estimated to take. It accounts for multithreading by assuming uniform progress by all threads.
Define the global macro
NOPROGRESS
disables the progress bar, which may speed up the program.The progress bar looks like this:
[=================================== ] (70%  0.2s  1 threads)
Richard Barnes (rbarnes@umn.edu), 2015

file
random.cpp
 #include “random.hpp”#include <cassert>#include <random>#include <cstdint>#include <algorithm>#include <iostream>#include <functional>#include <limits>#include <sstream>

file
random.hpp
 #include <random>#include <string>

file
timer.hpp
 #include <sys/time.h>#include <stdexcept>
Defines the Timer class, which is used for timing code.
Richard Barnes (rbarnes@umn.edu), 2015

file
version.hpp
 #include <string>#include <iostream>
Defines RichDEM version, git hash, compilation time. Used for program/app headers and for processing history entries.
Richard Barnes (rbarnes@umn.edu), 2015

file
Barnes2014.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/common/grid_cell.hpp”#include “richdem/flowmet/d8_flowdirs.hpp”#include <queue>#include <limits>#include <iostream>#include <cstdlib>
Defines all the PriorityFlood algorithms described by Barnes (2014) “PriorityFlood: An Optimal DepressionFilling and WatershedLabeling Algorithm for Digital Elevation Models”.
Richard Barnes (rbarnes@umn.edu), 2015

file
depressions.hpp
 #include <richdem/depressions/Barnes2014.hpp>#include <richdem/depressions/Zhou2016.hpp>

file
Lindsay2016.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/common/grid_cell.hpp”#include “richdem/common/ProgressBar.hpp”#include “richdem/common/timer.hpp”#include <limits>

file
main.cpp
 #include “richdem/common/interface.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/depressions/Barnes2014.hpp”#include <string>#include <iostream>#include <cstdint>

file
main.cpp
 #include <iostream>#include <cstdint>#include <string>#include “richdem/common/Array2D.hpp”#include “richdem/flats/flat_resolution.hpp”#include “richdem/flats/garbrecht.hpp”

file
README.md

file
README.md

file
Zhou2016.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/common/timer.hpp”#include <queue>#include <vector>#include <map>#include <iostream>
Defines the PriorityFlood algorithm described by Zhou, G., Sun, Z., Fu, S., 2016. An efficient variant of the PriorityFlood algorithm for filling depressions in raster digital elevation models. Computers & Geosciences 90, Part A, 87 – 96. doi:http://dx.doi.org/10.1016/j.cageo.2016.02.021.
The code herein has been extensive modified by Richard Barnes (rbarnes@umn.edu) for inclusion with RichDEM.
Richard Barnes (rbarnes@umn.edu), 2015

file
flat_resolution.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/ProgressBar.hpp”#include “richdem/common/grid_cell.hpp”#include “richdem/flowmet/d8_flowdirs.hpp”#include <deque>#include <vector>#include <queue>#include <cmath>#include <limits>
Resolve flats according to Barnes (2014)
Contains code to generate an elevation mask which is guaranteed to drain a flat using a convergent flow pattern (unless it’s a mesa)
 Author
 Richard Barnes (rbarnes@umn.edu), 2012

file
flat_resolution_dinf.hpp
 #include “richdem/flats/flat_resolution.hpp”#include “richdem/flowmet/dinf_flowdirs.hpp”#include “richdem/common/logger.hpp”
Couples the Barnes (2014) flat resolution algorithm with the Tarboton (1997) Dinfinity flow metric.
 Author
 Richard Barnes

file
garbrecht.hpp
 #include <deque>#include <cstdint>#include <iostream>#include “richdem/common/Array2D.hpp”#include “richdem/common/grid_cell.hpp”#include “richdem/flowmet/d8_flowdirs.hpp”#include “richdem/common/logger.hpp”

file
generate_square_grid.cpp
 #include <iostream>#include <fstream>#include <cstdlib>

file
d8_flowdirs.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/common/ProgressBar.hpp”
Functions for calculating D8 flow directions.
Richard Barnes (rbarnes@umn.edu), 2015

file
dinf_flowdirs.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/common/ProgressBar.hpp”
Defines the Dinfinite flow routing method described by Tarboton (1997)
This file implements the Dinfinite flow routing method originally described by Tarboton (1997). It incorporates minor alterations and additional safeguards described in Barnes (TODO).
Richard Barnes (rbarnes@umn.edu), 2015
Defines

dinf_NO_DATA
¶ Value used to indicate that a flow direction cell has no data.


file
Fairfield1991.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/common/ProgressBar.hpp”#include “richdem/common/random.hpp”

file
Freeman1991.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/common/ProgressBar.hpp”

file
Holmgren1994.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/common/ProgressBar.hpp”

file
OCallaghan1984.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/common/ProgressBar.hpp”

file
Orlandini2003.hpp

file
Quinn1991.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/common/ProgressBar.hpp”

file
Seibert2007.hpp

file
Tarboton1997.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/common/ProgressBar.hpp”

file
d8_methods.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/common/constants.hpp”#include “richdem/common/grid_cell.hpp”#include “richdem/common/ProgressBar.hpp”#include <queue>#include <stdexcept>
Defines a number of functions for calculating terrain attributes.
Richard Barnes (rbarnes@umn.edu), 2015

file
dinf_methods.hpp
 #include <cmath>#include <queue>#include “richdem/common/logger.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/common/constants.hpp”#include “richdem/common/ProgressBar.hpp”#include “richdem/common/grid_cell.hpp”
Terrain attributes that can only be calculated with Tarboton’s Dinfinity flow metric.
This file implements the Dinfinite flow routing method originally described by Tarboton (1997). It incorporates minor alterations and additional safeguards described in Barnes (2013, TODO).
 Author
 Richard Barnes (rbarnes@umn.edu), 2015

file
flow_accumulation.hpp
 #include <richdem/flowmet/Fairfield1991.hpp>#include <richdem/flowmet/Freeman1991.hpp>#include <richdem/flowmet/Holmgren1994.hpp>#include <richdem/flowmet/OCallaghan1984.hpp>#include <richdem/flowmet/Orlandini2003.hpp>#include <richdem/flowmet/Quinn1991.hpp>#include <richdem/flowmet/Seibert2007.hpp>#include <richdem/flowmet/Tarboton1997.hpp>#include <richdem/methods/flow_accumulation_generic.hpp>

file
flow_accumulation_generic.hpp
 #include <richdem/common/Array2D.hpp>#include <richdem/common/logger.hpp>#include <richdem/common/ProgressBar.hpp>#include <queue>

file
strahler.hpp

file
terrain_attributes.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/common/constants.hpp”#include “richdem/common/ProgressBar.hpp”

file
misc_methods.hpp
 #include <cmath>#include <queue>#include <stdexcept>#include <cassert>#include “richdem/common/Array2D.hpp”#include “richdem/common/constants.hpp”#include “richdem/common/ProgressBar.hpp”
Terrain attributes that can only be calculated with Tarboton’s Dinfinity flow metric.
This file implements the Dinfinite flow routing method originally described by Tarboton (1997). It incorporates minor alterations and additional safeguards described in Barnes (2013, TODO).
 Author
 Richard Barnes (rbarnes@umn.edu), 2015
Enums

enum
PerimType
¶ Calculates the perimeter of a digital elevation model.
Calculates the perimeter of a DEM in one of several ways: CELL_COUNT  # of cells bordering edges or NoData cells SQUARE_EDGE  Summation of all edges touch borders or NoData cells
 Author
 Richard Barnes (rbarnes@umn.edu)
 Return
 The surface area of the digital elevation model
 Parameters
&arr
:
Values:

CELL_COUNT
¶ Counts # of cells bordering DEM edges or NoData cells.

SQUARE_EDGE
¶ Adds all cell edges bordering DEM edges or NoData cells.
Functions
 template <class T>

double
dem_surface_area
(const Array2D<T> &elevations, const double zscale)¶ Calculate the surface of a digital elevation model.
Calculates the surface area of a digital elevation model by connecting the central points of cells with triangles and then calculating the area of the portion of each triangle which falls within the focal cell. The method is described in detail in Jenness (2004) <doi:10.2193/00917648(2004)032[0829:CLSAFD]2.0.CO;2>
 Author
 Jenness (2004), Richard Barnes (rbarnes@umn.edu)
 Return
 The surface area of the digital elevation model
 Parameters
&elevations
: A grid of elevationszscale
: DEM is scaled by this factor prior to calculation

file
richdem.cpp
 #include “richdem/common/logger.hpp”#include <map>#include <string>

file
richdem.hpp
 #include “common/Array2D.hpp”#include “common/constants.hpp”#include “common/grid_cell.hpp”#include “common/ManagedVector.hpp”#include “common/memory.hpp”#include “common/ProgressBar.hpp”#include “common/random.hpp”#include “common/timer.hpp”#include “common/version.hpp”#include “depressions/Barnes2014.hpp”#include “depressions/depressions.hpp”#include “depressions/Lindsay2016.hpp”#include “depressions/Zhou2016.hpp”#include “flats/flat_resolution.hpp”#include “flats/flat_resolution_dinf.hpp”#include “flowmet/d8_flowdirs.hpp”#include “flowmet/dinf_flowdirs.hpp”#include “flowmet/Fairfield1991.hpp”#include “flowmet/Freeman1991.hpp”#include “flowmet/Holmgren1994.hpp”#include “flowmet/OCallaghan1984.hpp”#include “flowmet/Orlandini2003.hpp”#include “flowmet/Quinn1991.hpp”#include “flowmet/Seibert2007.hpp”#include “flowmet/Tarboton1997.hpp”#include “methods/d8_methods.hpp”#include “methods/dinf_methods.hpp”#include “methods/flow_accumulation.hpp”#include “methods/flow_accumulation_generic.hpp”#include “methods/strahler.hpp”#include “methods/terrain_attributes.hpp”

file
A2Array2D.hpp
 #include “richdem/common/logger.hpp”#include “richdem/common/Layoutfile.hpp”#include “richdem/common/Array2D.hpp”#include “richdem/tiled/lru.hpp”#include “gdal_priv.h”
Experimental tile manager for large datasets (TODO)
 Author
 Richard Barnes

file
lru.hpp
 #include <list>#include <unordered_map>
Defines a LeastRecently Used (LRU) cache class.
Richard Barnes (rbarnes@umn.edu), 2016

page
md__home_rick_projects_watershed_richdem_include_richdem_depressions_README
Title of Manuscript: PriorityFlood: An Optimal DepressionFilling and WatershedLabeling Algorithm for Digital Elevation Models
Authors: Richard Barnes, Clarence Lehman, David Mulla
Corresponding Author: Richard Barnes (rbarnes@umn.edu)
DOI Number of Manuscript 10.1016/j.cageo.2013.04.024
Code Repositories
This repository contains a reference implementation of the algorithms presented in the manuscript above. These implementations were used in performing the tests described in the manuscript.
There is source code for every pseudocode algorithm presented in the manuscript. All the code can be compiled simply by running make. The result is a program called priority_flood.exe.
This program reads in a DEM file specified on the command line. The file may be any ArcGrid ASCII file. The program will run one of the algorithms described in the manuscript (and below), store the result in an output file, and report how long this took.
The program is run by typing:
./priority_flood.exe <ALGORITHM NUMBER> <INPUT DEM> ./priority_flood.exe 3 inputdata.asc
The algorithms available are described briefly below and in greater detail in the manuscript.
 Algorithm 1: PriorityFlood This algorithm alters the input DEM to produce an output with no depressions or digital dams. Every cell which would have been in a depression is increased to the level of that depression’s outlet, leaving a flat region in its place. It runs slower than Algorithm 2, but is otherwise the same. The result is saved to outpforiginal.
 Algorithm 2: Improved PriorityFlood This algorithm alters the input DEM to produce an output with no depressions or digital dams. Every cell which would have been in a depression is increased to the level of that depression’s outlet, leaving a flat region in its place. It runs faster than Algorithm 1, but is otherwise the same. The result is saved to outpfimproved.
 Algorithm 3: PriorityFlood+Epsilon This algorithm alters the input DEM to produce an output with no depressions or digital dams. Every cell which would have been in a depression is increased to the level of that depression’s output, plus an additional increment which is sufficient to direct flow to the periphery of the DEM. The result is saved to outpfepsilon.
 Algorithm 4: PriorityFlood+FlowDirs This algorithm determines a D8 flow direction for every cell in the DEM by implicitly filling depressions and eliminating digital dams. Though all depressions are guaranteed to drain, local elevation information is still used to determine flow directions within a depression. It is, in essence, a depressioncarving algorithm. The result is saved to outpfflowdirs.
 Algorithm 5: PriorityFlood+Watershed Labels For each cell c in a DEM, this algorithm determines which cell on the DEM’s periphery c will drain to. c is then given a label which corresponds to the peripheral cell. All cells bearing a common label belong to the same watershed. The result is saved to outpfwlabels.
Algorithm 4: PriorityFlood+FlowDirs and its output, outpfflowdirs, use the D8 neighbour system to indicate flow directions. In this system all the flow from a central cell is directed to a single neighbour which is represented by a number according to the following system where 0 indicates the central cell.
234 105 876
The directory src/ contains the source code for the reference implementations. All the source code is drawn from the RichDEM hydroanalysis package. At the time of writing, the entire RichDEM code base could be downloaded from: https://github.com/rbarnes
Assumptions
All of the algorithms assume that cells marked as having NoData will have extremely negative numerical values: less than the value of any of the actual data. NaN is considered to be less than all values, including negative infinity.
Notes on the Manuscript
Work by Cris Luengo on the speed of various priority queue algorithms is discussed in the manuscript. His website providing code for his implementatations is here.
Updates
Commit 51f9a7838d3e88628ef6c74846edd0cb18e7ffe6 (020150925) introduced a number of changes to the code versus what was originally published with the manuscript. The old codebase uses ASCIIformatted data for input and output; the new codebase uses GDAL to handle many kinds of data.
The old codebase had the advantage of not relying on external libraries and being readily accessible to all parties. It had the disadvantage of being a slow, clumsy, and limited way to work with the data. As of 020150925, the code requires the use of the GDAL library greatly expanding the data formats and data types which can be worked with, as well as greatly speeding up I/O.
Note that using the aforementioned 51f9a7838d directly will result in silent casting of your data to the
float
type; commit 8b11f535af23368d3bd26609cc88df3dbb7111f1 (020150928) fixes this issue.Additionally, the library now uses C++ for all streaming operations except the progress bar.

page
md__home_rick_projects_watershed_richdem_include_richdem_flats_README
Title of Manuscript: An Efficient Assignment of Drainage Direction Over Flat Surfaces in Raster Digital Elevation Models
Authors: Richard Barnes, Clarence Lehman, David Mulla
Corresponding Author: Richard Barnes (rbarnes@umn.edu)
DOI Number of Manuscript 10.1016/j.cageo.2013.01.009
Code Repositories
This repository contains a reference implementation of the algorithms presented in the manuscript above. It also contains a reference implementation of the algorithm presented by Garbrecht and Martz (1997). These implementations were used in performing speed comparison tests in the manuscript.
All the programs can be produced simply by running make.
The program generate_square_grid.exe makes a square DEM with a single outlet near the bottomleft corner. The grid size is specified as a commandline argument.
The two reference implementations use the D8 neighbour system to indicate flow directions. In this system all the flow from a central cell is directed to a single neighbour which is represented by a number according to the following system where 0 indicates the central cell.
234 105 876
The program barnes_algorithm.exe reads in a DEM file specified on the command line. The file may be generated by generate_square_grid.exe, but may also be any ArcGrid ASCII file. The program will time itself and report the results back. The program will print the determined flow directions for the DEM to a file named out_barnes. The determined flow directions are also printed as a matrix of arrows to out_barnes_arrows.
The program garbrecht_algorithm.exe attempts to reproduce the algorithm described by Garbrecht and Martz (1997). It accepts an ArcGRID ASCII file as a command line input. The input file may also be generated with generate_square_grid.exe. Note that this implementation does not apply itself iteratively, meaning that some flats will be unresolvable. It writes the determined flow directions to out_garbrecht. The determined flow directions are also printed as a matrix of arrows to out_garbrecht_arrows.
The directory src/ contains the source code for reference implementations. The source for the improved algorithm is drawn from the RichDEM hydroanalysis package. All code can be compiled by running the makefile included in the root directory. Running the BASH script FIRST_RUN will compile everything and run the programs.
At the time of writing, the entire RichDEM code base could be downloaded from: https://github.com/rbarnes

page
todo

dir
/home/rick/projects/watershed/richdem/include/richdem/common

dir
/home/rick/projects/watershed/richdem/include/richdem/depressions

dir
/home/rick/projects/watershed/richdem/include/richdem/flats

dir
/home/rick/projects/watershed/richdem/include/richdem/flowmet

dir
/home/rick/projects/watershed/richdem/include

dir
/home/rick/projects/watershed/richdem/include/richdem/methods

dir
/home/rick/projects/watershed/richdem/include/richdem/misc

dir
/home/rick/projects/watershed/richdem/include/richdem

dir
/home/rick/projects/watershed/richdem/include/richdem/tiled