diff --git a/Makefile b/Makefile index c01317ff6..8b1398a92 100644 --- a/Makefile +++ b/Makefile @@ -2,7 +2,7 @@ DATA_DIR := data GIT_SAMPLE_DATA_REPO := https://bitbucket.org/natcap/invest-sample-data.git GIT_SAMPLE_DATA_REPO_PATH := $(DATA_DIR)/invest-sample-data -GIT_SAMPLE_DATA_REPO_REV := 0f8b41557753dad3670ba8220f41650b51435a93 +GIT_SAMPLE_DATA_REPO_REV := 08b1a727598fcdbdc9a787f0cf1d47dc6bac3aec GIT_TEST_DATA_REPO := https://bitbucket.org/natcap/invest-test-data.git GIT_TEST_DATA_REPO_PATH := $(DATA_DIR)/invest-test-data diff --git a/src/natcap/invest/managed_raster/ManagedRaster.h b/src/natcap/invest/managed_raster/ManagedRaster.h index 86955e99c..ac40c7ee6 100644 --- a/src/natcap/invest/managed_raster/ManagedRaster.h +++ b/src/natcap/invest/managed_raster/ManagedRaster.h @@ -1,21 +1,14 @@ -#include +#ifndef NATCAP_INVEST_MANAGEDRASTER_H_ +#define NATCAP_INVEST_MANAGEDRASTER_H_ #include "gdal.h" #include "gdal_priv.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include // For std::ptrdiff_t +#include #include "LRUCache.h" -int MANAGED_RASTER_N_BLOCKS = pow(2, 6); +int MANAGED_RASTER_N_BLOCKS = static_cast(pow(2, 6)); // given the pixel neighbor numbering system // 3 2 1 // 4 x 0 @@ -32,6 +25,9 @@ int INFLOW_OFFSETS[8] = {4, 5, 6, 7, 0, 1, 2, 3}; typedef std::pair BlockBufferPair; +class D8 {}; +class MFD {}; + class NeighborTuple { public: int direction, x, y; @@ -70,28 +66,27 @@ class ManagedRaster { int write_mode; int closed; double nodata; + double* geotransform; + int hasNodata; ManagedRaster() { } + // Creates new instance of ManagedRaster. Opens the raster with GDAL, + // stores important information about the dataset, and creates a cache + // that will be used to efficiently read blocks from the raster. + // Args: + // raster_path: path to raster that has block sizes that are + // powers of 2. If not, an exception is raised. + // band_id: which band in `raster_path` to index. Uses GDAL + // notation that starts at 1. + // write_mode: if true, this raster is writable and dirty + // memory blocks will be written back to the raster as blocks + // are swapped out of the cache or when the object deconstructs. ManagedRaster(char* raster_path, int band_id, bool write_mode) : raster_path { raster_path } , band_id { band_id } , write_mode { write_mode } { - // """Create new instance of Managed Raster. - - // Args: - // raster_path (char*): path to raster that has block sizes that are - // powers of 2. If not, an exception is raised. - // band_id (int): which band in `raster_path` to index. Uses GDAL - // notation that starts at 1. - // write_mode (boolean): if true, this raster is writable and dirty - // memory blocks will be written back to the raster as blocks - // are swapped out of the cache or when the object deconstructs. - - // Returns: - // None. - // """ GDALAllRegister(); dataset = (GDALDataset *) GDALOpen( raster_path, GA_Update ); @@ -110,7 +105,10 @@ class ManagedRaster { block_xmod = block_xsize - 1; block_ymod = block_ysize - 1; - nodata = band->GetNoDataValue(); + nodata = band->GetNoDataValue( &hasNodata ); + + geotransform = (double *) CPLMalloc(sizeof(double) * 6); + dataset->GetGeoTransform(geotransform); if (((block_xsize & (block_xsize - 1)) != 0) or ( (block_ysize & (block_ysize - 1)) != 0)) { @@ -119,8 +117,8 @@ class ManagedRaster { "This error is happening in the ManagedRaster.h extension."); } - block_xbits = log2(block_xsize); - block_ybits = log2(block_ysize); + block_xbits = static_cast(log2(block_xsize)); + block_ybits = static_cast(log2(block_ysize)); // integer floor division block_nx = (raster_x_size + block_xsize - 1) / block_xsize; @@ -141,8 +139,8 @@ class ManagedRaster { closed = 0; } + // Sets the pixel at `xi,yi` to `value` void inline set(long xi, long yi, double value) { - // Set the pixel at `xi,yi` to `value` int block_xi = xi / block_xsize; int block_yi = yi / block_ysize; @@ -163,8 +161,8 @@ class ManagedRaster { } } + // Returns the value of the pixel at `xi,yi`. double inline get(long xi, long yi) { - // Return the value of the pixel at `xi,yi` int block_xi = xi / block_xsize; int block_yi = yi / block_ysize; @@ -176,6 +174,8 @@ class ManagedRaster { } double* block = lru_cache->get(block_index); + + // Using the property n % 2^i = n & (2^i - 1) // to efficienty compute the modulo: yi % block_xsize int idx = ((yi & block_ymod) * actualBlockWidths[block_index]) + (xi & block_xmod); @@ -184,6 +184,9 @@ class ManagedRaster { return value; } + // Reads a block from the raster and saves it to the cache. + // Args: + // block_index: Index of the block to read, counted from the top-left void _load_block(int block_index) { int block_xi = block_index % block_nx; int block_yi = block_index / block_nx; @@ -259,16 +262,12 @@ class ManagedRaster { } } + // Closes the ManagedRaster and frees up resources. + // This call writes any dirty blocks to disk, frees up the memory + // allocated as part of the cache, and frees all GDAL references. + // Any subsequent calls to any other functions in _ManagedRaster will + // have undefined behavior. void close() { - // """Close the _ManagedRaster and free up resources. - - // This call writes any dirty blocks to disk, frees up the memory - // allocated as part of the cache, and frees all GDAL references. - - // Any subsequent calls to any other functions in _ManagedRaster will - // have undefined behavior. - // """ - if (closed) { return; } @@ -292,6 +291,7 @@ class ManagedRaster { // write the changed value back if desired CPLFree(it->second); } + GDALClose( (GDALDatasetH) dataset ); return; } @@ -337,29 +337,26 @@ class ManagedRaster { } }; +// Represents a flow direction raster, which may be of type MFD or D8 +template class ManagedFlowDirRaster: public ManagedRaster { - public: - ManagedFlowDirRaster() {} ManagedFlowDirRaster(char* raster_path, int band_id, bool write_mode) : ManagedRaster(raster_path, band_id, write_mode) {} + // Checks if a given pixel is a local high point. (MFD implementation) + // Args: + // xi: x coord in pixel space of the pixel to consider + // yi: y coord in pixel space of the pixel to consider + // Returns: + // true if the pixel is a local high point, i.e. it has no + // upslope neighbors; false otherwise. + template::value>* = nullptr> bool is_local_high_point(int xi, int yi) { - // """Check if a given pixel is a local high point. - - // Args: - // xi (int): x coord in pixel space of the pixel to consider - // yi (int): y coord in pixel space of the pixel to consider - - // Returns: - // True if the pixel is a local high point, i.e. it has no - // upslope neighbors; False otherwise. - // """ - int flow_dir_j; + int flow_dir_j, flow_ji; long xj, yj; - float flow_ji; for (int n_dir = 0; n_dir < 8; n_dir++) { xj = xi + COL_OFFSETS[n_dir]; @@ -367,7 +364,7 @@ class ManagedFlowDirRaster: public ManagedRaster { if (xj < 0 or xj >= raster_x_size or yj < 0 or yj >= raster_y_size) { continue; } - flow_dir_j = get(xj, yj); + flow_dir_j = static_cast(get(xj, yj)); flow_ji = (0xF & (flow_dir_j >> (4 * FLOW_DIR_REVERSE_DIRECTION[n_dir]))); if (flow_ji) { @@ -375,55 +372,83 @@ class ManagedFlowDirRaster: public ManagedRaster { } } return true; - } + // Checks if a given pixel is a local high point. (D8 implementation) + // Args: + // xi: x coord in pixel space of the pixel to consider + // yi: y coord in pixel space of the pixel to consider + // Returns: + // true if the pixel is a local high point, i.e. it has no + // upslope neighbors; false otherwise. + template::value>* = nullptr> + bool is_local_high_point(int xi, int yi) { + int flow_dir_j; + long xj, yj; + + for (int n_dir = 0; n_dir < 8; n_dir++) { + xj = xi + COL_OFFSETS[n_dir]; + yj = yi + ROW_OFFSETS[n_dir]; + if (xj < 0 or xj >= raster_x_size or yj < 0 or yj >= raster_y_size) { + continue; + } + flow_dir_j = static_cast(get(xj, yj)); + if (flow_dir_j == FLOW_DIR_REVERSE_DIRECTION[n_dir]) { + return false; + } + } + return true; + } }; -struct Pixel { - ManagedFlowDirRaster raster; +// Represents a pixel in a ManagedFlowDirectionRaster +template +class Pixel { +public: + ManagedFlowDirRaster raster; int x; int y; int val; Pixel() {} - Pixel(ManagedFlowDirRaster raster, int x, int y) : raster(raster), x(x), y(y) { + Pixel(ManagedFlowDirRaster raster, int x, int y) : raster(raster), x(x), y(y) { double v = raster.get(x, y); val = static_cast(v); } }; +// Returned by the `.end()` method of Neighbor iterable classes static inline NeighborTuple endVal = NeighborTuple(8, -1, -1, -1); - -struct NeighborIterator { +// Iterates over all eight neighboring pixels of a given pixel +// This and subsequent iterator classes were written with a lot of help from +// https://internalpointers.com/post/writing-custom-iterators-modern-cpp +template +class NeighborIterator { +public: using iterator_category = std::forward_iterator_tag; using difference_type = std::ptrdiff_t; using value_type = NeighborTuple; - using pointer = NeighborTuple*; // or also value_type* - using reference = NeighborTuple&; // or also value_type& + using pointer = NeighborTuple*; + using reference = NeighborTuple&; - Pixel pixel; + Pixel pixel; pointer m_ptr = nullptr; int i = 0; NeighborIterator() {} - NeighborIterator(NeighborTuple* n) { - m_ptr = n; - } - NeighborIterator(Pixel pixel) : pixel(pixel) { - next(); - } + NeighborIterator(NeighborTuple* n) { m_ptr = n; } + NeighborIterator(const Pixel pixel) : pixel(pixel) { next(); } reference operator*() const { return *m_ptr; } pointer operator->() { return m_ptr; } // Prefix increment - NeighborIterator& operator++() { this->next(); return *this; } + NeighborIterator& operator++() { next(); return *this; } // Postfix increment - NeighborIterator operator++(int) { NeighborIterator tmp = *this; ++(*this); return tmp; } + NeighborIterator operator++(int) { NeighborIterator tmp = *this; ++(*this); return tmp; } friend bool operator== (const NeighborIterator& a, const NeighborIterator& b) { return a.m_ptr == b.m_ptr; @@ -432,6 +457,7 @@ struct NeighborIterator { return a.m_ptr != b.m_ptr; }; + // Increments the pointer to the next neighbor virtual void next() { long xj, yj, flow; if (i == 8) { @@ -441,204 +467,343 @@ struct NeighborIterator { xj = pixel.x + COL_OFFSETS[i]; yj = pixel.y + ROW_OFFSETS[i]; flow = (pixel.val >> (i * 4)) & 0xF; - m_ptr = new NeighborTuple(i, xj, yj, flow); + m_ptr = new NeighborTuple(i, xj, yj, static_cast(flow)); i++; } }; -class DownslopeNeighborIterator: public NeighborIterator { +// Iterates over neighbor pixels that are downslope of a given pixel, +// in either MFD or D8 mode +template +class DownslopeNeighborIterator: public NeighborIterator { public: - DownslopeNeighborIterator(): NeighborIterator() {} - DownslopeNeighborIterator(NeighborTuple* n): NeighborIterator(n) {} - DownslopeNeighborIterator(Pixel p) { - pixel = p; + + DownslopeNeighborIterator(): NeighborIterator() {} + DownslopeNeighborIterator(NeighborTuple* n): NeighborIterator(n) {} + DownslopeNeighborIterator(const Pixel p) { + this->pixel = p; next(); } + DownslopeNeighborIterator& operator++() { next(); return *this; } + DownslopeNeighborIterator operator++(int) { DownslopeNeighborIterator tmp = *this; ++(*this); return tmp; } + // Increments the pointer to the next downslope neighbor (MFD) + template::value>* = nullptr> void next() { long xj, yj, flow; - delete m_ptr; - m_ptr = nullptr; - if (i == 8) { - m_ptr = &endVal; + delete this->m_ptr; + this->m_ptr = nullptr; + if (this->i == 8) { + this->m_ptr = &endVal; return; } - xj = pixel.x + COL_OFFSETS[i]; - yj = pixel.y + ROW_OFFSETS[i]; - if (xj < 0 or xj >= pixel.raster.raster_x_size or - yj < 0 or yj >= pixel.raster.raster_y_size) { - i++; + xj = this->pixel.x + COL_OFFSETS[this->i]; + yj = this->pixel.y + ROW_OFFSETS[this->i]; + if (xj < 0 or xj >= this->pixel.raster.raster_x_size or + yj < 0 or yj >= this->pixel.raster.raster_y_size) { + this->i++; next(); return; } - flow = (pixel.val >> (i * 4)) & 0xF; + flow = (this->pixel.val >> (this->i * 4)) & 0xF; if (flow) { - m_ptr = new NeighborTuple(i, xj, yj, flow); - i++; + this->m_ptr = new NeighborTuple(this->i, xj, yj, static_cast(flow)); + this->i++; return; } else { - i++; + this->i++; next(); } } + + // Increments the pointer to the next downslope neighbor (D8) + template::value>* = nullptr> + void next() { + long xj, yj; + delete this->m_ptr; + this->m_ptr = nullptr; + + if (this->i == 8) { + this->m_ptr = &endVal; + return; + } + xj = this->pixel.x + COL_OFFSETS[this->pixel.val]; + yj = this->pixel.y + ROW_OFFSETS[this->pixel.val]; + if (xj < 0 or xj >= this->pixel.raster.raster_x_size or + yj < 0 or yj >= this->pixel.raster.raster_y_size) { + this->m_ptr = &endVal; + return; + } + this->i = 8; + this->m_ptr = new NeighborTuple(this->pixel.val, xj, yj, 1); + return; + } }; -class DownslopeNeighborNoSkipIterator: public NeighborIterator { +// Iterates over neighbor pixels that are downslope of a given pixel, +// without skipping pixels that are out-of-bounds of the raster, +// in either MFD or D8 mode +template +class DownslopeNeighborNoSkipIterator: public NeighborIterator { public: - DownslopeNeighborNoSkipIterator(): NeighborIterator() {} - DownslopeNeighborNoSkipIterator(NeighborTuple* n): NeighborIterator(n) {} - DownslopeNeighborNoSkipIterator(Pixel p) { - pixel = p; + DownslopeNeighborNoSkipIterator(): NeighborIterator() {} + DownslopeNeighborNoSkipIterator(NeighborTuple* n): NeighborIterator(n) {} + DownslopeNeighborNoSkipIterator(const Pixel p) { + this->pixel = p; next(); } + DownslopeNeighborNoSkipIterator& operator++() { next(); return *this; } + DownslopeNeighborNoSkipIterator operator++(int) { DownslopeNeighborNoSkipIterator tmp = *this; ++(*this); return tmp; } + // Increments the pointer to the next downslope neighbor (MFD) + template::value>* = nullptr> void next() { long xj, yj, flow; - delete m_ptr; - m_ptr = nullptr; - if (i == 8) { - m_ptr = &endVal; + delete this->m_ptr; + this->m_ptr = nullptr; + if (this->i == 8) { + this->m_ptr = &endVal; return; } - xj = pixel.x + COL_OFFSETS[i]; - yj = pixel.y + ROW_OFFSETS[i]; - flow = (pixel.val >> (i * 4)) & 0xF; + xj = this->pixel.x + COL_OFFSETS[this->i]; + yj = this->pixel.y + ROW_OFFSETS[this->i]; + flow = (this->pixel.val >> (this->i * 4)) & 0xF; if (flow) { - m_ptr = new NeighborTuple(i, xj, yj, flow); - i++; + this->m_ptr = new NeighborTuple(this->i, xj, yj, static_cast(flow)); + this->i++; return; } else { - i++; + this->i++; next(); } } + + // Increments the pointer to the next downslope neighbor (D8) + template::value>* = nullptr> + void next() { + long xj, yj; + delete this->m_ptr; + this->m_ptr = nullptr; + + if (this->i == 8) { + this->m_ptr = &endVal; + return; + } + xj = this->pixel.x + COL_OFFSETS[this->pixel.val]; + yj = this->pixel.y + ROW_OFFSETS[this->pixel.val]; + this->i = 8; + this->m_ptr = new NeighborTuple(this->pixel.val, xj, yj, 1); + return; + } }; -class UpslopeNeighborIterator: public NeighborIterator { +// Iterates over neighbor pixels that are upslope of a given pixel, +// in either MFD or D8 mode +template +class UpslopeNeighborIterator: public NeighborIterator { public: - UpslopeNeighborIterator(): NeighborIterator() {} - UpslopeNeighborIterator(NeighborTuple* n): NeighborIterator(n) {} - UpslopeNeighborIterator(Pixel p) { - pixel = p; + + UpslopeNeighborIterator(): NeighborIterator() {} + UpslopeNeighborIterator(NeighborTuple* n): NeighborIterator(n) {} + UpslopeNeighborIterator(const Pixel p) { + this->pixel = p; next(); } + UpslopeNeighborIterator& operator++() { next(); return *this; } + UpslopeNeighborIterator operator++(int) { UpslopeNeighborIterator tmp = *this; ++(*this); return tmp; } + // Increments the pointer to the next upslope neighbor (MFD) + template::value>* = nullptr> void next() { long xj, yj; int flow_dir_j; int flow_ji; long flow_dir_j_sum; - delete m_ptr; - m_ptr = nullptr; - if (i == 8) { - m_ptr = &endVal; + delete this->m_ptr; + this->m_ptr = nullptr; + if (this->i == 8) { + this->m_ptr = &endVal; return; } - xj = pixel.x + COL_OFFSETS[i]; - yj = pixel.y + ROW_OFFSETS[i]; - if (xj < 0 or xj >= pixel.raster.raster_x_size or - yj < 0 or yj >= pixel.raster.raster_y_size) { - i++; + xj = this->pixel.x + COL_OFFSETS[this->i]; + yj = this->pixel.y + ROW_OFFSETS[this->i]; + if (xj < 0 or xj >= this->pixel.raster.raster_x_size or + yj < 0 or yj >= this->pixel.raster.raster_y_size) { + this->i++; next(); return; } - flow_dir_j = pixel.raster.get(xj, yj); - flow_ji = (0xF & (flow_dir_j >> (4 * FLOW_DIR_REVERSE_DIRECTION[i]))); + + flow_dir_j = static_cast(this->pixel.raster.get(xj, yj)); + flow_ji = (0xF & (flow_dir_j >> (4 * FLOW_DIR_REVERSE_DIRECTION[this->i]))); if (flow_ji) { flow_dir_j_sum = 0; for (int idx = 0; idx < 8; idx++) { flow_dir_j_sum += (flow_dir_j >> (idx * 4)) & 0xF; } - m_ptr = new NeighborTuple(i, xj, yj, static_cast(flow_ji) / static_cast(flow_dir_j_sum)); - i++; + this->m_ptr = new NeighborTuple( + this->i, xj, yj, + static_cast(flow_ji) / static_cast(flow_dir_j_sum)); + this->i++; return; } else { - i++; + this->i++; + next(); + } + } + + // Increments the pointer to the next upslope neighbor (D8) + template::value>* = nullptr> + void next() { + long xj, yj; + int flow_dir_j; + delete this->m_ptr; + this->m_ptr = nullptr; + if (this->i == 8) { + this->m_ptr = &endVal; + return; + } + xj = this->pixel.x + COL_OFFSETS[this->i]; + yj = this->pixel.y + ROW_OFFSETS[this->i]; + if (xj < 0 or xj >= this->pixel.raster.raster_x_size or + yj < 0 or yj >= this->pixel.raster.raster_y_size) { + this->i++; + next(); + return; + } + + flow_dir_j = static_cast(this->pixel.raster.get(xj, yj)); + if (flow_dir_j == FLOW_DIR_REVERSE_DIRECTION[this->i]) { + this->m_ptr = new NeighborTuple(this->i, xj, yj, 1); + this->i++; + return; + } else { + this->i++; next(); } } }; -class UpslopeNeighborNoDivideIterator: public NeighborIterator { + +// Iterates over neighbor pixels that are upslope of a given pixel, +// without dividing the flow_proportion, in either MFD or D8 mode +template +class UpslopeNeighborNoDivideIterator: public NeighborIterator { public: - UpslopeNeighborNoDivideIterator(): NeighborIterator() {} - UpslopeNeighborNoDivideIterator(NeighborTuple* n): NeighborIterator(n) {} - UpslopeNeighborNoDivideIterator(Pixel p) { - pixel = p; + + UpslopeNeighborNoDivideIterator(): NeighborIterator() {} + UpslopeNeighborNoDivideIterator(NeighborTuple* n): NeighborIterator(n) {} + UpslopeNeighborNoDivideIterator(const Pixel p) { + this->pixel = p; next(); } + UpslopeNeighborNoDivideIterator& operator++() { next(); return *this; } + UpslopeNeighborNoDivideIterator operator++(int) { UpslopeNeighborNoDivideIterator tmp = *this; ++(*this); return tmp; } + // Increments the pointer to the next upslope neighbor (MFD) + template::value>* = nullptr> void next() { long xj, yj; int flow_dir_j; int flow_ji; - long flow_dir_j_sum; - delete m_ptr; - m_ptr = nullptr; - if (i == 8) { - m_ptr = &endVal; + delete this->m_ptr; + this->m_ptr = nullptr; + if (this->i == 8) { + this->m_ptr = &endVal; return; } - xj = pixel.x + COL_OFFSETS[i]; - yj = pixel.y + ROW_OFFSETS[i]; - if (xj < 0 or xj >= pixel.raster.raster_x_size or - yj < 0 or yj >= pixel.raster.raster_y_size) { - i++; + xj = this->pixel.x + COL_OFFSETS[this->i]; + yj = this->pixel.y + ROW_OFFSETS[this->i]; + if (xj < 0 or xj >= this->pixel.raster.raster_x_size or + yj < 0 or yj >= this->pixel.raster.raster_y_size) { + this->i++; next(); return; } - flow_dir_j = pixel.raster.get(xj, yj); - flow_ji = (0xF & (flow_dir_j >> (4 * FLOW_DIR_REVERSE_DIRECTION[i]))); + flow_dir_j = static_cast(this->pixel.raster.get(xj, yj)); + flow_ji = (0xF & (flow_dir_j >> (4 * FLOW_DIR_REVERSE_DIRECTION[this->i]))); if (flow_ji) { - flow_dir_j_sum = 0; - for (int idx = 0; idx < 8; idx++) { - flow_dir_j_sum += (flow_dir_j >> (idx * 4)) & 0xF; - } - m_ptr = new NeighborTuple(i, xj, yj, flow_ji); - i++; + this->m_ptr = new NeighborTuple(this->i, xj, yj, static_cast(flow_ji)); + this->i++; + return; + } else { + this->i++; + next(); + } + } + + // Increments the pointer to the next upslope neighbor (D8) + template::value>* = nullptr> + void next() { + long xj, yj; + int flow_dir_j; + delete this->m_ptr; + this->m_ptr = nullptr; + if (this->i == 8) { + this->m_ptr = &endVal; + return; + } + xj = this->pixel.x + COL_OFFSETS[this->i]; + yj = this->pixel.y + ROW_OFFSETS[this->i]; + if (xj < 0 or xj >= this->pixel.raster.raster_x_size or + yj < 0 or yj >= this->pixel.raster.raster_y_size) { + this->i++; + next(); + return; + } + + flow_dir_j = static_cast(this->pixel.raster.get(xj, yj)); + if (flow_dir_j == FLOW_DIR_REVERSE_DIRECTION[this->i]) { + this->m_ptr = new NeighborTuple(this->i, xj, yj, 1); + this->i++; return; } else { - i++; + this->i++; next(); } } }; +template class Neighbors { public: - Pixel pixel; + Pixel pixel; Neighbors() {} - Neighbors(Pixel pixel): pixel(pixel) {} - NeighborIterator begin() { return NeighborIterator(pixel); } - NeighborIterator end() { return NeighborIterator(&endVal); } + Neighbors(const Pixel pixel): pixel(pixel) {} + NeighborIterator begin() { return NeighborIterator(pixel); } + NeighborIterator end() { return NeighborIterator(&endVal); } }; -class DownslopeNeighbors: public Neighbors { +template +class DownslopeNeighbors: public Neighbors { public: - using Neighbors::Neighbors; - DownslopeNeighborIterator begin() { return DownslopeNeighborIterator(pixel); } - DownslopeNeighborIterator end() { return DownslopeNeighborIterator(&endVal); } + using Neighbors::Neighbors; + DownslopeNeighborIterator begin() { return DownslopeNeighborIterator(this->pixel); } + DownslopeNeighborIterator end() { return DownslopeNeighborIterator(&endVal); } }; -class DownslopeNeighborsNoSkip: public Neighbors { +template +class DownslopeNeighborsNoSkip: public Neighbors { public: - using Neighbors::Neighbors; - DownslopeNeighborNoSkipIterator begin() { return DownslopeNeighborNoSkipIterator(pixel); } - DownslopeNeighborNoSkipIterator end() { return DownslopeNeighborNoSkipIterator(&endVal); } + using Neighbors::Neighbors; + DownslopeNeighborNoSkipIterator begin() { return DownslopeNeighborNoSkipIterator(this->pixel); } + DownslopeNeighborNoSkipIterator end() { return DownslopeNeighborNoSkipIterator(&endVal); } }; -class UpslopeNeighbors: public Neighbors { +template +class UpslopeNeighbors: public Neighbors { public: - using Neighbors::Neighbors; - UpslopeNeighborIterator begin() { return UpslopeNeighborIterator(pixel); } - UpslopeNeighborIterator end() { return UpslopeNeighborIterator(&endVal); } + using Neighbors::Neighbors; + UpslopeNeighborIterator begin() { return UpslopeNeighborIterator(this->pixel); } + UpslopeNeighborIterator end() { return UpslopeNeighborIterator(&endVal); } }; -class UpslopeNeighborsNoDivide: public Neighbors { +template +class UpslopeNeighborsNoDivide: public Neighbors { public: - using Neighbors::Neighbors; - UpslopeNeighborNoDivideIterator begin() { return UpslopeNeighborNoDivideIterator(pixel); } - UpslopeNeighborNoDivideIterator end() { return UpslopeNeighborNoDivideIterator(&endVal); } + using Neighbors::Neighbors; + UpslopeNeighborNoDivideIterator begin() { return UpslopeNeighborNoDivideIterator(this->pixel); } + UpslopeNeighborNoDivideIterator end() { return UpslopeNeighborNoDivideIterator(&endVal); } }; inline bool is_close(double x, double y) { @@ -647,3 +812,5 @@ inline bool is_close(double x, double y) { } return abs(x - y) <= (pow(10, -8) + pow(10, -05) * abs(y)); } + +#endif // NATCAP_INVEST_MANAGEDRASTER_H_ diff --git a/src/natcap/invest/managed_raster/managed_raster.pxd b/src/natcap/invest/managed_raster/managed_raster.pxd index f2d7b6d37..cb5a5c13c 100644 --- a/src/natcap/invest/managed_raster/managed_raster.pxd +++ b/src/natcap/invest/managed_raster/managed_raster.pxd @@ -3,6 +3,7 @@ from libcpp.list cimport list as clist from libcpp.pair cimport pair from libcpp.set cimport set as cset +from libcpp.stack cimport stack from libcpp.string cimport string from libcpp.vector cimport vector from libc.math cimport isnan @@ -36,6 +37,7 @@ cdef extern from "ManagedRaster.h": string raster_path int band_id int closed + double nodata ManagedRaster() except + ManagedRaster(char*, int, bool) except + @@ -44,7 +46,7 @@ cdef extern from "ManagedRaster.h": void _load_block(int block_index) except * void close() - cdef cppclass ManagedFlowDirRaster: + cdef cppclass ManagedFlowDirRaster[T]: LRUCache[int, double*]* lru_cache cset[int] dirty_blocks int block_xsize @@ -61,6 +63,7 @@ cdef extern from "ManagedRaster.h": string raster_path int band_id int closed + double nodata bint is_local_high_point(int xi, int yi) @@ -70,106 +73,60 @@ cdef extern from "ManagedRaster.h": double get(long xi, long yi) void close() + cdef cppclass D8 + cdef cppclass MFD + cdef cppclass NeighborTuple: NeighborTuple() except + NeighborTuple(int, int, int, float) except + int direction, x, y float flow_proportion - cdef cppclass UpslopeNeighborIteratorSkip: - ManagedFlowDirRaster raster + cdef cppclass DownslopeNeighborIterator[T]: + ManagedFlowDirRaster[T] raster int col int row int n_dir int flow_dir + int flow_dir_sum - UpslopeNeighborIteratorSkip() - UpslopeNeighborIteratorSkip(ManagedFlowDirRaster, int, int, int) - NeighborTuple next() - - cdef cppclass Pixel: - ManagedFlowDirRaster raster - int x - int y - int val - - Pixel() - Pixel(ManagedFlowDirRaster, int, int) - - cdef cppclass NeighborIterator: - NeighborIterator() - NeighborIterator(NeighborTuple* n) - NeighborIterator(Pixel) - NeighborTuple operator*() - NeighborIterator operator++() - bint operator==(NeighborIterator) - bint operator!=(NeighborIterator) - - cdef cppclass DownslopeNeighborIterator: DownslopeNeighborIterator() - DownslopeNeighborIterator(NeighborTuple* n) - DownslopeNeighborIterator(Pixel) - NeighborTuple operator*() - DownslopeNeighborIterator operator++() - bint operator==(DownslopeNeighborIterator) - bint operator!=(DownslopeNeighborIterator) - - cdef cppclass DownslopeNeighborNoSkipIterator: + DownslopeNeighborIterator(ManagedFlowDirRaster[T], int, int) + void next() + + cdef cppclass DownslopeNeighborNoSkipIterator[T]: + ManagedFlowDirRaster[T] raster + int col + int row + int n_dir + int flow_dir + int flow_dir_sum + DownslopeNeighborNoSkipIterator() - DownslopeNeighborNoSkipIterator(NeighborTuple* n) - DownslopeNeighborNoSkipIterator(Pixel) - NeighborTuple operator*() - DownslopeNeighborNoSkipIterator operator++() - bint operator==(DownslopeNeighborNoSkipIterator) - bint operator!=(DownslopeNeighborNoSkipIterator) - - cdef cppclass UpslopeNeighborIterator: + DownslopeNeighborNoSkipIterator(ManagedFlowDirRaster[T], int, int) + void next() + + cdef cppclass UpslopeNeighborIterator[T]: + ManagedFlowDirRaster[T] raster + int col + int row + int n_dir + int flow_dir + UpslopeNeighborIterator() - UpslopeNeighborIterator(NeighborTuple* n) - UpslopeNeighborIterator(Pixel) - NeighborTuple operator*() - UpslopeNeighborIterator operator++() - bint operator==(UpslopeNeighborIterator) - bint operator!=(UpslopeNeighborIterator) - - cdef cppclass UpslopeNeighborNoDivideIterator: + UpslopeNeighborIterator(ManagedFlowDirRaster[T], int, int) + void next() + + cdef cppclass UpslopeNeighborNoDivideIterator[T]: + ManagedFlowDirRaster[T] raster + int col + int row + int n_dir + int flow_dir + UpslopeNeighborNoDivideIterator() - UpslopeNeighborNoDivideIterator(NeighborTuple* n) - UpslopeNeighborNoDivideIterator(Pixel) - NeighborTuple operator*() - UpslopeNeighborNoDivideIterator operator++() - bint operator==(UpslopeNeighborNoDivideIterator) - bint operator!=(UpslopeNeighborNoDivideIterator) - - cdef cppclass Neighbors: - Neighbors() - Neighbors(Pixel) - NeighborIterator begin() - NeighborIterator end() - - cdef cppclass DownslopeNeighbors: - DownslopeNeighbors() - DownslopeNeighbors(Pixel) - DownslopeNeighborIterator begin() - DownslopeNeighborIterator end() - - cdef cppclass DownslopeNeighborsNoSkip: - DownslopeNeighborsNoSkip() - DownslopeNeighborsNoSkip(Pixel) - DownslopeNeighborNoSkipIterator begin() - DownslopeNeighborNoSkipIterator end() - - cdef cppclass UpslopeNeighbors: - UpslopeNeighbors() - UpslopeNeighbors(Pixel) - UpslopeNeighborIterator begin() - UpslopeNeighborIterator end() - - cdef cppclass UpslopeNeighborsNoDivide: - UpslopeNeighborsNoDivide() - UpslopeNeighborsNoDivide(Pixel) - UpslopeNeighborNoDivideIterator begin() - UpslopeNeighborNoDivideIterator end() + UpslopeNeighborNoDivideIterator(ManagedFlowDirRaster[T], int, int) + void next() bint is_close(double, double) diff --git a/src/natcap/invest/ndr/effective_retention.h b/src/natcap/invest/ndr/effective_retention.h new file mode 100644 index 000000000..bb5780f7b --- /dev/null +++ b/src/natcap/invest/ndr/effective_retention.h @@ -0,0 +1,262 @@ +#include "ManagedRaster.h" +#include +#include +#include // std::setprecision, std::setw +#include +#include + +// cdef extern from "time.h" nogil: +// ctypedef int time_t +// time_t time(time_t*) + +// LOGGER = logging.getLogger(__name__) + + +// Calculate flow downhill effective_retention to the channel. +// Args: +// flow_direction_path: a path to a flow direction raster (MFD or D8) +// stream_path: a path to a raster where 1 indicates a +// stream all other values ignored must be same dimensions and +// projection as flow_direction_path. +// retention_eff_lulc_path: a path to a raster indicating +// the maximum retention efficiency that the landcover on that +// pixel can accumulate. +// crit_len_path: a path to a raster indicating the critical +// length of the retention efficiency that the landcover on this +// pixel. +// effective_retention_path: path to a raster that is +// created by this call that contains a per-pixel effective +// sediment retention to the stream. +template +void run_effective_retention( + char* flow_direction_path, + char* stream_path, + char* retention_eff_lulc_path, + char* crit_len_path, + char* to_process_flow_directions_path, + char* effective_retention_path) { + // Within a stream, the effective retention is 0 + int STREAM_EFFECTIVE_RETENTION = 0; + float effective_retention_nodata = -1; + stack processing_stack; + + ManagedFlowDirRaster flow_dir_raster = ManagedFlowDirRaster( + flow_direction_path, 1, false); + ManagedRaster stream_raster = ManagedRaster(stream_path, 1, false); + ManagedRaster retention_eff_lulc_raster = ManagedRaster( + retention_eff_lulc_path, 1, false); + ManagedRaster crit_len_raster = ManagedRaster(crit_len_path, 1, false); + ManagedRaster to_process_flow_directions_raster = ManagedRaster( + to_process_flow_directions_path, 1, true); + ManagedRaster effective_retention_raster = ManagedRaster( + effective_retention_path, 1, true); + + long n_cols = flow_dir_raster.raster_x_size; + long n_rows = flow_dir_raster.raster_y_size; + // cell sizes must be square, so no reason to test at this point. + double cell_size = stream_raster.geotransform[1]; + + double crit_len_nodata = crit_len_raster.nodata; + double retention_eff_nodata = retention_eff_lulc_raster.nodata; + + long win_xsize, win_ysize, xoff, yoff; + long global_col, global_row; + unsigned long flat_index; + long flow_dir, neighbor_flow_dirs; + double current_step_factor, step_size, crit_len, retention_eff_lulc; + long neighbor_row, neighbor_col; + int neighbor_outflow_dir, neighbor_outflow_dir_mask, neighbor_process_flow_dir; + int outflow_dirs, dir_mask; + NeighborTuple neighbor; + bool should_seed; + double working_retention_eff; + DownslopeNeighborsNoSkip dn_neighbors; + UpslopeNeighbors up_neighbors; + bool has_outflow; + double neighbor_effective_retention; + double intermediate_retention; + string s; + long flow_dir_sum; + + // efficient way to calculate ceiling division: + // a divided by b rounded up = (a + (b - 1)) / b + // note that / represents integer floor division + // https://stackoverflow.com/a/62032709/14451410 + int n_col_blocks = (flow_dir_raster.raster_x_size + (flow_dir_raster.block_xsize - 1)) / flow_dir_raster.block_xsize; + int n_row_blocks = (flow_dir_raster.raster_y_size + (flow_dir_raster.block_ysize - 1)) / flow_dir_raster.block_ysize; + + for (int row_block_index = 0; row_block_index < n_row_blocks; row_block_index++) { + yoff = row_block_index * flow_dir_raster.block_ysize; + win_ysize = flow_dir_raster.raster_y_size - yoff; + if (win_ysize > flow_dir_raster.block_ysize) { + win_ysize = flow_dir_raster.block_ysize; + } + for (int col_block_index = 0; col_block_index < n_col_blocks; col_block_index++) { + xoff = col_block_index * flow_dir_raster.block_xsize; + win_xsize = flow_dir_raster.raster_x_size - xoff; + if (win_xsize > flow_dir_raster.block_xsize) { + win_xsize = flow_dir_raster.block_xsize; + } + + // if ctime(NULL) - last_log_time > 5.0: + // last_log_time = ctime(NULL) + // LOGGER.info('Sediment deposition %.2f%% complete', 100 * ( + // n_pixels_processed / float(flow_dir_raster.raster_x_size * flow_dir_raster.raster_y_size))) + + for (int row_index = 0; row_index < win_ysize; row_index++) { + global_row = yoff + row_index; + for (int col_index = 0; col_index < win_xsize; col_index++) { + global_col = xoff + col_index; + outflow_dirs = int(to_process_flow_directions_raster.get( + global_col, global_row)); + should_seed = false; + // # see if this pixel drains to nodata or the edge, if so it's + // # a drain + for (int i = 0; i < 8; i++) { + dir_mask = 1 << i; + if ((outflow_dirs & dir_mask) > 0) { + neighbor_col = COL_OFFSETS[i] + global_col; + neighbor_row = ROW_OFFSETS[i] + global_row; + if (neighbor_col < 0 or neighbor_col >= n_cols or + neighbor_row < 0 or neighbor_row >= n_rows) { + should_seed = true; + outflow_dirs &= ~dir_mask; + } else { + // Only consider neighbor flow directions if the + // neighbor index is within the raster. + neighbor_flow_dirs = long( + to_process_flow_directions_raster.get( + neighbor_col, neighbor_row)); + if (neighbor_flow_dirs == 0) { + should_seed = true; + outflow_dirs &= ~dir_mask; + } + } + } + } + + if (should_seed) { + // mark all outflow directions processed + to_process_flow_directions_raster.set( + global_col, global_row, outflow_dirs); + processing_stack.push(global_row * n_cols + global_col); + } + } + } + + while (processing_stack.size() > 0) { + // loop invariant, we don't push a cell on the stack that + // hasn't already been set for processing. + flat_index = processing_stack.top(); + processing_stack.pop(); + global_row = flat_index / n_cols; // integer floor division + global_col = flat_index % n_cols; + + crit_len = crit_len_raster.get(global_col, global_row); + retention_eff_lulc = retention_eff_lulc_raster.get(global_col, global_row); + flow_dir = int(flow_dir_raster.get(global_col, global_row)); + if (stream_raster.get(global_col, global_row) == 1) { + // if it's a stream, effective retention is 0. + effective_retention_raster.set(global_col, global_row, STREAM_EFFECTIVE_RETENTION); + } else if (is_close(crit_len, crit_len_nodata) or + is_close(retention_eff_lulc, retention_eff_nodata) or + flow_dir == 0) { + // if it's nodata, effective retention is nodata. + effective_retention_raster.set( + global_col, global_row, effective_retention_nodata); + } else { + working_retention_eff = 0; + + dn_neighbors = DownslopeNeighborsNoSkip( + Pixel(flow_dir_raster, global_col, global_row)); + has_outflow = false; + flow_dir_sum = 0; + for (auto neighbor: dn_neighbors) { + has_outflow = true; + flow_dir_sum += static_cast(neighbor.flow_proportion); + if (neighbor.x < 0 or neighbor.x >= n_cols or + neighbor.y < 0 or neighbor.y >= n_rows) { + continue; + } + if (neighbor.direction % 2 == 1) { + step_size = cell_size * 1.41421356237; + } else { + step_size = cell_size; + } + // guard against a critical length factor that's 0 + if (crit_len > 0) { + current_step_factor = exp(-5 * step_size / crit_len); + } else { + current_step_factor = 0; + } + + neighbor_effective_retention = ( + effective_retention_raster.get( + neighbor.x, neighbor.y)); + + // Case 1: downslope neighbor is a stream pixel + if (neighbor_effective_retention == STREAM_EFFECTIVE_RETENTION) { + intermediate_retention = ( + retention_eff_lulc * (1 - current_step_factor)); + // Case 2: the current LULC's retention exceeds the neighbor's retention. + } else if (retention_eff_lulc > neighbor_effective_retention) { + intermediate_retention = ( + (neighbor_effective_retention * current_step_factor) + + (retention_eff_lulc * (1 - current_step_factor))); + // Case 3: the other 2 cases have not been hit. + } else { + intermediate_retention = neighbor_effective_retention; + } + + working_retention_eff += ( + intermediate_retention * neighbor.flow_proportion); + } + + if (has_outflow) { + double v = working_retention_eff / flow_dir_sum; + effective_retention_raster.set( + global_col, global_row, v); + } else { + throw std::logic_error( + "got to a cell that has no outflow! This error is happening" + "in effective_retention.h"); + } + } + // search upslope to see if we need to push a cell on the stack + // for i in range(8): + up_neighbors = UpslopeNeighbors(Pixel(flow_dir_raster, global_col, global_row)); + for (auto neighbor: up_neighbors) { + neighbor_outflow_dir = INFLOW_OFFSETS[neighbor.direction]; + neighbor_outflow_dir_mask = 1 << neighbor_outflow_dir; + neighbor_process_flow_dir = int( + to_process_flow_directions_raster.get( + neighbor.x, neighbor.y)); + if (neighbor_process_flow_dir == 0) { + // skip, due to loop invariant this must be a nodata pixel + continue; + } + if ((neighbor_process_flow_dir & neighbor_outflow_dir_mask )== 0) { + // no outflow + continue; + } + // mask out the outflow dir that this iteration processed + neighbor_process_flow_dir &= ~neighbor_outflow_dir_mask; + to_process_flow_directions_raster.set( + neighbor.x, neighbor.y, neighbor_process_flow_dir); + if (neighbor_process_flow_dir == 0) { + // if 0 then all downslope have been processed, + // push on stack, otherwise another downslope pixel will + // pick it up + processing_stack.push(neighbor.y * n_cols + neighbor.x); + } + } + } + } + } + stream_raster.close(); + crit_len_raster.close(); + retention_eff_lulc_raster.close(); + effective_retention_raster.close(); + flow_dir_raster.close(); + to_process_flow_directions_raster.close(); +} diff --git a/src/natcap/invest/ndr/effective_retention.pxd b/src/natcap/invest/ndr/effective_retention.pxd new file mode 100644 index 000000000..1d2798ccb --- /dev/null +++ b/src/natcap/invest/ndr/effective_retention.pxd @@ -0,0 +1,8 @@ +cdef extern from "effective_retention.h": + void run_effective_retention[T]( + char*, + char*, + char*, + char*, + char*, + char*) diff --git a/src/natcap/invest/ndr/ndr.py b/src/natcap/invest/ndr/ndr.py index 4958a6252..f9c9cd0ca 100644 --- a/src/natcap/invest/ndr/ndr.py +++ b/src/natcap/invest/ndr/ndr.py @@ -162,6 +162,21 @@ "reached through subsurface flow. This characterizes the " "retention due to biochemical degradation in soils. Required " "if Calculate Nitrogen is selected.") + }, + "algorithm": { + "type": "option_string", + "options": { + "D8": { + "display_name": gettext("D8"), + "description": "D8 flow direction" + }, + "MFD": { + "display_name": gettext("MFD"), + "description": "Multiple flow direction" + } + }, + "about": gettext("Flow direction algorithm to use."), + "name": gettext("flow direction algorithm") } }, "outputs": { @@ -668,35 +683,6 @@ def execute(args): target_path_list=[f_reg['filled_dem_path']], task_name='fill pits') - flow_dir_task = task_graph.add_task( - func=pygeoprocessing.routing.flow_dir_mfd, - args=( - (f_reg['filled_dem_path'], 1), f_reg['flow_direction_path']), - kwargs={'working_dir': intermediate_output_dir}, - dependent_task_list=[fill_pits_task], - target_path_list=[f_reg['flow_direction_path']], - task_name='flow dir') - - flow_accum_task = task_graph.add_task( - func=pygeoprocessing.routing.flow_accumulation_mfd, - args=( - (f_reg['flow_direction_path'], 1), - f_reg['flow_accumulation_path']), - target_path_list=[f_reg['flow_accumulation_path']], - dependent_task_list=[flow_dir_task], - task_name='flow accum') - - stream_extraction_task = task_graph.add_task( - func=pygeoprocessing.routing.extract_streams_mfd, - args=( - (f_reg['flow_accumulation_path'], 1), - (f_reg['flow_direction_path'], 1), - float(args['threshold_flow_accumulation']), - f_reg['stream_path']), - target_path_list=[f_reg['stream_path']], - dependent_task_list=[flow_accum_task], - task_name='stream extraction') - calculate_slope_task = task_graph.add_task( func=pygeoprocessing.calculate_slope, args=((f_reg['filled_dem_path'], 1), f_reg['slope_path']), @@ -714,6 +700,81 @@ def execute(args): dependent_task_list=[calculate_slope_task], task_name='threshold slope') + if args['algorithm'] == 'MFD': + flow_dir_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_dir_mfd, + args=( + (f_reg['filled_dem_path'], 1), f_reg['flow_direction_path']), + kwargs={'working_dir': intermediate_output_dir}, + dependent_task_list=[fill_pits_task], + target_path_list=[f_reg['flow_direction_path']], + task_name='flow dir') + + flow_accum_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_accumulation_mfd, + args=( + (f_reg['flow_direction_path'], 1), + f_reg['flow_accumulation_path']), + target_path_list=[f_reg['flow_accumulation_path']], + dependent_task_list=[flow_dir_task], + task_name='flow accum') + + stream_extraction_task = task_graph.add_task( + func=pygeoprocessing.routing.extract_streams_mfd, + args=( + (f_reg['flow_accumulation_path'], 1), + (f_reg['flow_direction_path'], 1), + float(args['threshold_flow_accumulation']), + f_reg['stream_path']), + target_path_list=[f_reg['stream_path']], + dependent_task_list=[flow_accum_task], + task_name='stream extraction') + s_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_accumulation_mfd, + args=((f_reg['flow_direction_path'], 1), f_reg['s_accumulation_path']), + kwargs={ + 'weight_raster_path_band': (f_reg['thresholded_slope_path'], 1)}, + target_path_list=[f_reg['s_accumulation_path']], + dependent_task_list=[flow_dir_task, threshold_slope_task], + task_name='route s') + else: # D8 + flow_dir_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_dir_d8, + args=( + (f_reg['filled_dem_path'], 1), f_reg['flow_direction_path']), + kwargs={'working_dir': intermediate_output_dir}, + dependent_task_list=[fill_pits_task], + target_path_list=[f_reg['flow_direction_path']], + task_name='flow dir') + + flow_accum_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_accumulation_d8, + args=( + (f_reg['flow_direction_path'], 1), + f_reg['flow_accumulation_path']), + target_path_list=[f_reg['flow_accumulation_path']], + dependent_task_list=[flow_dir_task], + task_name='flow accum') + + stream_extraction_task = task_graph.add_task( + func=pygeoprocessing.routing.extract_streams_d8, + kwargs=dict( + flow_accum_raster_path_band=(f_reg['flow_accumulation_path'], 1), + flow_threshold=float(args['threshold_flow_accumulation']), + target_stream_raster_path=f_reg['stream_path']), + target_path_list=[f_reg['stream_path']], + dependent_task_list=[flow_accum_task], + task_name='stream extraction') + + s_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_accumulation_d8, + args=((f_reg['flow_direction_path'], 1), f_reg['s_accumulation_path']), + kwargs={ + 'weight_raster_path_band': (f_reg['thresholded_slope_path'], 1)}, + target_path_list=[f_reg['s_accumulation_path']], + dependent_task_list=[flow_dir_task, threshold_slope_task], + task_name='route s') + runoff_proxy_index_task = task_graph.add_task( func=_normalize_raster, args=((f_reg['masked_runoff_proxy_path'], 1), @@ -722,15 +783,6 @@ def execute(args): dependent_task_list=[align_raster_task, mask_runoff_proxy_task], task_name='runoff proxy mean') - s_task = task_graph.add_task( - func=pygeoprocessing.routing.flow_accumulation_mfd, - args=((f_reg['flow_direction_path'], 1), f_reg['s_accumulation_path']), - kwargs={ - 'weight_raster_path_band': (f_reg['thresholded_slope_path'], 1)}, - target_path_list=[f_reg['s_accumulation_path']], - dependent_task_list=[flow_dir_task, threshold_slope_task], - task_name='route s') - s_bar_task = task_graph.add_task( func=pygeoprocessing.raster_map, kwargs=dict( @@ -762,27 +814,50 @@ def execute(args): dependent_task_list=[threshold_slope_task], task_name='s inv') - d_dn_task = task_graph.add_task( - func=pygeoprocessing.routing.distance_to_channel_mfd, - args=( - (f_reg['flow_direction_path'], 1), - (f_reg['stream_path'], 1), - f_reg['d_dn_path']), - kwargs={'weight_raster_path_band': ( - f_reg['s_factor_inverse_path'], 1)}, - dependent_task_list=[stream_extraction_task, s_inv_task], - target_path_list=[f_reg['d_dn_path']], - task_name='d dn') - - dist_to_channel_task = task_graph.add_task( - func=pygeoprocessing.routing.distance_to_channel_mfd, - args=( - (f_reg['flow_direction_path'], 1), - (f_reg['stream_path'], 1), - f_reg['dist_to_channel_path']), - dependent_task_list=[stream_extraction_task], - target_path_list=[f_reg['dist_to_channel_path']], - task_name='dist to channel') + if args['algorithm'] == 'MFD': + d_dn_task = task_graph.add_task( + func=pygeoprocessing.routing.distance_to_channel_mfd, + args=( + (f_reg['flow_direction_path'], 1), + (f_reg['stream_path'], 1), + f_reg['d_dn_path']), + kwargs={'weight_raster_path_band': ( + f_reg['s_factor_inverse_path'], 1)}, + dependent_task_list=[stream_extraction_task, s_inv_task], + target_path_list=[f_reg['d_dn_path']], + task_name='d dn') + + dist_to_channel_task = task_graph.add_task( + func=pygeoprocessing.routing.distance_to_channel_mfd, + args=( + (f_reg['flow_direction_path'], 1), + (f_reg['stream_path'], 1), + f_reg['dist_to_channel_path']), + dependent_task_list=[stream_extraction_task], + target_path_list=[f_reg['dist_to_channel_path']], + task_name='dist to channel') + else: # D8 + d_dn_task = task_graph.add_task( + func=pygeoprocessing.routing.distance_to_channel_d8, + args=( + (f_reg['flow_direction_path'], 1), + (f_reg['stream_path'], 1), + f_reg['d_dn_path']), + kwargs={'weight_raster_path_band': ( + f_reg['s_factor_inverse_path'], 1)}, + dependent_task_list=[stream_extraction_task, s_inv_task], + target_path_list=[f_reg['d_dn_path']], + task_name='d dn') + + dist_to_channel_task = task_graph.add_task( + func=pygeoprocessing.routing.distance_to_channel_d8, + args=( + (f_reg['flow_direction_path'], 1), + (f_reg['stream_path'], 1), + f_reg['dist_to_channel_path']), + dependent_task_list=[stream_extraction_task], + target_path_list=[f_reg['dist_to_channel_path']], + task_name='dist to channel') _ = task_graph.add_task( func=sdr._calculate_what_drains_to_stream, @@ -869,7 +944,8 @@ def execute(args): args=( f_reg['flow_direction_path'], f_reg['stream_path'], eff_path, - crit_len_path, effective_retention_path), + crit_len_path, effective_retention_path, + args['algorithm']), target_path_list=[effective_retention_path], dependent_task_list=[ stream_extraction_task, eff_task, crit_len_task], diff --git a/src/natcap/invest/ndr/ndr_core.pyx b/src/natcap/invest/ndr/ndr_core.pyx index 2f67113be..6cdffb2df 100644 --- a/src/natcap/invest/ndr/ndr_core.pyx +++ b/src/natcap/invest/ndr/ndr_core.pyx @@ -1,7 +1,6 @@ import tempfile import logging import os -import collections import numpy import pygeoprocessing @@ -9,20 +8,9 @@ cimport numpy cimport cython from osgeo import gdal -from libcpp.stack cimport stack -from libc.math cimport exp -from ..managed_raster.managed_raster cimport ManagedRaster -from ..managed_raster.managed_raster cimport ManagedFlowDirRaster -from ..managed_raster.managed_raster cimport DownslopeNeighborsNoSkip -from ..managed_raster.managed_raster cimport Pixel -from ..managed_raster.managed_raster cimport UpslopeNeighbors -from ..managed_raster.managed_raster cimport NeighborTuple -from ..managed_raster.managed_raster cimport is_close -from ..managed_raster.managed_raster cimport INFLOW_OFFSETS -from ..managed_raster.managed_raster cimport COL_OFFSETS -from ..managed_raster.managed_raster cimport ROW_OFFSETS - -from cython.operator import dereference +from ..managed_raster.managed_raster cimport D8 +from ..managed_raster.managed_raster cimport MFD +from .effective_retention cimport run_effective_retention cdef extern from "time.h" nogil: ctypedef int time_t @@ -34,16 +22,16 @@ LOGGER = logging.getLogger(__name__) cdef int STREAM_EFFECTIVE_RETENTION = 0 def ndr_eff_calculation( - mfd_flow_direction_path, stream_path, retention_eff_lulc_path, - crit_len_path, effective_retention_path): + flow_direction_path, stream_path, retention_eff_lulc_path, + crit_len_path, effective_retention_path, algorithm): """Calculate flow downhill effective_retention to the channel. Args: - mfd_flow_direction_path (string): a path to a raster with - pygeoprocessing.routing MFD flow direction values. + flow_direction_path (string): a path to a raster with + pygeoprocessing.routing flow direction values (MFD or D8). stream_path (string): a path to a raster where 1 indicates a stream all other values ignored must be same dimensions and - projection as mfd_flow_direction_path. + projection as flow_direction_path. retention_eff_lulc_path (string): a path to a raster indicating the maximum retention efficiency that the landcover on that pixel can accumulate. @@ -60,38 +48,13 @@ def ndr_eff_calculation( """ cdef float effective_retention_nodata = -1.0 pygeoprocessing.new_raster_from_base( - mfd_flow_direction_path, effective_retention_path, gdal.GDT_Float32, + flow_direction_path, effective_retention_path, gdal.GDT_Float32, [effective_retention_nodata]) fp, to_process_flow_directions_path = tempfile.mkstemp( suffix='.tif', prefix='flow_to_process', dir=os.path.dirname(effective_retention_path)) os.close(fp) - cdef long n_cols, n_rows - flow_dir_info = pygeoprocessing.get_raster_info(mfd_flow_direction_path) - n_cols, n_rows = flow_dir_info['raster_size'] - - cdef stack[long] processing_stack - stream_info = pygeoprocessing.get_raster_info(stream_path) - # cell sizes must be square, so no reason to test at this point. - cdef float cell_size = abs(stream_info['pixel_size'][0]) - - cdef ManagedRaster stream_raster = ManagedRaster( - stream_path.encode('utf-8'), 1, False) - cdef ManagedRaster crit_len_raster = ManagedRaster( - crit_len_path.encode('utf-8'), 1, False) - cdef ManagedRaster retention_eff_lulc_raster = ManagedRaster( - retention_eff_lulc_path.encode('utf-8'), 1, False) - cdef ManagedRaster effective_retention_raster = ManagedRaster( - effective_retention_path.encode('utf-8'), 1, True) - cdef ManagedFlowDirRaster mfd_flow_direction_raster = ManagedFlowDirRaster( - mfd_flow_direction_path.encode('utf-8'), 1, False) - - cdef float crit_len_nodata = pygeoprocessing.get_raster_info( - crit_len_path)['nodata'][0] - cdef float retention_eff_nodata = pygeoprocessing.get_raster_info( - retention_eff_lulc_path)['nodata'][0] - # create direction raster in bytes def _mfd_to_flow_dir_op(mfd_array): result = numpy.zeros(mfd_array.shape, dtype=numpy.uint8) @@ -99,172 +62,32 @@ def ndr_eff_calculation( result[:] |= ((((mfd_array >> (i*4)) & 0xF) > 0) << i).astype(numpy.uint8) return result + # create direction raster in bytes + def _d8_to_flow_dir_op(d8_array): + result = numpy.zeros(d8_array.shape, dtype=numpy.uint8) + for i in range(8): + result[d8_array == i] = 1 << i + return result + + flow_dir_op = _mfd_to_flow_dir_op if algorithm == 'MFD' else _d8_to_flow_dir_op + # convert mfd raster to binary mfd # each value is an 8-digit binary number # where 1 indicates that the pixel drains in that direction # and 0 indicates that it does not drain in that direction pygeoprocessing.raster_calculator( - [(mfd_flow_direction_path, 1)], _mfd_to_flow_dir_op, + [(flow_direction_path, 1)], flow_dir_op, to_process_flow_directions_path, gdal.GDT_Byte, None) - cdef ManagedRaster to_process_flow_directions_raster = ManagedRaster( - to_process_flow_directions_path.encode('utf-8'), 1, True) - - cdef long col_index, row_index, win_xsize, win_ysize, xoff, yoff - cdef long global_col, global_row - cdef unsigned long flat_index - cdef long outflow_weight, flow_dir - cdef long ds_col, ds_row, i - cdef float current_step_factor, step_size, crit_len, flow_dir_sum - cdef long neighbor_row, neighbor_col - cdef int neighbor_outflow_dir, neighbor_outflow_dir_mask, neighbor_process_flow_dir - cdef int outflow_dirs, dir_mask - cdef NeighborTuple neighbor - - for offset_dict in pygeoprocessing.iterblocks( - (mfd_flow_direction_path, 1), offset_only=True, largest_block=0): - # use cython variables to avoid python overhead of dict values - win_xsize = offset_dict['win_xsize'] - win_ysize = offset_dict['win_ysize'] - xoff = offset_dict['xoff'] - yoff = offset_dict['yoff'] - for row_index in range(win_ysize): - global_row = yoff + row_index - for col_index in range(win_xsize): - global_col = xoff + col_index - outflow_dirs = to_process_flow_directions_raster.get( - global_col, global_row) - should_seed = 0 - # see if this pixel drains to nodata or the edge, if so it's - # a drain - for i in range(8): - dir_mask = 1 << i - if outflow_dirs & dir_mask > 0: - neighbor_col = COL_OFFSETS[i] + global_col - neighbor_row = ROW_OFFSETS[i] + global_row - if (neighbor_col < 0 or neighbor_col >= n_cols or - neighbor_row < 0 or neighbor_row >= n_rows): - should_seed = 1 - outflow_dirs &= ~dir_mask - - # Only consider neighbor flow directions if the - # neighbor index is within the raster. - else: - neighbor_flow_dirs = ( - to_process_flow_directions_raster.get( - neighbor_col, neighbor_row)) - if neighbor_flow_dirs == 0: - should_seed = 1 - outflow_dirs &= ~dir_mask - - if should_seed: - # mark all outflow directions processed - to_process_flow_directions_raster.set( - global_col, global_row, outflow_dirs) - processing_stack.push(global_row * n_cols + global_col) - - while processing_stack.size() > 0: - # loop invariant, we don't push a cell on the stack that - # hasn't already been set for processing. - flat_index = processing_stack.top() - processing_stack.pop() - global_row = flat_index // n_cols - global_col = flat_index % n_cols - - crit_len = crit_len_raster.get(global_col, global_row) - retention_eff_lulc = retention_eff_lulc_raster.get( - global_col, global_row) - flow_dir = mfd_flow_direction_raster.get( - global_col, global_row) - if stream_raster.get(global_col, global_row) == 1: - # if it's a stream, effective retention is 0. - effective_retention_raster.set(global_col, global_row, STREAM_EFFECTIVE_RETENTION) - elif (is_close(crit_len, crit_len_nodata) or - is_close(retention_eff_lulc, retention_eff_nodata) or - flow_dir == 0): - # if it's nodata, effective retention is nodata. - effective_retention_raster.set( - global_col, global_row, effective_retention_nodata) - else: - working_retention_eff = 0.0 - dn_neighbors = DownslopeNeighborsNoSkip( - Pixel(mfd_flow_direction_raster, global_col, global_row)) - has_outflow = False - flow_dir_sum = 0 - for neighbor in dn_neighbors: - has_outflow = True - flow_dir_sum += neighbor.flow_proportion - if (neighbor.x < 0 or neighbor.x >= n_cols or - neighbor.y < 0 or neighbor.y >= n_rows): - continue - if neighbor.direction % 2 == 1: - step_size = (cell_size * 1.41421356237) - else: - step_size = cell_size - # guard against a critical length factor that's 0 - if crit_len > 0: - current_step_factor = ( - exp(-5 * step_size / crit_len)) - else: - current_step_factor = 0.0 - - neighbor_effective_retention = ( - effective_retention_raster.get( - neighbor.x, neighbor.y)) - - # Case 1: downslope neighbor is a stream pixel - if neighbor_effective_retention == STREAM_EFFECTIVE_RETENTION: - intermediate_retention = ( - retention_eff_lulc * (1 - current_step_factor)) - - # Case 2: the current LULC's retention exceeds the neighbor's retention. - elif retention_eff_lulc > neighbor_effective_retention: - intermediate_retention = ( - (neighbor_effective_retention * current_step_factor) + - (retention_eff_lulc * (1 - current_step_factor))) - - # Case 3: the other 2 cases have not been hit. - else: - intermediate_retention = neighbor_effective_retention - - working_retention_eff += ( - intermediate_retention * neighbor.flow_proportion) - - if has_outflow: - working_retention_eff /= flow_dir_sum - effective_retention_raster.set( - global_col, global_row, working_retention_eff) - else: - raise Exception("got to a cell that has no outflow!") - # search upslope to see if we need to push a cell on the stack - # for i in range(8): - up_neighbors = UpslopeNeighbors( - Pixel(mfd_flow_direction_raster, global_col, global_row)) - for neighbor in up_neighbors: - neighbor_outflow_dir = INFLOW_OFFSETS[neighbor.direction] - neighbor_outflow_dir_mask = 1 << neighbor_outflow_dir - neighbor_process_flow_dir = ( - to_process_flow_directions_raster.get( - neighbor.x, neighbor.y)) - if neighbor_process_flow_dir == 0: - # skip, due to loop invariant this must be a nodata pixel - continue - if neighbor_process_flow_dir & neighbor_outflow_dir_mask == 0: - # no outflow - continue - # mask out the outflow dir that this iteration processed - neighbor_process_flow_dir &= ~neighbor_outflow_dir_mask - to_process_flow_directions_raster.set( - neighbor.x, neighbor.y, neighbor_process_flow_dir) - if neighbor_process_flow_dir == 0: - # if 0 then all downslope have been processed, - # push on stack, otherwise another downslope pixel will - # pick it up - processing_stack.push(neighbor.y * n_cols + neighbor.x) - - stream_raster.close() - crit_len_raster.close() - retention_eff_lulc_raster.close() - effective_retention_raster.close() - mfd_flow_direction_raster.close() - to_process_flow_directions_raster.close() + args = [ + flow_direction_path.encode('utf-8'), + stream_path.encode('utf-8'), + retention_eff_lulc_path.encode('utf-8'), + crit_len_path.encode('utf-8'), + to_process_flow_directions_path.encode('utf-8'), + effective_retention_path.encode('utf-8')] + + if algorithm == 'MFD': + run_effective_retention[MFD](*args) + else: # D8 + run_effective_retention[D8](*args) diff --git a/src/natcap/invest/sdr/sdr.py b/src/natcap/invest/sdr/sdr.py index f1540982b..6a5f5c144 100644 --- a/src/natcap/invest/sdr/sdr.py +++ b/src/natcap/invest/sdr/sdr.py @@ -140,6 +140,21 @@ "watershed. Pixels with 1 are drainages and are treated like " "streams. Pixels with 0 are not drainages."), "name": gettext("drainages") + }, + "algorithm": { + "type": "option_string", + "options": { + "D8": { + "display_name": gettext("D8"), + "description": "D8 flow direction" + }, + "MFD": { + "display_name": gettext("MFD"), + "description": "Multiple flow direction" + } + }, + "about": gettext("Flow direction algorithm to use."), + "name": gettext("flow direction algorithm") } }, "outputs": { @@ -681,23 +696,68 @@ def execute(args): dependent_task_list=[slope_task], task_name='threshold slope') - flow_dir_task = task_graph.add_task( - func=pygeoprocessing.routing.flow_dir_mfd, - args=( - (f_reg['pit_filled_dem_path'], 1), - f_reg['flow_direction_path']), - target_path_list=[f_reg['flow_direction_path']], - dependent_task_list=[pit_fill_task], - task_name='flow direction calculation') + if args['algorithm'] == 'MFD': + flow_dir_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_dir_mfd, + args=( + (f_reg['pit_filled_dem_path'], 1), + f_reg['flow_direction_path']), + target_path_list=[f_reg['flow_direction_path']], + dependent_task_list=[pit_fill_task], + task_name='flow direction calculation') + + flow_accumulation_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_accumulation_mfd, + args=( + (f_reg['flow_direction_path'], 1), + f_reg['flow_accumulation_path']), + target_path_list=[f_reg['flow_accumulation_path']], + dependent_task_list=[flow_dir_task], + task_name='flow accumulation calculation') + + stream_task = task_graph.add_task( + func=pygeoprocessing.routing.extract_streams_mfd, + args=( + (f_reg['flow_accumulation_path'], 1), + (f_reg['flow_direction_path'], 1), + float(args['threshold_flow_accumulation']), + f_reg['stream_path']), + kwargs={'trace_threshold_proportion': 0.7}, + target_path_list=[f_reg['stream_path']], + dependent_task_list=[flow_accumulation_task], + task_name='extract streams') + + d_dn_func = pygeoprocessing.routing.distance_to_channel_mfd + else: - flow_accumulation_task = task_graph.add_task( - func=pygeoprocessing.routing.flow_accumulation_mfd, - args=( - (f_reg['flow_direction_path'], 1), - f_reg['flow_accumulation_path']), - target_path_list=[f_reg['flow_accumulation_path']], - dependent_task_list=[flow_dir_task], - task_name='flow accumulation calculation') + flow_dir_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_dir_d8, + args=( + (f_reg['pit_filled_dem_path'], 1), + f_reg['flow_direction_path']), + target_path_list=[f_reg['flow_direction_path']], + dependent_task_list=[pit_fill_task], + task_name='flow direction calculation') + + flow_accumulation_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_accumulation_d8, + args=( + (f_reg['flow_direction_path'], 1), + f_reg['flow_accumulation_path']), + target_path_list=[f_reg['flow_accumulation_path']], + dependent_task_list=[flow_dir_task], + task_name='flow accumulation calculation') + + stream_task = task_graph.add_task( + func=pygeoprocessing.routing.extract_streams_d8, + kwargs=dict( + flow_accum_raster_path_band=(f_reg['flow_accumulation_path'], 1), + flow_threshold=float(args['threshold_flow_accumulation']), + target_stream_raster_path=f_reg['stream_path']), + target_path_list=[f_reg['stream_path']], + dependent_task_list=[flow_accumulation_task], + task_name='extract streams') + d_dn_func = pygeoprocessing.routing.distance_to_channel_d8 ls_factor_task = task_graph.add_task( func=_calculate_ls_factor, @@ -711,18 +771,6 @@ def execute(args): flow_accumulation_task, slope_task], task_name='ls factor calculation') - stream_task = task_graph.add_task( - func=pygeoprocessing.routing.extract_streams_mfd, - args=( - (f_reg['flow_accumulation_path'], 1), - (f_reg['flow_direction_path'], 1), - float(args['threshold_flow_accumulation']), - f_reg['stream_path']), - kwargs={'trace_threshold_proportion': 0.7}, - target_path_list=[f_reg['stream_path']], - dependent_task_list=[flow_accumulation_task], - task_name='extract streams') - if drainage_present: drainage_task = task_graph.add_task( func=pygeoprocessing.raster_map, @@ -796,14 +844,17 @@ def execute(args): 's_bar')]: bar_task = task_graph.add_task( func=_calculate_bar_factor, - args=( - f_reg['flow_direction_path'], factor_path, - f_reg['flow_accumulation_path'], - accumulation_path, out_bar_path), + kwargs=dict( + flow_direction_path=f_reg['flow_direction_path'], + factor_path=factor_path, + flow_accumulation_path=f_reg['flow_accumulation_path'], + accumulation_path=accumulation_path, + out_bar_path=out_bar_path, + algorithm=args['algorithm']), target_path_list=[accumulation_path, out_bar_path], dependent_task_list=[ factor_task, flow_accumulation_task, flow_dir_task], - task_name='calculate %s' % bar_id) + task_name=f'calculate {bar_id}') bar_task_map[bar_id] = bar_task d_up_task = task_graph.add_task( @@ -829,7 +880,7 @@ def execute(args): task_name='calculate inverse ws factor') d_dn_task = task_graph.add_task( - func=pygeoprocessing.routing.distance_to_channel_mfd, + func=d_dn_func, args=( (f_reg['flow_direction_path'], 1), (drainage_raster_path_task[0], 1), @@ -880,10 +931,13 @@ def execute(args): sed_deposition_task = task_graph.add_task( func=sdr_core.calculate_sediment_deposition, - args=( - f_reg['flow_direction_path'], f_reg['e_prime_path'], - f_reg['f_path'], f_reg['sdr_path'], - f_reg['sed_deposition_path']), + kwargs=dict( + flow_direction_path=f_reg['flow_direction_path'], + e_prime_path=f_reg['e_prime_path'], + f_path=f_reg['f_path'], + sdr_path=f_reg['sdr_path'], + target_sediment_deposition_path=f_reg['sed_deposition_path'], + algorithm=args['algorithm']), dependent_task_list=[e_prime_task, sdr_task, flow_dir_task], target_path_list=[f_reg['sed_deposition_path'], f_reg['f_path']], task_name='sediment deposition') @@ -982,55 +1036,53 @@ def inverse_ws_op(w_factor, s_factor): return 1 / (w_factor * s_factor) def _calculate_what_drains_to_stream( - flow_dir_mfd_path, dist_to_channel_mfd_path, target_mask_path): + flow_dir_path, dist_to_channel_path, target_mask_path): """Create a mask indicating regions that do or do not drain to a stream. - This is useful because ``pygeoprocessing.distance_to_stream_mfd`` may leave + This is useful because the distance-to-stream functions may leave some unexpected regions as nodata if they do not drain to a stream. This may be confusing behavior, so this mask is intended to locate what drains to a stream and what does not. A pixel doesn't drain to a stream if it has a defined flow direction but undefined distance to stream. Args: - flow_dir_mfd_path (string): The path to an MFD flow direction raster. - This raster must have a nodata value defined. - dist_to_channel_mfd_path (string): The path to an MFD - distance-to-channel raster. This raster must have a nodata value - defined. + flow_dir_path (string): The path to a flow direction raster + (MFD or D8). This raster must have a nodata value defined. + dist_to_channel_path (string): The path to a distance-to-channel + raster. This raster must have a nodata value defined. target_mask_path (string): The path to where the mask raster should be written. Returns: ``None`` """ - flow_dir_mfd_nodata = pygeoprocessing.get_raster_info( - flow_dir_mfd_path)['nodata'][0] + flow_dir_nodata = pygeoprocessing.get_raster_info( + flow_dir_path)['nodata'][0] dist_to_channel_nodata = pygeoprocessing.get_raster_info( - dist_to_channel_mfd_path)['nodata'][0] + dist_to_channel_path)['nodata'][0] - def _what_drains_to_stream(flow_dir_mfd, dist_to_channel): + def _what_drains_to_stream(flow_dir, dist_to_channel): """Determine which pixels do and do not drain to a stream. Args: - flow_dir_mfd (numpy.array): A numpy array of MFD flow direction - values. + flow_dir (numpy.array): A numpy array of flow direction values. dist_to_channel (numpy.array): A numpy array of calculated distances to the nearest channel. Returns: A ``numpy.array`` of dtype ``numpy.uint8`` with pixels where: - * ``255`` where ``flow_dir_mfd`` is nodata (and thus + * ``255`` where ``flow_dir`` is nodata (and thus ``dist_to_channel`` is also nodata). - * ``0`` where ``flow_dir_mfd`` has data and ``dist_to_channel`` + * ``0`` where ``flow_dir`` has data and ``dist_to_channel`` does not - * ``1`` where ``flow_dir_mfd`` has data, and + * ``1`` where ``flow_dir`` has data, and ``dist_to_channel`` also has data. """ drains_to_stream = numpy.full( - flow_dir_mfd.shape, _BYTE_NODATA, dtype=numpy.uint8) + flow_dir.shape, _BYTE_NODATA, dtype=numpy.uint8) valid_flow_dir = ~pygeoprocessing.array_equals_nodata( - flow_dir_mfd, flow_dir_mfd_nodata) + flow_dir, flow_dir_nodata) valid_dist_to_channel = ( ~pygeoprocessing.array_equals_nodata( dist_to_channel, dist_to_channel_nodata) & @@ -1044,7 +1096,7 @@ def _what_drains_to_stream(flow_dir_mfd, dist_to_channel): return drains_to_stream pygeoprocessing.raster_calculator( - [(flow_dir_mfd_path, 1), (dist_to_channel_mfd_path, 1)], + [(flow_dir_path, 1), (dist_to_channel_path, 1)], _what_drains_to_stream, target_mask_path, gdal.GDT_Byte, _BYTE_NODATA) @@ -1338,7 +1390,7 @@ def _calculate_cp(lulc_to_cp, lulc_path, cp_factor_path): def _calculate_bar_factor( flow_direction_path, factor_path, flow_accumulation_path, - accumulation_path, out_bar_path): + accumulation_path, out_bar_path, algorithm): """Route user defined source across DEM. Used for calculating S and W bar in the SDR operation. @@ -1353,15 +1405,21 @@ def _calculate_bar_factor( out_bar_path (string): path to output raster that is the result of the factor accumulation raster divided by the flow accumulation raster. + algorithm (string): flow direction algorithm, 'D8' or 'MFD' Returns: None. """ - LOGGER.debug("doing flow accumulation mfd on %s", factor_path) + LOGGER.debug(f"doing flow accumulation on {factor_path}") + + if algorithm == 'D8': + flow_accum_func = pygeoprocessing.routing.flow_accumulation_d8 + else: # MFD + flow_accum_func = pygeoprocessing.routing.flow_accumulation_mfd # manually setting compression to DEFLATE because we got some LZW # errors when testing with large data. - pygeoprocessing.routing.flow_accumulation_mfd( + flow_accum_func( (flow_direction_path, 1), accumulation_path, weight_raster_path_band=(factor_path, 1), raster_driver_creation_tuple=('GTIFF', [ diff --git a/src/natcap/invest/sdr/sdr_core.pyx b/src/natcap/invest/sdr/sdr_core.pyx index 4cfb961a7..90df9407b 100644 --- a/src/natcap/invest/sdr/sdr_core.pyx +++ b/src/natcap/invest/sdr/sdr_core.pyx @@ -1,30 +1,19 @@ import logging -import os import pygeoprocessing cimport cython from osgeo import gdal -from libc.time cimport time as ctime -from libcpp.stack cimport stack -from ..managed_raster.managed_raster cimport ManagedRaster -from ..managed_raster.managed_raster cimport ManagedFlowDirRaster -from ..managed_raster.managed_raster cimport NeighborTuple -from ..managed_raster.managed_raster cimport is_close -from ..managed_raster.managed_raster cimport DownslopeNeighbors -from ..managed_raster.managed_raster cimport UpslopeNeighbors -from ..managed_raster.managed_raster cimport Pixel -from ..managed_raster.managed_raster cimport INFLOW_OFFSETS +from ..managed_raster.managed_raster cimport D8 +from ..managed_raster.managed_raster cimport MFD +from .sediment_deposition cimport run_sediment_deposition -cdef extern from "time.h" nogil: - ctypedef int time_t - time_t time(time_t*) LOGGER = logging.getLogger(__name__) def calculate_sediment_deposition( - mfd_flow_direction_path, e_prime_path, f_path, sdr_path, - target_sediment_deposition_path): + flow_direction_path, e_prime_path, f_path, sdr_path, + target_sediment_deposition_path, algorithm): """Calculate sediment deposition layer. This algorithm outputs both sediment deposition (t_i) and flux (f_i):: @@ -62,8 +51,8 @@ def calculate_sediment_deposition( will come from the SDR model and have nodata in the same places. Args: - mfd_flow_direction_path (string): a path to a raster with - pygeoprocessing.routing MFD flow direction values. + flow_direction_path (string): a path to a flow direction raster, + in either MFD or D8 format. Specify with the ``algorithm`` arg. e_prime_path (string): path to a raster that shows sources of sediment that wash off a pixel but do not reach the stream. f_path (string): path to a raster that shows the sediment flux @@ -71,6 +60,7 @@ def calculate_sediment_deposition( sdr_path (string): path to Sediment Delivery Ratio raster. target_sediment_deposition_path (string): path to created that shows where the E' sources end up across the landscape. + algorithm (string): MFD or D8 Returns: None. @@ -79,195 +69,19 @@ def calculate_sediment_deposition( LOGGER.info('Calculate sediment deposition') cdef float target_nodata = -1 pygeoprocessing.new_raster_from_base( - mfd_flow_direction_path, target_sediment_deposition_path, + flow_direction_path, target_sediment_deposition_path, gdal.GDT_Float32, [target_nodata]) pygeoprocessing.new_raster_from_base( - mfd_flow_direction_path, f_path, + flow_direction_path, f_path, gdal.GDT_Float32, [target_nodata]) - cdef ManagedFlowDirRaster mfd_flow_direction_raster = ManagedFlowDirRaster( - mfd_flow_direction_path.encode('utf-8'), 1, False) - cdef ManagedRaster e_prime_raster = ManagedRaster( - e_prime_path.encode('utf-8'), 1, False) - cdef ManagedRaster sdr_raster = ManagedRaster(sdr_path.encode('utf-8'), 1, False) - cdef ManagedRaster f_raster = ManagedRaster(f_path.encode('utf-8'), 1, True) - cdef ManagedRaster sediment_deposition_raster = ManagedRaster( - target_sediment_deposition_path.encode('utf-8'), 1, True) - - cdef long n_cols, n_rows - flow_dir_info = pygeoprocessing.get_raster_info(mfd_flow_direction_path) - n_cols, n_rows = flow_dir_info['raster_size'] - cdef int mfd_nodata = 0 - cdef stack[long] processing_stack - cdef float sdr_nodata = pygeoprocessing.get_raster_info( - sdr_path)['nodata'][0] - cdef float e_prime_nodata = pygeoprocessing.get_raster_info( - e_prime_path)['nodata'][0] - cdef long col_index, row_index, win_xsize, win_ysize, xoff, yoff - cdef long global_col, global_row, j, k - cdef int xs, ys - cdef long flat_index - cdef long seed_col = 0 - cdef long seed_row = 0 - cdef long neighbor_row, neighbor_col, ds_neighbor_row, ds_neighbor_col - cdef int flow_val, neighbor_flow_val, ds_neighbor_flow_val - cdef int flow_weight, neighbor_flow_weight - cdef float flow_sum, neighbor_flow_sum - cdef float downslope_sdr_weighted_sum, sdr_i, sdr_j - cdef float p_j, p_val - cdef unsigned long n_pixels_processed = 0 - cdef time_t last_log_time = ctime(NULL) - cdef float f_j_weighted_sum - cdef NeighborTuple neighbor - cdef UpslopeNeighbors up_neighbors - - for offset_dict in pygeoprocessing.iterblocks( - (mfd_flow_direction_path, 1), offset_only=True, largest_block=0): - # use cython variables to avoid python overhead of dict values - win_xsize = offset_dict['win_xsize'] - win_ysize = offset_dict['win_ysize'] - xoff = offset_dict['xoff'] - yoff = offset_dict['yoff'] - if ctime(NULL) - last_log_time > 5.0: - last_log_time = ctime(NULL) - LOGGER.info('Sediment deposition %.2f%% complete', 100 * ( - n_pixels_processed / float(n_cols * n_rows))) - - for row_index in range(win_ysize): - ys = yoff + row_index - for col_index in range(win_xsize): - xs = xoff + col_index - - - if mfd_flow_direction_raster.get(xs, ys) == mfd_nodata: - continue - - # if this can be a seed pixel and hasn't already been - # calculated, put it on the stack - if (mfd_flow_direction_raster.is_local_high_point(xs, ys) and - is_close(sediment_deposition_raster.get(xs, ys), target_nodata)): - processing_stack.push(ys * n_cols + xs) - - while processing_stack.size() > 0: - # loop invariant: cell has all upslope neighbors - # processed. this is true for seed pixels because they - # have no upslope neighbors. - flat_index = processing_stack.top() - processing_stack.pop() - global_row = flat_index // n_cols - global_col = flat_index % n_cols - - # (sum over j ∈ J of f_j * p(i,j) in the equation for t_i) - # calculate the upslope f_j contribution to this pixel, - # the weighted sum of flux flowing onto this pixel from - # all neighbors - f_j_weighted_sum = 0 - up_neighbors = UpslopeNeighbors( - Pixel(mfd_flow_direction_raster, global_col, global_row)) - for neighbor in up_neighbors: - f_j = f_raster.get(neighbor.x, neighbor.y) - if is_close(f_j, target_nodata): - continue - # add the neighbor's flux value, weighted by the - # flow proportion - f_j_weighted_sum += neighbor.flow_proportion * f_j - - # calculate sum of SDR values of immediate downslope - # neighbors, weighted by proportion of flow into each - # neighbor - # (sum over k ∈ K of SDR_k * p(i,k) in the equation above) - downslope_sdr_weighted_sum = 0 - dn_neighbors = DownslopeNeighbors( - Pixel(mfd_flow_direction_raster, global_col, global_row)) - - flow_dir_sum = 0 - for neighbor in dn_neighbors: - flow_dir_sum += neighbor.flow_proportion - sdr_j = sdr_raster.get(neighbor.x, neighbor.y) - if is_close(sdr_j, sdr_nodata): - continue - if sdr_j == 0: - # this means it's a stream, for SDR deposition - # purposes, we set sdr to 1 to indicate this - # is the last step on which to retain sediment - sdr_j = 1 - - downslope_sdr_weighted_sum += ( - sdr_j * neighbor.flow_proportion) - # check if we can add neighbor j to the stack yet - # - # if there is a downslope neighbor it - # couldn't have been pushed on the processing - # stack yet, because the upslope was just - # completed - upslope_neighbors_processed = 1 - # iterate over each neighbor-of-neighbor - up_neighbors = UpslopeNeighbors( - Pixel(mfd_flow_direction_raster, neighbor.x, neighbor.y)) - for neighbor_of_neighbor in up_neighbors: - if INFLOW_OFFSETS[neighbor_of_neighbor.direction] == neighbor.direction: - continue - if is_close(sediment_deposition_raster.get( - neighbor_of_neighbor.x, neighbor_of_neighbor.y - ), target_nodata): - upslope_neighbors_processed = 0 - break - # if all upslope neighbors of neighbor j are - # processed, we can push j onto the stack. - if upslope_neighbors_processed: - processing_stack.push( - neighbor.y * n_cols + neighbor.x) - - # nodata pixels should propagate to the results - sdr_i = sdr_raster.get(global_col, global_row) - if is_close(sdr_i, sdr_nodata): - continue - e_prime_i = e_prime_raster.get(global_col, global_row) - if is_close(e_prime_i, e_prime_nodata): - continue - - if flow_dir_sum: - downslope_sdr_weighted_sum /= flow_dir_sum - - # This condition reflects property A in the user's guide. - if downslope_sdr_weighted_sum < sdr_i: - # i think this happens because of our low resolution - # flow direction, it's okay to zero out. - downslope_sdr_weighted_sum = sdr_i - - # these correspond to the full equations for - # dr_i, t_i, and f_i given in the docstring - if sdr_i == 1: - # This reflects property B in the user's guide and is - # an edge case to avoid division-by-zero. - dt_i = 1 - else: - dt_i = (downslope_sdr_weighted_sum - sdr_i) / (1 - sdr_i) - - # Lisa's modified equations - t_i = dt_i * f_j_weighted_sum # deposition, a.k.a trapped sediment - f_i = (1 - dt_i) * f_j_weighted_sum + e_prime_i # flux - - # On large flow paths, it's possible for dt_i, f_i and t_i - # to have very small negative values that are numerically - # equivalent to 0. These negative values were raising - # questions on the forums and it's easier to clamp the - # values here than to explain IEEE 754. - if dt_i < 0: - dt_i = 0 - if t_i < 0: - t_i = 0 - if f_i < 0: - f_i = 0 - - sediment_deposition_raster.set(global_col, global_row, t_i) - f_raster.set(global_col, global_row, f_i) - n_pixels_processed += win_xsize * win_ysize - - sediment_deposition_raster.close() - mfd_flow_direction_raster.close() - e_prime_raster.close() - sdr_raster.close() - f_raster.close() - - LOGGER.info('Sediment deposition 100% complete') + if algorithm == 'D8': + run_sediment_deposition[D8]( + flow_direction_path.encode('utf-8'), e_prime_path.encode('utf-8'), + f_path.encode('utf-8'), sdr_path.encode('utf-8'), + target_sediment_deposition_path.encode('utf-8')) + else: + run_sediment_deposition[MFD]( + flow_direction_path.encode('utf-8'), e_prime_path.encode('utf-8'), + f_path.encode('utf-8'), sdr_path.encode('utf-8'), + target_sediment_deposition_path.encode('utf-8')) diff --git a/src/natcap/invest/sdr/sediment_deposition.h b/src/natcap/invest/sdr/sediment_deposition.h new file mode 100644 index 000000000..2924adb87 --- /dev/null +++ b/src/natcap/invest/sdr/sediment_deposition.h @@ -0,0 +1,270 @@ +#include "ManagedRaster.h" + + +// Calculate sediment deposition layer. +// +// This algorithm outputs both sediment deposition (t_i) and flux (f_i):: +// +// t_i = dt_i * (sum over j ∈ J of f_j * p(j,i)) +// +// f_i = (1 - dt_i) * (sum over j ∈ J of f_j * p(j,i)) + E'_i +// +// +// (sum over k ∈ K of SDR_k * p(i,k)) - SDR_i +// dt_i = -------------------------------------------- +// (1 - SDR_i) +// +// where: +// +// - ``p(i,j)`` is the proportion of flow from pixel ``i`` into pixel ``j`` +// - ``J`` is the set of pixels that are immediate upslope neighbors of +// pixel ``i`` +// - ``K`` is the set of pixels that are immediate downslope neighbors of +// pixel ``i`` +// - ``E'`` is ``USLE * (1 - SDR)``, the amount of sediment loss from pixel +// ``i`` that doesn't reach a stream (``e_prime_path``) +// - ``SDR`` is the sediment delivery ratio (``sdr_path``) +// +// ``f_i`` is recursively defined in terms of ``i``'s upslope neighbors. +// The algorithm begins from seed pixels that are local high points and so +// have no upslope neighbors. It works downslope from each seed pixel, +// only adding a pixel to the stack when all its upslope neighbors are +// already calculated. +// +// Note that this function is designed to be used in the context of the SDR +// model. Because the algorithm is recursive upslope and downslope of each +// pixel, nodata values in the SDR input would propagate along the flow path. +// This case is not handled because we assume the SDR and flow dir inputs +// will come from the SDR model and have nodata in the same places. +// +// Args: +// flow_direction_path: a path to a flow direction raster, +// in either MFD or D8 format. Specify with the ``algorithm`` arg. +// e_prime_path: path to a raster that shows sources of +// sediment that wash off a pixel but do not reach the stream. +// f_path: path to a raster that shows the sediment flux +// on a pixel for sediment that does not reach the stream. +// sdr_path: path to Sediment Delivery Ratio raster. +// target_sediment_deposition_path: path to created that +// shows where the E' sources end up across the landscape. +template +void run_sediment_deposition( + char* flow_direction_path, + char* e_prime_path, + char* f_path, + char* sdr_path, + char* sediment_deposition_path) { + + ManagedFlowDirRaster flow_dir_raster = ManagedFlowDirRaster( + flow_direction_path, 1, false); + + ManagedRaster e_prime_raster = ManagedRaster(e_prime_path, 1, false); + ManagedRaster sdr_raster = ManagedRaster(sdr_path, 1, false); + ManagedRaster f_raster = ManagedRaster(f_path, 1, true); + ManagedRaster sediment_deposition_raster = ManagedRaster( + sediment_deposition_path, 1, true); + + stack processing_stack; + float target_nodata = -1; + long win_xsize, win_ysize, xoff, yoff; + long global_col, global_row; + int xs, ys; + long flat_index; + double downslope_sdr_weighted_sum; + double sdr_i, e_prime_i, sdr_j, f_j; + long flow_dir_sum; + + // unsigned long n_pixels_processed = 0; + bool upslope_neighbors_processed; + // time_t last_log_time = ctime(NULL) + double f_j_weighted_sum; + NeighborTuple neighbor; + NeighborTuple neighbor_of_neighbor; + double dr_i, t_i, f_i; + + UpslopeNeighbors up_neighbors; + DownslopeNeighbors dn_neighbors; + + // efficient way to calculate ceiling division: + // a divided by b rounded up = (a + (b - 1)) / b + // note that / represents integer floor division + // https://stackoverflow.com/a/62032709/14451410 + int n_col_blocks = (flow_dir_raster.raster_x_size + (flow_dir_raster.block_xsize - 1)) / flow_dir_raster.block_xsize; + int n_row_blocks = (flow_dir_raster.raster_y_size + (flow_dir_raster.block_ysize - 1)) / flow_dir_raster.block_ysize; + + for (int row_block_index = 0; row_block_index < n_row_blocks; row_block_index++) { + yoff = row_block_index * flow_dir_raster.block_ysize; + win_ysize = flow_dir_raster.raster_y_size - yoff; + if (win_ysize > flow_dir_raster.block_ysize) { + win_ysize = flow_dir_raster.block_ysize; + } + for (int col_block_index = 0; col_block_index < n_col_blocks; col_block_index++) { + xoff = col_block_index * flow_dir_raster.block_xsize; + win_xsize = flow_dir_raster.raster_x_size - xoff; + if (win_xsize > flow_dir_raster.block_xsize) { + win_xsize = flow_dir_raster.block_xsize; + } + + // if ctime(NULL) - last_log_time > 5.0: + // last_log_time = ctime(NULL) + // LOGGER.info('Sediment deposition %.2f%% complete', 100 * ( + // n_pixels_processed / float(flow_dir_raster.raster_x_size * flow_dir_raster.raster_y_size))) + + for (int row_index = 0; row_index < win_ysize; row_index++) { + ys = yoff + row_index; + for (int col_index = 0; col_index < win_xsize; col_index++) { + xs = xoff + col_index; + + if (flow_dir_raster.get(xs, ys) == flow_dir_raster.nodata) { + continue; + } + + // if this can be a seed pixel and hasn't already been + // calculated, put it on the stack + if (flow_dir_raster.is_local_high_point(xs, ys) and + is_close(sediment_deposition_raster.get(xs, ys), target_nodata)) { + processing_stack.push(ys * flow_dir_raster.raster_x_size + xs); + } + + while (processing_stack.size() > 0) { + // # loop invariant: cell has all upslope neighbors + // # processed. this is true for seed pixels because they + // # have no upslope neighbors. + flat_index = processing_stack.top(); + processing_stack.pop(); + global_row = flat_index / flow_dir_raster.raster_x_size; + global_col = flat_index % flow_dir_raster.raster_x_size; + + // # (sum over j ∈ J of f_j * p(i,j) in the equation for t_i) + // # calculate the upslope f_j contribution to this pixel, + // # the weighted sum of flux flowing onto this pixel from + // # all neighbors + f_j_weighted_sum = 0; + up_neighbors = UpslopeNeighbors( + Pixel(flow_dir_raster, global_col, global_row)); + for (auto neighbor: up_neighbors) { + f_j = f_raster.get(neighbor.x, neighbor.y); + if (is_close(f_j, target_nodata)) { + continue; + } + // add the neighbor's flux value, weighted by the + // flow proportion + f_j_weighted_sum += neighbor.flow_proportion * f_j; + } + + // # calculate sum of SDR values of immediate downslope + // # neighbors, weighted by proportion of flow into each + // # neighbor + // # (sum over k ∈ K of SDR_k * p(i,k) in the equation above) + downslope_sdr_weighted_sum = 0; + dn_neighbors = DownslopeNeighbors( + Pixel(flow_dir_raster, global_col, global_row)); + flow_dir_sum = 0; + for (auto neighbor: dn_neighbors) { + flow_dir_sum += static_cast(neighbor.flow_proportion); + sdr_j = sdr_raster.get(neighbor.x, neighbor.y); + if (is_close(sdr_j, sdr_raster.nodata)) { + continue; + } + if (sdr_j == 0) { + // # this means it's a stream, for SDR deposition + // # purposes, we set sdr to 1 to indicate this + // # is the last step on which to retain sediment + sdr_j = 1; + } + + downslope_sdr_weighted_sum += ( + sdr_j * neighbor.flow_proportion); + // # check if we can add neighbor j to the stack yet + // # + // # if there is a downslope neighbor it + // # couldn't have been pushed on the processing + // # stack yet, because the upslope was just + // # completed + upslope_neighbors_processed = true; + // # iterate over each neighbor-of-neighbor + up_neighbors = UpslopeNeighbors( + Pixel(flow_dir_raster, neighbor.x, neighbor.y)); + for (auto neighbor_of_neighbor: up_neighbors) { + if (INFLOW_OFFSETS[neighbor_of_neighbor.direction] == neighbor.direction) { + continue; + } + if (is_close(sediment_deposition_raster.get( + neighbor_of_neighbor.x, neighbor_of_neighbor.y + ), target_nodata)) { + upslope_neighbors_processed = false; + break; + } + } + // # if all upslope neighbors of neighbor j are + // # processed, we can push j onto the stack. + if (upslope_neighbors_processed) { + processing_stack.push( + neighbor.y * flow_dir_raster.raster_x_size + neighbor.x); + } + } + + // # nodata pixels should propagate to the results + sdr_i = sdr_raster.get(global_col, global_row); + if (is_close(sdr_i, sdr_raster.nodata)) { + continue; + } + e_prime_i = e_prime_raster.get(global_col, global_row); + if (is_close(e_prime_i, e_prime_raster.nodata)) { + continue; + } + + if (flow_dir_sum) { + downslope_sdr_weighted_sum /= flow_dir_sum; + } + + // # This condition reflects property A in the user's guide. + if (downslope_sdr_weighted_sum < sdr_i) { + // # i think this happens because of our low resolution + // # flow direction, it's okay to zero out. + downslope_sdr_weighted_sum = sdr_i; + } + + // # these correspond to the full equations for + // # dr_i, t_i, and f_i given in the docstring + if (sdr_i == 1) { + // # This reflects property B in the user's guide and is + // # an edge case to avoid division-by-zero. + dr_i = 1; + } else { + dr_i = (downslope_sdr_weighted_sum - sdr_i) / (1 - sdr_i); + } + + // # Lisa's modified equations + t_i = dr_i * f_j_weighted_sum; // deposition, a.k.a trapped sediment + f_i = (1 - dr_i) * f_j_weighted_sum + e_prime_i; // flux + + // # On large flow paths, it's possible for dr_i, f_i and t_i + // # to have very small negative values that are numerically + // # equivalent to 0. These negative values were raising + // # questions on the forums and it's easier to clamp the + // # values here than to explain IEEE 754. + if (dr_i < 0) { + dr_i = 0; + } + if (t_i < 0) { + t_i = 0; + } + if (f_i < 0) { + f_i = 0; + } + sediment_deposition_raster.set(global_col, global_row, t_i); + f_raster.set(global_col, global_row, f_i); + } + } + } + // n_pixels_processed += win_xsize * win_ysize; + } + } + sediment_deposition_raster.close(); + flow_dir_raster.close(); + e_prime_raster.close(); + sdr_raster.close(); + f_raster.close(); +// // LOGGER.info('Sediment deposition 100% complete') +} diff --git a/src/natcap/invest/sdr/sediment_deposition.pxd b/src/natcap/invest/sdr/sediment_deposition.pxd new file mode 100644 index 000000000..41143986c --- /dev/null +++ b/src/natcap/invest/sdr/sediment_deposition.pxd @@ -0,0 +1,7 @@ +cdef extern from "sediment_deposition.h": + void run_sediment_deposition[T]( + char*, + char*, + char*, + char*, + char*) diff --git a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py index e5c2ba57e..ed5eb2981 100644 --- a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py +++ b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py @@ -283,6 +283,21 @@ "Table of alpha values for each month. " "Required if Use Monthly Alpha Table is selected."), "name": gettext("monthly alpha table") + }, + "algorithm": { + "type": "option_string", + "options": { + "D8": { + "display_name": gettext("D8"), + "description": "D8 flow direction" + }, + "MFD": { + "display_name": gettext("MFD"), + "description": "Multiple flow direction" + } + }, + "about": gettext("Flow direction algorithm to use."), + "name": gettext("flow direction algorithm") } }, "outputs": { @@ -398,10 +413,10 @@ "units": u.millimeter }} }, - "flow_dir_mfd.tif": { + "flow_dir.tif": { "about": gettext( - "Map of multiple flow direction. Values are encoded in " - "a binary format and should not be used directly."), + "Map of flow direction, in either D8 or MFD format " + "according to the option selected."), "bands": {1: {"type": "integer"}} }, "qf_[MONTH].tif": { @@ -497,7 +512,7 @@ _INTERMEDIATE_BASE_FILES = { 'aet_path': 'aet.tif', 'aetm_path_list': ['aetm_%d.tif' % (x+1) for x in range(N_MONTHS)], - 'flow_dir_mfd_path': 'flow_dir_mfd.tif', + 'flow_dir_path': 'flow_dir.tif', 'qfm_path_list': ['qf_%d.tif' % (x+1) for x in range(N_MONTHS)], 'stream_path': 'stream.tif', 'si_path': 'Si.tif', @@ -701,35 +716,67 @@ def execute(args): dependent_task_list=[align_task], task_name='fill dem pits') - flow_dir_task = task_graph.add_task( - func=pygeoprocessing.routing.flow_dir_mfd, - args=( - (file_registry['dem_pit_filled_path'], 1), - file_registry['flow_dir_mfd_path']), - kwargs={'working_dir': intermediate_output_dir}, - target_path_list=[file_registry['flow_dir_mfd_path']], - dependent_task_list=[fill_pit_task], - task_name='flow dir mfd') + if args['algorithm'] == 'MFD': + flow_dir_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_dir_mfd, + args=( + (file_registry['dem_pit_filled_path'], 1), + file_registry['flow_dir_path']), + kwargs={'working_dir': intermediate_output_dir}, + target_path_list=[file_registry['flow_dir_path']], + dependent_task_list=[fill_pit_task], + task_name='flow direction - MFD') + + flow_accum_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_accumulation_mfd, + args=( + (file_registry['flow_dir_path'], 1), + file_registry['flow_accum_path']), + target_path_list=[file_registry['flow_accum_path']], + dependent_task_list=[flow_dir_task], + task_name='flow accumulation - MFD') + + stream_threshold_task = task_graph.add_task( + func=pygeoprocessing.routing.extract_streams_mfd, + args=( + (file_registry['flow_accum_path'], 1), + (file_registry['flow_dir_path'], 1), + threshold_flow_accumulation, + file_registry['stream_path']), + target_path_list=[file_registry['stream_path']], + dependent_task_list=[flow_accum_task], + task_name='stream threshold - MFD') + else: # D8 + flow_dir_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_dir_d8, + args=( + (file_registry['dem_pit_filled_path'], 1), + file_registry['flow_dir_path']), + kwargs={'working_dir': intermediate_output_dir}, + target_path_list=[file_registry['flow_dir_path']], + dependent_task_list=[fill_pit_task], + task_name='flow direction - D8') + + flow_accum_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_accumulation_d8, + args=( + (file_registry['flow_dir_path'], 1), + file_registry['flow_accum_path']), + target_path_list=[file_registry['flow_accum_path']], + dependent_task_list=[flow_dir_task], + task_name='flow accumulation - D8') + + stream_threshold_task = task_graph.add_task( + func=pygeoprocessing.routing.extract_streams_d8, + kwargs=dict( + flow_accum_raster_path_band=(file_registry['flow_accum_path'], 1), + flow_threshold=threshold_flow_accumulation, + target_stream_raster_path=file_registry['stream_path']), + target_path_list=[file_registry['stream_path']], + dependent_task_list=[flow_accum_task], + task_name='stream threshold - D8') + - flow_accum_task = task_graph.add_task( - func=pygeoprocessing.routing.flow_accumulation_mfd, - args=( - (file_registry['flow_dir_mfd_path'], 1), - file_registry['flow_accum_path']), - target_path_list=[file_registry['flow_accum_path']], - dependent_task_list=[flow_dir_task], - task_name='flow accum task') - - stream_threshold_task = task_graph.add_task( - func=pygeoprocessing.routing.extract_streams_mfd, - args=( - (file_registry['flow_accum_path'], 1), - (file_registry['flow_dir_mfd_path'], 1), - threshold_flow_accumulation, - file_registry['stream_path']), - target_path_list=[file_registry['stream_path']], - dependent_task_list=[flow_accum_task], - task_name='stream threshold') LOGGER.info('quick flow') if args['user_defined_local_recharge']: @@ -862,7 +909,7 @@ def execute(args): file_registry['precip_path_aligned_list'], file_registry['et0_path_aligned_list'], file_registry['qfm_path_list'], - file_registry['flow_dir_mfd_path'], + file_registry['flow_dir_path'], file_registry['kc_path_list'], alpha_month_map, beta_i, gamma, file_registry['stream_path'], @@ -870,7 +917,8 @@ def execute(args): file_registry['l_avail_path'], file_registry['l_sum_avail_path'], file_registry['aet_path'], - file_registry['annual_precip_path']), + file_registry['annual_precip_path'], + args['algorithm']), target_path_list=[ file_registry['l_path'], file_registry['l_avail_path'], @@ -907,16 +955,28 @@ def execute(args): task_name='aggregate recharge') LOGGER.info('calculate L_sum') # Eq. [12] - l_sum_task = task_graph.add_task( - func=pygeoprocessing.routing.flow_accumulation_mfd, - args=( - (file_registry['flow_dir_mfd_path'], 1), - file_registry['l_sum_path']), - kwargs={'weight_raster_path_band': (file_registry['l_path'], 1)}, - target_path_list=[file_registry['l_sum_path']], - dependent_task_list=vri_dependent_task_list + [ - fill_pit_task, flow_dir_task, stream_threshold_task], - task_name='calculate l sum') + if args['algorithm'] == 'MFD': + l_sum_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_accumulation_mfd, + args=( + (file_registry['flow_dir_path'], 1), + file_registry['l_sum_path']), + kwargs={'weight_raster_path_band': (file_registry['l_path'], 1)}, + target_path_list=[file_registry['l_sum_path']], + dependent_task_list=vri_dependent_task_list + [ + fill_pit_task, flow_dir_task, stream_threshold_task], + task_name='calculate l sum - MFD') + else: # D8 + l_sum_task = task_graph.add_task( + func=pygeoprocessing.routing.flow_accumulation_d8, + args=( + (file_registry['flow_dir_path'], 1), + file_registry['l_sum_path']), + kwargs={'weight_raster_path_band': (file_registry['l_path'], 1)}, + target_path_list=[file_registry['l_sum_path']], + dependent_task_list=vri_dependent_task_list + [ + fill_pit_task, flow_dir_task, stream_threshold_task], + task_name='calculate l sum - D8') if args['user_defined_local_recharge']: b_sum_dependent_task_list = [l_avail_task] @@ -926,14 +986,14 @@ def execute(args): b_sum_task = task_graph.add_task( func=seasonal_water_yield_core.route_baseflow_sum, args=( - file_registry['flow_dir_mfd_path'], + file_registry['flow_dir_path'], file_registry['l_path'], file_registry['l_avail_path'], file_registry['l_sum_path'], file_registry['stream_path'], file_registry['b_path'], - file_registry['b_sum_path']), - + file_registry['b_sum_path'], + args['algorithm']), target_path_list=[ file_registry['b_sum_path'], file_registry['b_path']], dependent_task_list=b_sum_dependent_task_list + [l_sum_task], diff --git a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield_core.pyx b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield_core.pyx index 3e37c25e5..253b2e723 100644 --- a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield_core.pyx +++ b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield_core.pyx @@ -12,33 +12,18 @@ from osgeo import gdal from osgeo import ogr from osgeo import osr -from libcpp.pair cimport pair -from libcpp.stack cimport stack -from libcpp.queue cimport queue from libcpp.vector cimport vector -from libc.time cimport time as ctime -from ..managed_raster.managed_raster cimport ManagedRaster -from ..managed_raster.managed_raster cimport ManagedFlowDirRaster -from ..managed_raster.managed_raster cimport Pixel -from ..managed_raster.managed_raster cimport DownslopeNeighbors -from ..managed_raster.managed_raster cimport UpslopeNeighbors -from ..managed_raster.managed_raster cimport UpslopeNeighborsNoDivide -from ..managed_raster.managed_raster cimport NeighborTuple -from ..managed_raster.managed_raster cimport is_close -cdef extern from "time.h" nogil: - ctypedef int time_t - time_t time(time_t*) +from ..managed_raster.managed_raster cimport D8, MFD +from .swy cimport run_route_baseflow_sum, run_calculate_local_recharge LOGGER = logging.getLogger(__name__) -cdef int N_MONTHS = 12 - cpdef calculate_local_recharge( precip_path_list, et0_path_list, qf_m_path_list, flow_dir_mfd_path, kc_path_list, alpha_month_map, float beta_i, float gamma, stream_path, target_li_path, target_li_avail_path, target_l_sum_avail_path, - target_aet_path, target_pi_path): + target_aet_path, target_pi_path, algorithm): """ Calculate the rasters defined by equations [3]-[7]. @@ -78,255 +63,67 @@ cpdef calculate_local_recharge( None. """ - cdef int i_n, flow_dir_nodata, flow_dir_mfd - cdef int peak_pixel - cdef long xs, ys, xs_root, ys_root, xoff, yoff - cdef int flow_dir_s - cdef long xi, yi, xj, yj - cdef int flow_dir_j - cdef int n_dir - cdef long raster_x_size, raster_y_size, win_xsize, win_ysize - cdef double pet_m, p_m, qf_m, et0_m, aet_i, p_i, qf_i, l_i, l_avail_i - cdef float qf_nodata, kc_nodata - - cdef float mfd_direction_array[8] - - cdef queue[pair[long, long]] work_queue - cdef ManagedRaster et0_m_raster, qf_m_raster, kc_m_raster - - cdef numpy.ndarray[numpy.npy_float32, ndim=1] alpha_month_array = ( - numpy.array( - [x[1] for x in sorted(alpha_month_map.items())], - dtype=numpy.float32)) - - # used for time-delayed logging - cdef time_t last_log_time - last_log_time = ctime(NULL) - - # we know the PyGeoprocessing MFD raster flow dir type is a 32 bit int. - flow_dir_raster_info = pygeoprocessing.get_raster_info(flow_dir_mfd_path) - flow_dir_nodata = flow_dir_raster_info['nodata'][0] - raster_x_size, raster_y_size = flow_dir_raster_info['raster_size'] - - cdef ManagedFlowDirRaster flow_raster = ManagedFlowDirRaster( - flow_dir_mfd_path.encode('utf-8'), 1, 0) - cdef NeighborTuple neighbor - - # make sure that user input nodata values are defined - # set to -1 if not defined - # precipitation and evapotranspiration data should - # always be non-negative - cdef vector[ManagedRaster] et0_m_rasters - et0_m_nodata_list = [] - for et0_path in et0_path_list: - et0_m_rasters.push_back(ManagedRaster(et0_path.encode('utf-8'), 1, 0)) - nodata = pygeoprocessing.get_raster_info(et0_path)['nodata'][0] - if nodata is None: - nodata = -1 - et0_m_nodata_list.append(nodata) - - cdef vector[ManagedRaster] precip_m_rasters - precip_m_nodata_list = [] - for precip_m_path in precip_path_list: - precip_m_rasters.push_back(ManagedRaster(precip_m_path.encode('utf-8'), 1, 0)) - nodata = pygeoprocessing.get_raster_info(precip_m_path)['nodata'][0] - if nodata is None: - nodata = -1 - precip_m_nodata_list.append(nodata) - - cdef vector[ManagedRaster] qf_m_rasters - qf_m_nodata_list = [] - for qf_m_path in qf_m_path_list: - qf_m_rasters.push_back(ManagedRaster(qf_m_path.encode('utf-8'), 1, 0)) - qf_m_nodata_list.append( - pygeoprocessing.get_raster_info(qf_m_path)['nodata'][0]) - - cdef vector[ManagedRaster] kc_m_rasters - kc_m_nodata_list = [] - for kc_m_path in kc_path_list: - kc_m_rasters.push_back(ManagedRaster(kc_m_path.encode('utf-8'), 1, 0)) - kc_m_nodata_list.append( - pygeoprocessing.get_raster_info(kc_m_path)['nodata'][0]) + cdef vector[float] alpha_values + cdef vector[char*] et0_paths + cdef vector[char*] precip_paths + cdef vector[char*] qf_paths + cdef vector[char*] kc_paths + encoded_et0_paths = [p.encode('utf-8') for p in et0_path_list] + encoded_precip_paths = [p.encode('utf-8') for p in precip_path_list] + encoded_qf_paths = [p.encode('utf-8') for p in qf_m_path_list] + encoded_kc_paths = [p.encode('utf-8') for p in kc_path_list] + for i in range(12): + et0_paths.push_back(encoded_et0_paths[i]) + precip_paths.push_back(encoded_precip_paths[i]) + qf_paths.push_back(encoded_qf_paths[i]) + kc_paths.push_back(encoded_kc_paths[i]) + alpha_values.push_back(alpha_month_map[i + 1]) target_nodata = -1e32 pygeoprocessing.new_raster_from_base( flow_dir_mfd_path, target_li_path, gdal.GDT_Float32, [target_nodata], fill_value_list=[target_nodata]) - cdef ManagedRaster target_li_raster = ManagedRaster( - target_li_path.encode('utf-8'), 1, 1) - pygeoprocessing.new_raster_from_base( flow_dir_mfd_path, target_li_avail_path, gdal.GDT_Float32, [target_nodata], fill_value_list=[target_nodata]) - cdef ManagedRaster target_li_avail_raster = ManagedRaster( - target_li_avail_path.encode('utf-8'), 1, 1) - pygeoprocessing.new_raster_from_base( flow_dir_mfd_path, target_l_sum_avail_path, gdal.GDT_Float32, [target_nodata], fill_value_list=[target_nodata]) - cdef ManagedRaster target_l_sum_avail_raster = ManagedRaster( - target_l_sum_avail_path.encode('utf-8'), 1, 1) - pygeoprocessing.new_raster_from_base( flow_dir_mfd_path, target_aet_path, gdal.GDT_Float32, [target_nodata], fill_value_list=[target_nodata]) - cdef ManagedRaster target_aet_raster = ManagedRaster( - target_aet_path.encode('utf-8'), 1, 1) - pygeoprocessing.new_raster_from_base( flow_dir_mfd_path, target_pi_path, gdal.GDT_Float32, [target_nodata], fill_value_list=[target_nodata]) - cdef ManagedRaster target_pi_raster = ManagedRaster( - target_pi_path.encode('utf-8'), 1, 1) - - for offset_dict in pygeoprocessing.iterblocks( - (flow_dir_mfd_path, 1), offset_only=True, largest_block=0): - # use cython variables to avoid python overhead of dict values - win_xsize = offset_dict['win_xsize'] - win_ysize = offset_dict['win_ysize'] - xoff = offset_dict['xoff'] - yoff = offset_dict['yoff'] - if ctime(NULL) - last_log_time > 5.0: - last_log_time = ctime(NULL) - current_pixel = xoff + yoff * raster_x_size - LOGGER.info( - 'peak point detection %.2f%% complete', - 100.0 * current_pixel / ( - raster_x_size * raster_y_size)) - - # search block for a peak pixel where no other pixel drains to it. - for ys in xrange(win_ysize): - ys_root = yoff + ys - for xs in xrange(win_xsize): - xs_root = xoff + xs - - if flow_raster.is_local_high_point(xs_root, ys_root): - work_queue.push(pair[long, long](xs_root, ys_root)) - - while work_queue.size() > 0: - xi = work_queue.front().first - yi = work_queue.front().second - work_queue.pop() - - l_sum_avail_i = target_l_sum_avail_raster.get(xi, yi) - if not is_close(l_sum_avail_i, target_nodata): - # already defined - continue - - # Equation 7, calculate L_sum_avail_i if possible, skip - # otherwise - upslope_defined = 1 - # initialize to 0 so we indicate we haven't tracked any - # mfd values yet - l_sum_avail_i = 0.0 - mfd_dir_sum = 0 - up_neighbors = UpslopeNeighborsNoDivide(Pixel(flow_raster, xi, yi)) - for neighbor in up_neighbors: - # pixel flows inward, check upslope - l_sum_avail_j = target_l_sum_avail_raster.get( - neighbor.x, neighbor.y) - if is_close(l_sum_avail_j, target_nodata): - upslope_defined = 0 - break - l_avail_j = target_li_avail_raster.get( - neighbor.x, neighbor.y) - # A step of Equation 7 - l_sum_avail_i += ( - l_sum_avail_j + l_avail_j) * neighbor.flow_proportion - mfd_dir_sum += neighbor.flow_proportion - # calculate l_sum_avail_i by summing all the valid - # directions then normalizing by the sum of the mfd - # direction weights (Equation 8) - if upslope_defined: - # Equation 7 - if mfd_dir_sum > 0: - l_sum_avail_i /= mfd_dir_sum - target_l_sum_avail_raster.set(xi, yi, l_sum_avail_i) - else: - # if not defined, we'll get it on another pass - continue - - aet_i = 0 - p_i = 0 - qf_i = 0 - - for m_index in range(12): - precip_m_raster = ( - precip_m_rasters[m_index]) - qf_m_raster = ( - qf_m_rasters[m_index]) - et0_m_raster = ( - et0_m_rasters[m_index]) - kc_m_raster = ( - kc_m_rasters[m_index]) - - et0_nodata = et0_m_nodata_list[m_index] - precip_nodata = precip_m_nodata_list[m_index] - qf_nodata = qf_m_nodata_list[m_index] - kc_nodata = kc_m_nodata_list[m_index] - - p_m = precip_m_raster.get(xi, yi) - if not is_close(p_m, precip_nodata): - p_i += p_m - else: - p_m = 0 - - qf_m = qf_m_raster.get(xi, yi) - if not is_close(qf_m, qf_nodata): - qf_i += qf_m - else: - qf_m = 0 - - kc_m = kc_m_raster.get(xi, yi) - pet_m = 0 - et0_m = et0_m_raster.get(xi, yi) - if not ( - is_close(kc_m, kc_nodata) or - is_close(et0_m, et0_nodata)): - # Equation 6 - pet_m = kc_m * et0_m - - # Equation 4/5 - aet_i += min( - pet_m, - p_m - qf_m + - alpha_month_array[m_index]*beta_i*l_sum_avail_i) - - l_i = (p_i - qf_i - aet_i) - l_avail_i = min(gamma * l_i, l_i) - - target_pi_raster.set(xi, yi, p_i) - target_aet_raster.set(xi, yi, aet_i) - target_li_raster.set(xi, yi, l_i) - target_li_avail_raster.set(xi, yi, l_avail_i) - - dn_neighbors = DownslopeNeighbors(Pixel(flow_raster, xi, yi)) - for neighbor in dn_neighbors: - work_queue.push(pair[long, long](neighbor.x, neighbor.y)) - - flow_raster.close() - target_li_raster.close() - target_li_avail_raster.close() - target_l_sum_avail_raster.close() - target_aet_raster.close() - target_pi_raster.close() - for raster in et0_m_rasters: - raster.close() - for raster in precip_m_rasters: - raster.close() - for raster in qf_m_rasters: - raster.close() - for raster in kc_m_rasters: - raster.close() + args = [ + precip_paths, + et0_paths, + qf_paths, + flow_dir_mfd_path.encode('utf-8'), + kc_paths, + alpha_values, + beta_i, + gamma, + stream_path.encode('utf-8'), + target_li_path.encode('utf-8'), + target_li_avail_path.encode('utf-8'), + target_l_sum_avail_path.encode('utf-8'), + target_aet_path.encode('utf-8'), + target_pi_path.encode('utf-8')] + + if algorithm == 'MFD': + run_calculate_local_recharge[MFD](*args) + else: # D8 + run_calculate_local_recharge[D8](*args) def route_baseflow_sum( - flow_dir_mfd_path, l_path, l_avail_path, l_sum_path, - stream_path, target_b_path, target_b_sum_path): + flow_dir_path, l_path, l_avail_path, l_sum_path, + stream_path, target_b_path, target_b_sum_path, algorithm): """Route Baseflow through MFD as described in Equation 11. Args: - flow_dir_mfd_path (string): path to a pygeoprocessing multiple flow + flow_dir_path (string): path to a pygeoprocessing multiple flow direction raster. l_path (string): path to local recharge raster. l_avail_path (string): path to local recharge raster that shows @@ -341,140 +138,24 @@ def route_baseflow_sum( Returns: None. """ - - # used for time-delayed logging - cdef time_t last_log_time - last_log_time = ctime(NULL) - cdef float target_nodata = -1e32 - cdef int stream_val, outlet - cdef float b_i, b_sum_i, l_j, l_avail_j, l_sum_j - cdef long xi, yi, xj, yj - cdef int flow_dir_i - cdef int flow_dir_nodata - cdef long raster_x_size, raster_y_size, xs_root, ys_root, xoff, yoff - cdef int n_dir - cdef int xs, ys, flow_dir_s, win_xsize, win_ysize - cdef int stream_nodata - cdef stack[pair[long, long]] work_stack - - # we know the PyGeoprocessing MFD raster flow dir type is a 32 bit int. - flow_dir_raster_info = pygeoprocessing.get_raster_info(flow_dir_mfd_path) - flow_dir_nodata = flow_dir_raster_info['nodata'][0] - raster_x_size, raster_y_size = flow_dir_raster_info['raster_size'] - stream_nodata = pygeoprocessing.get_raster_info(stream_path)['nodata'][0] pygeoprocessing.new_raster_from_base( - flow_dir_mfd_path, target_b_sum_path, gdal.GDT_Float32, + flow_dir_path, target_b_sum_path, gdal.GDT_Float32, [target_nodata], fill_value_list=[target_nodata]) pygeoprocessing.new_raster_from_base( - flow_dir_mfd_path, target_b_path, gdal.GDT_Float32, + flow_dir_path, target_b_path, gdal.GDT_Float32, [target_nodata], fill_value_list=[target_nodata]) - - cdef ManagedRaster target_b_sum_raster = ManagedRaster( - target_b_sum_path.encode('utf-8'), 1, 1) - cdef ManagedRaster target_b_raster = ManagedRaster( - target_b_path.encode('utf-8'), 1, 1) - cdef ManagedRaster l_raster = ManagedRaster(l_path.encode('utf-8'), 1, 0) - cdef ManagedRaster l_avail_raster = ManagedRaster(l_avail_path.encode('utf-8'), 1, 0) - cdef ManagedRaster l_sum_raster = ManagedRaster(l_sum_path.encode('utf-8'), 1, 0) - cdef ManagedFlowDirRaster flow_dir_mfd_raster = ManagedFlowDirRaster( - flow_dir_mfd_path.encode('utf-8'), 1, 0) - cdef ManagedRaster stream_raster = ManagedRaster(stream_path.encode('utf-8'), 1, 0) - - current_pixel = 0 - for offset_dict in pygeoprocessing.iterblocks( - (flow_dir_mfd_path, 1), offset_only=True, largest_block=0): - # use cython variables to avoid python overhead of dict values - win_xsize = offset_dict['win_xsize'] - win_ysize = offset_dict['win_ysize'] - xoff = offset_dict['xoff'] - yoff = offset_dict['yoff'] - for ys in xrange(win_ysize): - ys_root = yoff + ys - for xs in xrange(win_xsize): - xs_root = xoff + xs - flow_dir_s = flow_dir_mfd_raster.get(xs_root, ys_root) - if is_close(flow_dir_s, flow_dir_nodata): - current_pixel += 1 - continue - - # search for a pixel that has no downslope neighbors, - # or whose downslope neighbors all have nodata in the stream raster (?) - outlet = 1 - dn_neighbors = DownslopeNeighbors( - Pixel(flow_dir_mfd_raster, xs_root, ys_root)) - for neighbor in dn_neighbors: - stream_val = stream_raster.get(neighbor.x, neighbor.y) - if stream_val != stream_nodata: - outlet = 0 - break - if not outlet: - continue - work_stack.push(pair[long, long](xs_root, ys_root)) - - while work_stack.size() > 0: - xi = work_stack.top().first - yi = work_stack.top().second - work_stack.pop() - b_sum_i = target_b_sum_raster.get(xi, yi) - if not is_close(b_sum_i, target_nodata): - continue - - if ctime(NULL) - last_log_time > 5.0: - last_log_time = ctime(NULL) - LOGGER.info( - 'route base flow %.2f%% complete', - 100.0 * current_pixel / ( - raster_x_size * raster_y_size)) - - b_sum_i = 0.0 - downslope_defined = 1 - dn_neighbors = DownslopeNeighbors(Pixel(flow_dir_mfd_raster, xi, yi)) - for neighbor in dn_neighbors: - stream_val = stream_raster.get(neighbor.x, neighbor.y) - if stream_val: - b_sum_i += neighbor.flow_proportion - else: - b_sum_j = target_b_sum_raster.get(neighbor.x, neighbor.y) - if is_close(b_sum_j, target_nodata): - downslope_defined = 0 - break - l_j = l_raster.get(neighbor.x, neighbor.y) - l_avail_j = l_avail_raster.get(neighbor.x, neighbor.y) - l_sum_j = l_sum_raster.get(neighbor.x, neighbor.y) - - if l_sum_j != 0 and (l_sum_j - l_j) != 0: - b_sum_i += neighbor.flow_proportion * ( - (1 - l_avail_j / l_sum_j) * ( - b_sum_j / (l_sum_j - l_j))) - else: - b_sum_i += neighbor.flow_proportion - - if not downslope_defined: - continue - - l_i = l_raster.get(xi, yi) - l_sum_i = l_sum_raster.get(xi, yi) - b_sum_i = l_sum_i * b_sum_i - - if l_sum_i != 0: - b_i = max(b_sum_i * l_i / l_sum_i, 0.0) - else: - b_i = 0.0 - - target_b_raster.set(xi, yi, b_i) - target_b_sum_raster.set(xi, yi, b_sum_i) - - current_pixel += 1 - up_neighbors = UpslopeNeighbors(Pixel(flow_dir_mfd_raster, xi, yi)) - for neighbor in up_neighbors: - work_stack.push(pair[long, long](neighbor.x, neighbor.y)) - - target_b_sum_raster.close() - target_b_raster.close() - l_raster.close() - l_avail_raster.close() - l_sum_raster.close() - flow_dir_mfd_raster.close() - stream_raster.close() + args = [ + flow_dir_path.encode('utf-8'), + l_path.encode('utf-8'), + l_avail_path.encode('utf-8'), + l_sum_path.encode('utf-8'), + stream_path.encode('utf-8'), + target_b_path.encode('utf-8'), + target_b_sum_path.encode('utf-8')] + + if algorithm == 'MFD': + run_route_baseflow_sum[MFD](*args) + else: # D8 + run_route_baseflow_sum[D8](*args) diff --git a/src/natcap/invest/seasonal_water_yield/swy.h b/src/natcap/invest/seasonal_water_yield/swy.h new file mode 100644 index 000000000..a74821bd7 --- /dev/null +++ b/src/natcap/invest/seasonal_water_yield/swy.h @@ -0,0 +1,463 @@ +#include +#include +#include + +#include "ManagedRaster.h" + +// Calculate the rasters defined by equations [3]-[7]. +// +// Note all input rasters must be in the same coordinate system and +// have the same dimensions. +// +// Args: +// precip_paths: paths to monthly precipitation rasters. (model input) +// et0_paths: paths to monthly ET0 rasters. (model input) +// qf_m_paths: paths to monthly quickflow rasters calculated by +// Equation [1]. +// flow_dir_path: path to a flow direction raster (MFD or D8). Indicate MFD +// or D8 with the template argument. +// kc_paths: list of rasters of the monthly crop factor for the pixel. +// alpha_values: list of monthly alpha values (fraction of upslope annual +// available recharge that is available in each month) +// beta_i: fraction of the upgradient subsidy that is available +// for downgradient evapotranspiration. +// gamma: the fraction of pixel recharge that is available to +// downgradient pixels. +// stream_path: path to the stream raster where 1 is a stream, +// 0 is not, and nodata is outside of the DEM. +// target_li_path: created by this call, path to local recharge +// derived from the annual water budget. (Equation 3). +// target_li_avail_path: created by this call, path to raster +// indicating available recharge to a pixel. +// target_l_sum_avail_path: created by this call, the recursive +// upslope accumulation of target_li_avail_path. +// target_aet_path: created by this call, the annual actual +// evapotranspiration. +// target_pi_path: created by this call, the annual precipitation on +// a pixel. +template +void run_calculate_local_recharge( + vector precip_paths, + vector et0_paths, + vector qf_m_paths, + char* flow_dir_path, + vector kc_paths, + vector alpha_values, + float beta_i, + float gamma, + char* stream_path, + char* target_li_path, + char* target_li_avail_path, + char* target_l_sum_avail_path, + char* target_aet_path, + char* target_pi_path) { + long xs_root, ys_root, xoff, yoff; + long xi, yi, mfd_dir_sum; + long win_xsize, win_ysize; + double kc_m, pet_m, p_m, qf_m, et0_m, aet_i, p_i, qf_i, l_i; + double l_avail_i, l_avail_j, l_sum_avail_i, l_sum_avail_j; + bool upslope_defined; + + queue> work_queue; + + UpslopeNeighborsNoDivide up_neighbors; + DownslopeNeighbors dn_neighbors; + + // # used for time-delayed logging + // cdef time_t last_log_time + // last_log_time = ctime(NULL) + + ManagedFlowDirRaster flow_dir_raster = ManagedFlowDirRaster( + flow_dir_path, 1, 0); + NeighborTuple neighbor; + + // make sure that user input nodata values are defined + // set to -1 if not defined + // precipitation and evapotranspiration data should + // always be non-negative + vector et0_m_rasters; + vector et0_m_nodata_list; + for (auto et0_m_path: et0_paths) { + ManagedRaster et0_raster = ManagedRaster(et0_m_path, 1, 0); + et0_m_rasters.push_back(et0_raster); + if (et0_raster.hasNodata) { + et0_m_nodata_list.push_back(et0_raster.nodata); + } else { + et0_m_nodata_list.push_back(-1); + } + } + + vector precip_m_rasters; + vector precip_m_nodata_list; + for (auto precip_m_path: precip_paths) { + ManagedRaster precip_raster = ManagedRaster(precip_m_path, 1, 0); + precip_m_rasters.push_back(precip_raster); + if (precip_raster.hasNodata) { + precip_m_nodata_list.push_back(precip_raster.nodata); + } else { + precip_m_nodata_list.push_back(-1); + } + } + + vector qf_m_rasters; + for (auto qf_m_path: qf_m_paths) { + qf_m_rasters.push_back(ManagedRaster(qf_m_path, 1, 0)); + } + + vector kc_m_rasters; + for (auto kc_m_path: kc_paths) { + kc_m_rasters.push_back(ManagedRaster(kc_m_path, 1, 0)); + } + + ManagedRaster target_li_raster = ManagedRaster(target_li_path, 1, 1); + ManagedRaster target_li_avail_raster = ManagedRaster(target_li_avail_path, 1, 1); + ManagedRaster target_l_sum_avail_raster = ManagedRaster(target_l_sum_avail_path, 1, 1); + ManagedRaster target_aet_raster = ManagedRaster(target_aet_path, 1, 1); + ManagedRaster target_pi_raster = ManagedRaster(target_pi_path, 1, 1); + + double target_nodata = -1e32; + + // efficient way to calculate ceiling division: + // a divided by b rounded up = (a + (b - 1)) / b + // note that / represents integer floor division + // https://stackoverflow.com/a/62032709/14451410 + int n_col_blocks = (flow_dir_raster.raster_x_size + (flow_dir_raster.block_xsize - 1)) / flow_dir_raster.block_xsize; + int n_row_blocks = (flow_dir_raster.raster_y_size + (flow_dir_raster.block_ysize - 1)) / flow_dir_raster.block_ysize; + + for (int row_block_index = 0; row_block_index < n_row_blocks; row_block_index++) { + yoff = row_block_index * flow_dir_raster.block_ysize; + win_ysize = flow_dir_raster.raster_y_size - yoff; + if (win_ysize > flow_dir_raster.block_ysize) { + win_ysize = flow_dir_raster.block_ysize; + } + for (int col_block_index = 0; col_block_index < n_col_blocks; col_block_index++) { + xoff = col_block_index * flow_dir_raster.block_xsize; + win_xsize = flow_dir_raster.raster_x_size - xoff; + if (win_xsize > flow_dir_raster.block_xsize) { + win_xsize = flow_dir_raster.block_xsize; + } + + // if ctime(NULL) - last_log_time > 5.0: + // last_log_time = ctime(NULL) + // LOGGER.info('Sediment deposition %.2f%% complete', 100 * ( + // n_pixels_processed / float(flow_dir_raster.raster_x_size * flow_dir_raster.raster_y_size))) + + for (int row_index = 0; row_index < win_ysize; row_index++) { + ys_root = yoff + row_index; + for (int col_index = 0; col_index < win_xsize; col_index++) { + xs_root = xoff + col_index; + + if (flow_dir_raster.get(xs_root, ys_root) == flow_dir_raster.nodata) { + continue; + } + + if (flow_dir_raster.is_local_high_point(xs_root, ys_root)) { + work_queue.push(pair(xs_root, ys_root)); + } + + while (work_queue.size() > 0) { + xi = work_queue.front().first; + yi = work_queue.front().second; + work_queue.pop(); + + l_sum_avail_i = target_l_sum_avail_raster.get(xi, yi); + if (not is_close(l_sum_avail_i, target_nodata)) { + // already defined + continue; + } + + // Equation 7, calculate L_sum_avail_i if possible, skip + // otherwise + upslope_defined = true; + // initialize to 0 so we indicate we haven't tracked any + // mfd values yet + l_sum_avail_i = 0.0; + mfd_dir_sum = 0; + up_neighbors = UpslopeNeighborsNoDivide(Pixel(flow_dir_raster, xi, yi)); + for (auto neighbor: up_neighbors) { + // pixel flows inward, check upslope + l_sum_avail_j = target_l_sum_avail_raster.get( + neighbor.x, neighbor.y); + if (is_close(l_sum_avail_j, target_nodata)) { + upslope_defined = false; + break; + } + l_avail_j = target_li_avail_raster.get( + neighbor.x, neighbor.y); + // A step of Equation 7 + l_sum_avail_i += ( + l_sum_avail_j + l_avail_j) * neighbor.flow_proportion; + mfd_dir_sum += static_cast(neighbor.flow_proportion); + } + // calculate l_sum_avail_i by summing all the valid + // directions then normalizing by the sum of the mfd + // direction weights (Equation 8) + if (upslope_defined) { + // Equation 7 + if (mfd_dir_sum > 0) { + l_sum_avail_i /= static_cast(mfd_dir_sum); + } + target_l_sum_avail_raster.set(xi, yi, l_sum_avail_i); + } else { + // if not defined, we'll get it on another pass + continue; + } + + aet_i = 0; + p_i = 0; + qf_i = 0; + + for (int m_index = 0; m_index < 12; m_index++) { + p_m = precip_m_rasters[m_index].get(xi, yi); + if (not is_close(p_m, precip_m_rasters[m_index].nodata)) { + p_i += p_m; + } else { + p_m = 0; + } + + qf_m = qf_m_rasters[m_index].get(xi, yi); + if (not is_close(qf_m, qf_m_rasters[m_index].nodata)) { + qf_i += qf_m; + } else { + qf_m = 0; + } + + kc_m = kc_m_rasters[m_index].get(xi, yi); + pet_m = 0; + et0_m = et0_m_rasters[m_index].get(xi, yi); + if (not ( + is_close(kc_m, kc_m_rasters[m_index].nodata) or + is_close(et0_m, et0_m_rasters[m_index].nodata))) { + // Equation 6 + pet_m = kc_m * et0_m; + } + + // Equation 4/5 + aet_i += min( + pet_m, + p_m - qf_m + + alpha_values[m_index]*beta_i*l_sum_avail_i); + } + l_i = (p_i - qf_i - aet_i); + l_avail_i = min(gamma * l_i, l_i); + + target_pi_raster.set(xi, yi, p_i); + target_aet_raster.set(xi, yi, aet_i); + target_li_raster.set(xi, yi, l_i); + target_li_avail_raster.set(xi, yi, l_avail_i); + + dn_neighbors = DownslopeNeighbors(Pixel(flow_dir_raster, xi, yi)); + for (auto neighbor: dn_neighbors) { + work_queue.push(pair(neighbor.x, neighbor.y)); + } + } + } + } + } + } + + flow_dir_raster.close(); + target_li_raster.close(); + target_li_avail_raster.close(); + target_l_sum_avail_raster.close(); + target_aet_raster.close(); + target_pi_raster.close(); + for (int i = 0; i < 12; i++) { + et0_m_rasters[i].close(); + precip_m_rasters[i].close(); + qf_m_rasters[i].close(); + kc_m_rasters[i].close(); + } +} + +// Route Baseflow as described in Equation 11. +// Args: +// flow_dir_path: path to a MFD or D8 flow direction raster. +// l_path: path to local recharge raster. +// l_avail_path: path to local recharge raster that shows +// recharge available to the pixel. +// l_sum_path: path to upslope sum of l_path. +// stream_path: path to stream raster, 1 stream, 0 no stream, +// and nodata. +// target_b_path: path to created raster for per-pixel baseflow. +// target_b_sum_path: path to created raster for per-pixel +// upslope sum of baseflow. +template +void run_route_baseflow_sum( + char* flow_dir_path, + char* l_path, + char* l_avail_path, + char* l_sum_path, + char* stream_path, + char* target_b_path, + char* target_b_sum_path) { + + // used for time-delayed logging + // cdef time_t last_log_time + // last_log_time = ctime(NULL) + + float target_nodata = static_cast(-1e32); + double b_i, b_sum_i, b_sum_j, l_j, l_avail_j, l_sum_j; + double l_i, l_sum_i; + long xi, yi, flow_dir_sum; + long xs_root, ys_root, xoff, yoff; + int win_xsize, win_ysize; + stack> work_stack; + bool outlet, downslope_defined; + + ManagedRaster target_b_sum_raster = ManagedRaster(target_b_sum_path, 1, 1); + ManagedRaster target_b_raster = ManagedRaster(target_b_path, 1, 1); + ManagedRaster l_raster = ManagedRaster(l_path, 1, 0); + ManagedRaster l_avail_raster = ManagedRaster(l_avail_path, 1, 0); + ManagedRaster l_sum_raster = ManagedRaster(l_sum_path, 1, 0); + ManagedFlowDirRaster flow_dir_raster = ManagedFlowDirRaster(flow_dir_path, 1, 0); + ManagedRaster stream_raster = ManagedRaster(stream_path, 1, 0); + + UpslopeNeighbors up_neighbors; + DownslopeNeighbors dn_neighbors; + DownslopeNeighborsNoSkip dn_neighbors_no_skip; + NeighborTuple neighbor; + + // int current_pixel = 0; + + // efficient way to calculate ceiling division: + // a divided by b rounded up = (a + (b - 1)) / b + // note that / represents integer floor division + // https://stackoverflow.com/a/62032709/14451410 + int n_col_blocks = (flow_dir_raster.raster_x_size + (flow_dir_raster.block_xsize - 1)) / flow_dir_raster.block_xsize; + int n_row_blocks = (flow_dir_raster.raster_y_size + (flow_dir_raster.block_ysize - 1)) / flow_dir_raster.block_ysize; + + for (int row_block_index = 0; row_block_index < n_row_blocks; row_block_index++) { + yoff = row_block_index * flow_dir_raster.block_ysize; + win_ysize = flow_dir_raster.raster_y_size - yoff; + if (win_ysize > flow_dir_raster.block_ysize) { + win_ysize = flow_dir_raster.block_ysize; + } + for (int col_block_index = 0; col_block_index < n_col_blocks; col_block_index++) { + xoff = col_block_index * flow_dir_raster.block_xsize; + win_xsize = flow_dir_raster.raster_x_size - xoff; + if (win_xsize > flow_dir_raster.block_xsize) { + win_xsize = flow_dir_raster.block_xsize; + } + + // if ctime(NULL) - last_log_time > 5.0: + // last_log_time = ctime(NULL) + // LOGGER.info('Sediment deposition %.2f%% complete', 100 * ( + // n_pixels_processed / float(flow_dir_raster.raster_x_size * flow_dir_raster.raster_y_size))) + + for (int row_index = 0; row_index < win_ysize; row_index++) { + ys_root = yoff + row_index; + for (int col_index = 0; col_index < win_xsize; col_index++) { + xs_root = xoff + col_index; + + if (static_cast(flow_dir_raster.get(xs_root, ys_root)) == + static_cast(flow_dir_raster.nodata)) { + // current_pixel += 1; + continue; + } + + // search for a pixel that has no downslope neighbors, + // or whose downslope neighbors all have nodata in the stream raster (?) + outlet = true; + dn_neighbors = DownslopeNeighbors(Pixel(flow_dir_raster, xs_root, ys_root)); + for (auto neighbor: dn_neighbors) { + if (static_cast(stream_raster.get(neighbor.x, neighbor.y)) != + static_cast(stream_raster.nodata)) { + outlet = 0; + break; + } + } + if (not outlet) { + continue; + } + work_stack.push(pair(xs_root, ys_root)); + + while (work_stack.size() > 0) { + xi = work_stack.top().first; + yi = work_stack.top().second; + work_stack.pop(); + b_sum_i = target_b_sum_raster.get(xi, yi); + if (not is_close(b_sum_i, target_nodata)) { + continue; + } + + // if ctime(NULL) - last_log_time > 5.0: + // last_log_time = ctime(NULL) + // LOGGER.info( + // 'route base flow %.2f%% complete', + // 100.0 * current_pixel / ( + // flow_dir_raster.raster_x_size * flow_dir_raster.raster_y_size)) + + b_sum_i = 0; + downslope_defined = true; + dn_neighbors_no_skip = DownslopeNeighborsNoSkip(Pixel(flow_dir_raster, xi, yi)); + flow_dir_sum = 0; + for (auto neighbor: dn_neighbors_no_skip) { + flow_dir_sum += static_cast(neighbor.flow_proportion); + + if (neighbor.x < 0 or neighbor.x >= flow_dir_raster.raster_x_size or + neighbor.y < 0 or neighbor.y >= flow_dir_raster.raster_y_size) { + continue; + } + + if (static_cast(stream_raster.get(neighbor.x, neighbor.y))) { + b_sum_i += neighbor.flow_proportion; + } else { + b_sum_j = target_b_sum_raster.get(neighbor.x, neighbor.y); + if (is_close(b_sum_j, target_nodata)) { + downslope_defined = false; + break; + } + l_j = l_raster.get(neighbor.x, neighbor.y); + l_avail_j = l_avail_raster.get(neighbor.x, neighbor.y); + l_sum_j = l_sum_raster.get(neighbor.x, neighbor.y); + + if (l_sum_j != 0 and (l_sum_j - l_j) != 0) { + b_sum_i += neighbor.flow_proportion * ( + (1 - l_avail_j / l_sum_j) * ( + b_sum_j / (l_sum_j - l_j))); + } else { + b_sum_i += neighbor.flow_proportion; + } + } + } + + if (not downslope_defined) { + continue; + } + l_i = l_raster.get(xi, yi); + l_sum_i = l_sum_raster.get(xi, yi); + + if (flow_dir_sum > 0) { + b_sum_i = l_sum_i * b_sum_i / flow_dir_sum; + } + + + if (l_sum_i != 0) { + b_i = max(b_sum_i * l_i / l_sum_i, 0.0); + } else { + b_i = 0; + } + + target_b_raster.set(xi, yi, b_i); + target_b_sum_raster.set(xi, yi, b_sum_i); + + // current_pixel += 1; + up_neighbors = UpslopeNeighbors(Pixel(flow_dir_raster, xi, yi)); + for (auto neighbor: up_neighbors) { + work_stack.push(pair(neighbor.x, neighbor.y)); + } + } + } + } + } + } + + target_b_sum_raster.close(); + target_b_raster.close(); + l_raster.close(); + l_avail_raster.close(); + l_sum_raster.close(); + flow_dir_raster.close(); + stream_raster.close(); +} diff --git a/src/natcap/invest/seasonal_water_yield/swy.pxd b/src/natcap/invest/seasonal_water_yield/swy.pxd new file mode 100644 index 000000000..712ee45ab --- /dev/null +++ b/src/natcap/invest/seasonal_water_yield/swy.pxd @@ -0,0 +1,28 @@ +from libcpp.vector cimport vector + +cdef extern from "swy.h": + void run_calculate_local_recharge[T]( + vector[char*], # precip_path_list + vector[char*], # et0_path_list + vector[char*], # qf_m_path_list + char*, # flow_dir_mfd_path + vector[char*], # kc_path_list + vector[float], # alpha_values + float, # beta_i + float, # gamma + char*, # stream_path + char*, # target_li_path + char*, # target_li_avail_path + char*, # target_l_sum_avail_path + char*, # target_aet_path + char* # target_pi_path + ) + + void run_route_baseflow_sum[T]( + char*, + char*, + char*, + char*, + char*, + char*, + char*) diff --git a/tests/test_ndr.py b/tests/test_ndr.py index a719c7dd5..47f15bfb2 100644 --- a/tests/test_ndr.py +++ b/tests/test_ndr.py @@ -48,6 +48,7 @@ def generate_base_args(workspace_dir): 'watersheds_path': os.path.join(REGRESSION_DATA, 'input', 'watersheds.shp'), 'workspace_dir': workspace_dir, + 'algorithm': 'MFD' } return args.copy() @@ -247,6 +248,61 @@ def test_base_regression(self): args['workspace_dir'], 'intermediate_outputs', 'what_drains_to_stream.tif'))) + def test_base_regression_d8(self): + """NDR base regression test on sample data in D8 mode. + + Execute NDR with sample data and checks that the output files are + generated and that the aggregate shapefile fields are the same as the + regression case. + """ + from natcap.invest.ndr import ndr + + # use predefined directory so test can clean up files during teardown + args = NDRTests.generate_base_args(self.workspace_dir) + args['algorithm'] = 'D8' + # make an empty output shapefile on top of where the new output + # shapefile should reside to ensure the model overwrites it + with open( + os.path.join(self.workspace_dir, 'watershed_results_ndr.gpkg'), + 'wb') as f: + f.write(b'') + ndr.execute(args) + + result_vector = ogr.Open(os.path.join( + args['workspace_dir'], 'watershed_results_ndr.gpkg')) + result_layer = result_vector.GetLayer() + result_feature = result_layer.GetFeature(1) + result_layer = None + result_vector = None + mismatch_list = [] + # these values were generated by manual inspection of regression + # results + for field, expected_value in [ + ('p_surface_load', 41.921860), + ('p_surface_export', 4.985530), + ('n_surface_load', 2978.519775), + ('n_surface_export', 329.733093), + ('n_subsurface_load', 28.614094), + ('n_subsurface_export', 12.637195), + ('n_total_export', 339.975503)]: + val = result_feature.GetField(field) + if not numpy.isclose(val, expected_value): + mismatch_list.append( + (field, 'expected: %f' % expected_value, + 'actual: %f' % val)) + result_feature = None + if mismatch_list: + raise RuntimeError("results not expected: %s" % mismatch_list) + + # We only need to test that the drainage mask exists. Functionality + # for that raster is tested in SDR. + self.assertTrue( + os.path.exists( + os.path.join( + args['workspace_dir'], 'intermediate_outputs', + 'what_drains_to_stream.tif'))) + + def test_regression_undefined_nodata(self): """NDR test when DEM, LULC and runoff proxy have undefined nodata.""" from natcap.invest.ndr import ndr @@ -362,6 +418,7 @@ def test_validation(self): 'k_param', 'watersheds_path', 'subsurface_eff_n', + 'algorithm' ] self.assertEqual(set(invalid_args), set(expected_missing_args)) diff --git a/tests/test_sdr.py b/tests/test_sdr.py index 341076a08..2593913f4 100644 --- a/tests/test_sdr.py +++ b/tests/test_sdr.py @@ -71,6 +71,7 @@ def generate_base_args(workspace_dir): 'watersheds_path': os.path.join(SAMPLE_DATA, 'watersheds.shp'), 'workspace_dir': workspace_dir, 'n_workers': -1, + 'algorithm': 'MFD' } return args @@ -105,15 +106,6 @@ def test_sdr_validation_wrong_types(self): validate_result, ['GDAL raster', 'GDAL vector']): self.assertTrue(phrase in error_msg) - def test_sdr_validation_missing_key(self): - """SDR test validation that's missing keys.""" - from natcap.invest.sdr import sdr - - # use predefined directory so test can clean up files during teardown - args = {} - validation_warnings = sdr.validate(args, limit_to=None) - self.assertEqual(len(validation_warnings[0][0]), 12) - def test_sdr_validation_key_no_value(self): """SDR test validation that's missing a value on a key.""" from natcap.invest.sdr import sdr @@ -175,6 +167,57 @@ def test_base_regression(self): self.assertEqual( numpy.count_nonzero(sed_dep_array[negative_non_nodata_mask]), 0) + def test_base_regression_d8(self): + """SDR base regression test on sample data in D8 mode. + + Execute SDR with sample data and checks that the output files are + generated and that the aggregate shapefile fields are the same as the + regression case. + """ + from natcap.invest.sdr import sdr + + # use predefined directory so test can clean up files during teardown + args = SDRTests.generate_base_args(self.workspace_dir) + args['algorithm'] = 'D8' + args['threshold_flow_accumulation'] = 100 + # make args explicit that this is a base run of SWY + + sdr.execute(args) + expected_results = { + 'usle_tot': 2.520746, + 'sed_export': 0.187428, + 'sed_dep': 2.300645, + 'avoid_exp': 19283.767578, + 'avoid_eros': 263415, + } + + vector_path = os.path.join( + args['workspace_dir'], 'watershed_results_sdr.shp') + assert_expected_results_in_vector(expected_results, vector_path) + + # We only need to test that the drainage mask exists. Functionality + # for that raster is tested elsewhere + self.assertTrue( + os.path.exists( + os.path.join( + args['workspace_dir'], 'intermediate_outputs', + 'what_drains_to_stream.tif'))) + + # Check that sed_deposition does not have any negative, non-nodata + # values, even if they are very small. + sed_deposition_path = os.path.join(args['workspace_dir'], + 'sed_deposition.tif') + sed_dep_nodata = pygeoprocessing.get_raster_info( + sed_deposition_path)['nodata'][0] + sed_dep_array = pygeoprocessing.raster_to_numpy_array( + sed_deposition_path) + negative_non_nodata_mask = ( + (~numpy.isclose(sed_dep_array, sed_dep_nodata)) & + (sed_dep_array < 0)) + self.assertEqual( + numpy.count_nonzero(sed_dep_array[negative_non_nodata_mask]), 0) + + def test_regression_with_undefined_nodata(self): """SDR base regression test with undefined nodata values. diff --git a/tests/test_seasonal_water_yield_regression.py b/tests/test_seasonal_water_yield_regression.py index 2b4359cba..bdc5774d6 100644 --- a/tests/test_seasonal_water_yield_regression.py +++ b/tests/test_seasonal_water_yield_regression.py @@ -423,6 +423,7 @@ def test_ambiguous_precip_data(self): 'user_defined_climate_zones': False, 'user_defined_local_recharge': False, 'monthly_alpha': False, + 'algorithm': 'MFD' } watershed_shp_path = os.path.join(args['workspace_dir'], @@ -484,6 +485,7 @@ def test_precip_data_missing(self): 'user_defined_climate_zones': False, 'user_defined_local_recharge': False, 'monthly_alpha': False, + 'algorithm': 'MFD' } watershed_shp_path = os.path.join(args['workspace_dir'], @@ -584,6 +586,7 @@ def test_duplicate_aoi_assertion(self): 'user_defined_climate_zones': False, 'user_defined_local_recharge': False, 'monthly_alpha': False, + 'algorithm': 'MFD' } biophysical_csv_path = os.path.join(args['workspace_dir'], @@ -643,6 +646,7 @@ def generate_base_args(workspace_dir): 'results_suffix': '', 'threshold_flow_accumulation': '50', 'workspace_dir': workspace_dir, + 'algorithm': 'MFD' } watershed_shp_path = os.path.join(workspace_dir, 'watershed.shp') @@ -729,6 +733,58 @@ def test_base_regression(self): os.path.join(args['workspace_dir'], 'aggregated_results_swy.shp'), agg_results_csv_path) + def test_base_regression_d8(self): + """SWY base regression test on sample data in D8 mode. + + Executes SWY in default mode and checks that the output files are + generated and that the aggregate shapefile fields are the same as the + regression case. + """ + from natcap.invest.seasonal_water_yield import seasonal_water_yield + + # use predefined directory so test can clean up files during teardown + args = SeasonalWaterYieldRegressionTests.generate_base_args( + self.workspace_dir) + + # Ensure the model can pass when a nodata value is not defined. + size = 100 + lulc_array = numpy.zeros((size, size), dtype=numpy.int8) + lulc_array[size // 2:, :] = 1 + + driver = gdal.GetDriverByName('GTiff') + new_raster = driver.Create( + args['lulc_raster_path'], lulc_array.shape[0], + lulc_array.shape[1], 1, gdal.GDT_Byte) + band = new_raster.GetRasterBand(1) + band.WriteArray(lulc_array) + geotransform = [1180000, 1, 0, 690000, 0, -1] + new_raster.SetGeoTransform(geotransform) + band = None + new_raster = None + driver = None + + # make args explicit that this is a base run of SWY + args['user_defined_climate_zones'] = False + args['user_defined_local_recharge'] = False + args['monthly_alpha'] = False + args['results_suffix'] = '' + args['algorithm'] = 'D8' + + seasonal_water_yield.execute(args) + + result_vector = ogr.Open(os.path.join( + args['workspace_dir'], 'aggregated_results_swy.shp')) + result_layer = result_vector.GetLayer() + result_feature = result_layer.GetFeature(0) + mismatch_list = [] + for field, expected_value in [('vri_sum', 1), ('qb', 52.9128)]: + val = result_feature.GetField(field) + if not numpy.isclose(val, expected_value): + mismatch_list.append( + (field, f'expected: {expected_value}', f'actual: {val}')) + if mismatch_list: + raise RuntimeError(f'results not expected: {mismatch_list}') + def test_base_regression_nodata_inf(self): """SWY base regression test on sample data with really small nodata. @@ -1263,13 +1319,14 @@ def test_local_recharge_undefined_nodata(self): seasonal_water_yield_core.calculate_local_recharge( [precip_path for i in range(12)], [et0_path for i in range(12)], [quickflow_path for i in range(12)], flow_dir_path, - [kc_path for i in range(12)], {i: 0.5 for i in range(12)}, 0.5, + [kc_path for i in range(12)], {i: 0.5 for i in range(1, 13)}, 0.5, 0.5, stream_path, os.path.join(self.workspace_dir, 'target_li_path.tif'), os.path.join(self.workspace_dir, 'target_li_avail_path.tif'), os.path.join(self.workspace_dir, 'target_l_sum_avail_path.tif'), os.path.join(self.workspace_dir, 'target_aet_path.tif'), - os.path.join(self.workspace_dir, 'target_precip_path.tif')) + os.path.join(self.workspace_dir, 'target_precip_path.tif'), + algorithm='MFD') class SWYValidationTests(unittest.TestCase): @@ -1295,6 +1352,7 @@ def setUp(self): 'precip_dir', 'threshold_flow_accumulation', 'user_defined_local_recharge', + 'algorithm' ] def tearDown(self): diff --git a/workbench/src/renderer/ui_config.js b/workbench/src/renderer/ui_config.js index b9a7af3d6..452a8e180 100644 --- a/workbench/src/renderer/ui_config.js +++ b/workbench/src/renderer/ui_config.js @@ -202,7 +202,7 @@ const UI_SPEC = { ['dem_path', 'lulc_path', 'runoff_proxy_path', 'watersheds_path', 'biophysical_table_path'], ['calc_p'], ['calc_n', 'subsurface_critical_length_n', 'subsurface_eff_n'], - ['threshold_flow_accumulation', 'k_param'], + ['algorithm', 'threshold_flow_accumulation', 'k_param'], ], enabledFunctions: { subsurface_critical_length_n: isSufficient.bind(null, 'calc_n'), @@ -286,7 +286,7 @@ const UI_SPEC = { ['dem_path', 'erosivity_path', 'erodibility_path'], ['lulc_path', 'biophysical_table_path'], ['watersheds_path', 'drainage_path'], - ['threshold_flow_accumulation', 'k_param', 'sdr_max', 'ic_0_param', 'l_max'], + ['algorithm', 'threshold_flow_accumulation', 'k_param', 'sdr_max', 'ic_0_param', 'l_max'], ], }, seasonal_water_yield: { @@ -294,7 +294,7 @@ const UI_SPEC = { ['workspace_dir', 'results_suffix'], ['lulc_raster_path', 'biophysical_table_path'], ['dem_raster_path', 'aoi_path'], - ['threshold_flow_accumulation', 'beta_i', 'gamma'], + ['algorithm', 'threshold_flow_accumulation', 'beta_i', 'gamma'], ['user_defined_local_recharge', 'l_path', 'et0_dir', 'precip_dir', 'soil_group_path'], ['monthly_alpha', 'alpha_m', 'monthly_alpha_path'], ['user_defined_climate_zones', 'rain_events_table_path', 'climate_zone_table_path', 'climate_zone_raster_path'],