diff --git a/.github/workflows/build-and-test.yml b/.github/workflows/build-and-test.yml index 155b27add4..172fede3a8 100644 --- a/.github/workflows/build-and-test.yml +++ b/.github/workflows/build-and-test.yml @@ -119,7 +119,7 @@ jobs: run: | # uninstall InVEST if it was already in the restored cache python -m pip uninstall -y natcap.invest - python -m build --wheel + NATCAP_INVEST_GDAL_LIB_PATH="$CONDA_PREFIX/Library" python -m build --wheel ls -la dist pip install $(find dist -name "natcap.invest*.whl") @@ -180,12 +180,12 @@ jobs: # .cpp files in addition to the .pyx files. # # Elevating any python warnings to errors to catch build issues ASAP. - python -W error -m build --sdist + NATCAP_INVEST_GDAL_LIB_PATH="$CONDA_PREFIX/Library" python -W error -m build --sdist - name: Install from source distribution run : | # Install natcap.invest from the sdist in dist/ - pip install $(find dist -name "natcap[._-]invest*") + NATCAP_INVEST_GDAL_LIB_PATH="$CONDA_PREFIX/Library" pip install $(find dist -name "natcap[._-]invest*") # Model tests should cover model functionality, we just want # to be sure that we can import `natcap.invest` here. @@ -235,7 +235,7 @@ jobs: requirements: ${{ env.CONDA_DEFAULT_DEPENDENCIES }} pytest - name: Make install - run: make install + run: NATCAP_INVEST_GDAL_LIB_PATH="$CONDA_PREFIX/Library" make install - name: Validate sample data run: make validate_sampledata @@ -267,7 +267,7 @@ jobs: requirements: ${{ env.CONDA_DEFAULT_DEPENDENCIES }} - name: Make install - run: make install + run: NATCAP_INVEST_GDAL_LIB_PATH="$CONDA_PREFIX/Library" make install - name: Set up node uses: actions/setup-node@v3 @@ -328,7 +328,7 @@ jobs: requirements: ${{ env.CONDA_DEFAULT_DEPENDENCIES }} pandoc - name: Make install - run: make install + run: NATCAP_INVEST_GDAL_LIB_PATH="$CONDA_PREFIX/Library" make install # Not caching chocolatey packages because the cache may not be reliable # https://github.com/chocolatey/choco/issues/2134 diff --git a/README.rst b/README.rst index 0e98605a3d..12527c21ab 100644 --- a/README.rst +++ b/README.rst @@ -105,6 +105,11 @@ Building InVEST Distributions ----------------------------- Once the required tools and packages are available, we can build InVEST. +On Windows, you must indicate the location of the GDAL libraries with the environment +variable ``NATCAP_INVEST_GDAL_LIB_PATH``. If you are using conda to manage dependencies +as we recommend, you can add ``NATCAP_INVEST_GDAL_LIB_PATH="$CONDA_PREFIX/Library"`` to +the commands below. (On Mac and Linux, the gdal library path is determined for you +automatically using ``gdal-config``, which isn't available on Windows.) Building ``natcap.invest`` python package diff --git a/requirements-dev.txt b/requirements-dev.txt index 75b210cc1b..07ae8b13ab 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -14,7 +14,7 @@ virtualenv>=12.0.1 pytest -pytest-subtests +pytest-subtests<0.14.0 wheel>=0.27.0 pypiwin32; sys_platform == 'win32' # pip-only diff --git a/setup.py b/setup.py index 3c849b51df..531b8fa762 100644 --- a/setup.py +++ b/setup.py @@ -15,14 +15,25 @@ if (not req.startswith(('#', 'hg+', 'git+')) and len(req.strip()) > 0)] -# Since OSX Mavericks, the stdlib has been renamed. So if we're on OSX, we -# need to be sure to define which standard c++ library to use. I don't have -# access to a pre-Mavericks mac, so hopefully this won't break on someone's -# older system. Tested and it works on Mac OSX Catalina. compiler_and_linker_args = [] -if platform.system() == 'Darwin': - compiler_and_linker_args = ['-stdlib=libc++'] - +include_dirs = [numpy.get_include(), 'src/natcap/invest/managed_raster'] +if platform.system() == 'Windows': + compiler_args = ['/std:c++20'] + if 'NATCAP_INVEST_GDAL_LIB_PATH' not in os.environ: + raise RuntimeError( + 'env variable NATCAP_INVEST_GDAL_LIB_PATH is not defined. ' + 'This env variable is required when building on Windows. If ' + 'using conda to manage your gdal installation, you may set ' + 'NATCAP_INVEST_GDAL_LIB_PATH="$CONDA_PREFIX/Library".') + library_dirs = [f'{os.environ["NATCAP_INVEST_GDAL_LIB_PATH"]}/lib'] + include_dirs.append(f'{os.environ["NATCAP_INVEST_GDAL_LIB_PATH"]}/include') +else: + library_dirs = [] + compiler_args = [] + compiler_and_linker_args = ['-std=c++20'] + library_dirs = [subprocess.run( + ['gdal-config', '--libs'], capture_output=True, text=True + ).stdout.split()[0][2:]] # get the first argument which is the library path class build_py(_build_py): """Command to compile translation message catalogs before building.""" @@ -50,13 +61,14 @@ def run(self): Extension( name=f'natcap.invest.{package}.{module}', sources=[f'src/natcap/invest/{package}/{module}.pyx'], - include_dirs=[numpy.get_include()] + ['src/natcap/invest/managed_raster'], - extra_compile_args=compiler_args + compiler_and_linker_args, + include_dirs=include_dirs, + extra_compile_args=compiler_args + package_compiler_args + compiler_and_linker_args, extra_link_args=compiler_and_linker_args, language='c++', + libraries=['gdal'], + library_dirs=library_dirs, define_macros=[("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION")] - ) for package, module, compiler_args in [ - ('managed_raster', 'managed_raster', []), + ) for package, module, package_compiler_args in [ ('delineateit', 'delineateit_core', []), ('recreation', 'out_of_core_quadtree', []), # clang-14 defaults to -ffp-contract=on, which causes the diff --git a/src/natcap/invest/managed_raster/LRUCache.h b/src/natcap/invest/managed_raster/LRUCache.h index d8c93c4c88..26b22f2e31 100644 --- a/src/natcap/invest/managed_raster/LRUCache.h +++ b/src/natcap/invest/managed_raster/LRUCache.h @@ -8,66 +8,72 @@ using namespace std; template >::iterator, - typename MapIter = typename map::iterator > class LRUCache{ -private: - // item_list keeps track of the order of which elements have been accessed - // element at begin is most recent, element at end is least recent. - // first element in the pair is its key while the second is the element - list< pair > item_list; - // item_map maps an element's key to its location in the `item_list` - // used to make lookups O(log n) time - map item_map; - size_t cache_size; -private: - void clean(list< pair > &removed_value_list){ - while(item_map.size()>cache_size){ - ListIter last_it = item_list.end(); last_it --; - removed_value_list.push_back( - make_pair(last_it->first, last_it->second)); - item_map.erase(last_it->first); - item_list.pop_back(); - } - }; -public: - LRUCache(int cache_size_):cache_size(cache_size_){ - ; - }; + typename ListIter = typename list>::iterator, + typename MapIter = typename map::iterator> class LRUCache { + private: + // item_list keeps track of the order of which elements have been accessed + // element at begin is most recent, element at end is least recent. + // first element in the pair is its key while the second is the element + list> item_list; + // item_map maps an element's key to its location in the `item_list` + // used to make lookups O(log n) time + map item_map; + size_t cache_size; - ListIter begin() { - return item_list.begin(); + void clean(list> &removed_value_list) { + while(item_map.size() > cache_size) { + ListIter last_it = item_list.end(); + last_it--; + removed_value_list.push_back( + make_pair(last_it->first, last_it->second)); + item_map.erase(last_it->first); + item_list.pop_back(); } + }; + public: + LRUCache(int cache_size_):cache_size(cache_size_) { + ; + }; - ListIter end() { - return item_list.end(); + ListIter begin() { + return item_list.begin(); + } + + ListIter end() { + return item_list.end(); + } + + // Insert a new key-value pair into the cache. + void put( + const KEY_T &key, const VAL_T &val, + list> &removed_value_list) { + MapIter it = item_map.find(key); + if(it != item_map.end()){ + // it's already in the cache, delete the location in the item + // list and in the lookup map + item_list.erase(it->second); + item_map.erase(it); } + // insert a new item in the front since it's most recently used + item_list.push_front(make_pair(key, val)); + // record its iterator in the map + item_map.insert(make_pair(key, item_list.begin())); + // possibly remove any elements that have exceeded the cache size + return clean(removed_value_list); + }; + + // Return whether a key exists in the cache. + bool exist(const KEY_T &key) { + return (item_map.count(key) > 0); + }; - void put( - const KEY_T &key, const VAL_T &val, - list< pair > &removed_value_list) { - MapIter it = item_map.find(key); - if(it != item_map.end()){ - // it's already in the cache, delete the location in the item - // list and in the lookup map - item_list.erase(it->second); - item_map.erase(it); - } - // insert a new item in the front since it's most recently used - item_list.push_front(make_pair(key,val)); - // record its iterator in the map - item_map.insert(make_pair(key, item_list.begin())); - // possibly remove any elements that have exceeded the cache size - return clean(removed_value_list); - }; - bool exist(const KEY_T &key){ - return (item_map.count(key)>0); - }; - VAL_T& get(const KEY_T &key){ - MapIter it = item_map.find(key); - assert(it!=item_map.end()); - // move the element to the front of the list - item_list.splice(item_list.begin(), item_list, it->second); - return it->second->second; - }; + // Return the cached value associated with a key. + VAL_T& get(const KEY_T &key) { + MapIter it = item_map.find(key); + assert(it != item_map.end()); + // move the element to the front of the list + item_list.splice(item_list.begin(), item_list, it->second); + return it->second->second; + }; }; #endif diff --git a/src/natcap/invest/managed_raster/ManagedRaster.h b/src/natcap/invest/managed_raster/ManagedRaster.h new file mode 100644 index 0000000000..86955e99cb --- /dev/null +++ b/src/natcap/invest/managed_raster/ManagedRaster.h @@ -0,0 +1,649 @@ +#include + +#include "gdal.h" +#include "gdal_priv.h" +#include + +#include +#include +#include +#include +#include +#include +#include +#include // For std::ptrdiff_t + +#include "LRUCache.h" + +int MANAGED_RASTER_N_BLOCKS = pow(2, 6); +// given the pixel neighbor numbering system +// 3 2 1 +// 4 x 0 +// 5 6 7 +// These offsets are for the neighbor rows and columns +int ROW_OFFSETS[8] = {0, -1, -1, -1, 0, 1, 1, 1}; +int COL_OFFSETS[8] = {1, 1, 0, -1, -1, -1, 0, 1}; +int FLOW_DIR_REVERSE_DIRECTION[8] = {4, 5, 6, 7, 0, 1, 2, 3}; + +// if a pixel `x` has a neighbor `n` in position `i`, +// then `n`'s neighbor in position `inflow_offsets[i]` +// is the original pixel `x` +int INFLOW_OFFSETS[8] = {4, 5, 6, 7, 0, 1, 2, 3}; + +typedef std::pair BlockBufferPair; + +class NeighborTuple { +public: + int direction, x, y; + float flow_proportion; + + NeighborTuple () {} + + NeighborTuple (int direction, int x, int y, float flow_proportion) { + this->direction = direction; + this->x = x; + this->y = y; + this->flow_proportion = flow_proportion; + } +}; + + +class ManagedRaster { + public: + LRUCache* lru_cache; + std::set dirty_blocks; + int* actualBlockWidths; + int block_xsize; + int block_ysize; + int block_xmod; + int block_ymod; + int block_xbits; + int block_ybits; + long raster_x_size; + long raster_y_size; + int block_nx; + int block_ny; + char* raster_path; + int band_id; + GDALDataset* dataset; + GDALRasterBand* band; + int write_mode; + int closed; + double nodata; + + ManagedRaster() { } + + 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 ); + + raster_x_size = dataset->GetRasterXSize(); + raster_y_size = dataset->GetRasterYSize(); + + if (band_id < 1 or band_id > dataset->GetRasterCount()) { + throw std::invalid_argument( + "Error: band ID is not a valid band number. " + "This error is happening in the ManagedRaster.h extension."); + } + band = dataset->GetRasterBand(band_id); + band->GetBlockSize( &block_xsize, &block_ysize ); + + block_xmod = block_xsize - 1; + block_ymod = block_ysize - 1; + + nodata = band->GetNoDataValue(); + + if (((block_xsize & (block_xsize - 1)) != 0) or ( + (block_ysize & (block_ysize - 1)) != 0)) { + throw std::invalid_argument( + "Error: Block size is not a power of two. " + "This error is happening in the ManagedRaster.h extension."); + } + + block_xbits = log2(block_xsize); + block_ybits = log2(block_ysize); + + // integer floor division + block_nx = (raster_x_size + block_xsize - 1) / block_xsize; + block_ny = (raster_y_size + block_ysize - 1) / block_ysize; + + int actual_x = 0; + int actual_y = 0; + actualBlockWidths = (int *) CPLMalloc(sizeof(int) * block_nx * block_ny); + + for (int block_yi = 0; block_yi < block_ny; block_yi++) { + for (int block_xi = 0; block_xi < block_nx; block_xi++) { + band->GetActualBlockSize(block_xi, block_yi, &actual_x, &actual_y); + actualBlockWidths[block_yi * block_nx + block_xi] = actual_x; + } + } + + lru_cache = new LRUCache(MANAGED_RASTER_N_BLOCKS); + closed = 0; + } + + 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; + + // this is the flat index for the block + int block_index = block_yi * block_nx + block_xi; + + if (not lru_cache->exist(block_index)) { + _load_block(block_index); + } + + int idx = ((yi & block_ymod) * actualBlockWidths[block_index]) + (xi & block_xmod); + lru_cache->get(block_index)[idx] = value; + if (write_mode) { + std::set::iterator dirty_itr = dirty_blocks.find(block_index); + if (dirty_itr == dirty_blocks.end()) { + dirty_blocks.insert(block_index); + } + } + } + + 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; + + // this is the flat index for the block + int block_index = block_yi * block_nx + block_xi; + + if (not lru_cache->exist(block_index)) { + _load_block(block_index); + } + 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); + + double value = block[idx]; + return value; + } + + void _load_block(int block_index) { + int block_xi = block_index % block_nx; + int block_yi = block_index / block_nx; + + // we need the offsets to subtract from global indexes for cached array + int xoff = block_xi << block_xbits; + int yoff = block_yi << block_ybits; + + double *double_buffer; + list removed_value_list; + + // determine the block aligned xoffset for read as array + + // initially the win size is the same as the block size unless + // we're at the edge of a raster + int win_xsize = block_xsize; + int win_ysize = block_ysize; + + // load a new block + if ((xoff + win_xsize) > raster_x_size) { + win_xsize = win_xsize - (xoff + win_xsize - raster_x_size); + } + if ((yoff + win_ysize) > raster_y_size) { + win_ysize = win_ysize - (yoff + win_ysize - raster_y_size); + } + + double *pafScanline = (double *) CPLMalloc(sizeof(double) * win_xsize * win_ysize); + CPLErr err = band->RasterIO(GF_Read, xoff, yoff, win_xsize, win_ysize, + pafScanline, win_xsize, win_ysize, GDT_Float64, + 0, 0 ); + + if (err != CE_None) { + std::cerr << "Error reading block\n"; + } + lru_cache->put(block_index, pafScanline, removed_value_list); + while (not removed_value_list.empty()) { + // write the changed value back if desired + double_buffer = removed_value_list.front().second; + + if (write_mode) { + block_index = removed_value_list.front().first; + + // write back the block if it's dirty + std::set::iterator dirty_itr = dirty_blocks.find(block_index); + if (dirty_itr != dirty_blocks.end()) { + dirty_blocks.erase(dirty_itr); + + block_xi = block_index % block_nx; + block_yi = block_index / block_nx; + + xoff = block_xi << block_xbits; + yoff = block_yi << block_ybits; + + win_xsize = block_xsize; + win_ysize = block_ysize; + + if (xoff + win_xsize > raster_x_size) { + win_xsize = win_xsize - (xoff + win_xsize - raster_x_size); + } + if (yoff + win_ysize > raster_y_size) { + win_ysize = win_ysize - (yoff + win_ysize - raster_y_size); + } + err = band->RasterIO( GF_Write, xoff, yoff, win_xsize, win_ysize, + double_buffer, win_xsize, win_ysize, GDT_Float64, 0, 0 ); + if (err != CE_None) { + std::cerr << "Error writing block\n"; + } + } + } + + CPLFree(double_buffer); + removed_value_list.pop_front(); + } + } + + 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; + } + closed = 1; + + double *double_buffer; + int block_xi; + int block_yi; + int block_index; + // initially the win size is the same as the block size unless + // we're at the edge of a raster + int win_xsize; + int win_ysize; + + // we need the offsets to subtract from global indexes for cached array + int xoff; + int yoff; + + if (not write_mode) { + for (auto it = lru_cache->begin(); it != lru_cache->end(); it++) { + // write the changed value back if desired + CPLFree(it->second); + } + return; + } + + // if we get here, we're in write_mode + std::set::iterator dirty_itr; + for (auto it = lru_cache->begin(); it != lru_cache->end(); it++) { + double_buffer = it->second; + block_index = it->first; + + // write to disk if block is dirty + dirty_itr = dirty_blocks.find(block_index); + if (dirty_itr != dirty_blocks.end()) { + dirty_blocks.erase(dirty_itr); + block_xi = block_index % block_nx; + block_yi = block_index / block_nx; + + // we need the offsets to subtract from global indexes for + // cached array + xoff = block_xi << block_xbits; + yoff = block_yi << block_ybits; + + win_xsize = block_xsize; + win_ysize = block_ysize; + + // clip window sizes if necessary + if (xoff + win_xsize > raster_x_size) { + win_xsize = win_xsize - (xoff + win_xsize - raster_x_size); + } + if (yoff + win_ysize > raster_y_size) { + win_ysize = win_ysize - (yoff + win_ysize - raster_y_size); + } + CPLErr err = band->RasterIO( GF_Write, xoff, yoff, win_xsize, win_ysize, + double_buffer, win_xsize, win_ysize, GDT_Float64, 0, 0 ); + if (err != CE_None) { + std::cerr << "Error writing block\n"; + } + } + CPLFree(double_buffer); + } + GDALClose( (GDALDatasetH) dataset ); + delete lru_cache; + free(actualBlockWidths); + } +}; + +class ManagedFlowDirRaster: public ManagedRaster { + +public: + + ManagedFlowDirRaster() {} + + ManagedFlowDirRaster(char* raster_path, int band_id, bool write_mode) + : ManagedRaster(raster_path, band_id, write_mode) {} + + 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; + long xj, yj; + float flow_ji; + + 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 = get(xj, yj); + flow_ji = (0xF & (flow_dir_j >> (4 * FLOW_DIR_REVERSE_DIRECTION[n_dir]))); + + if (flow_ji) { + return false; + } + } + return true; + + } + +}; + +struct Pixel { + ManagedFlowDirRaster raster; + int x; + int y; + int val; + + Pixel() {} + + Pixel(ManagedFlowDirRaster raster, int x, int y) : raster(raster), x(x), y(y) { + double v = raster.get(x, y); + val = static_cast(v); + } +}; + +static inline NeighborTuple endVal = NeighborTuple(8, -1, -1, -1); + + +struct NeighborIterator { + 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& + + Pixel pixel; + pointer m_ptr = nullptr; + int i = 0; + + NeighborIterator() {} + NeighborIterator(NeighborTuple* n) { + m_ptr = n; + } + NeighborIterator(Pixel pixel) : pixel(pixel) { + next(); + } + + reference operator*() const { return *m_ptr; } + pointer operator->() { return m_ptr; } + + // Prefix increment + NeighborIterator& operator++() { this->next(); return *this; } + + // Postfix increment + 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; + }; + friend bool operator!= (const NeighborIterator& a, const NeighborIterator& b) { + return a.m_ptr != b.m_ptr; + }; + + virtual void next() { + long xj, yj, flow; + if (i == 8) { + m_ptr = &endVal; + return; + } + 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); + i++; + } +}; + +class DownslopeNeighborIterator: public NeighborIterator { +public: + DownslopeNeighborIterator(): NeighborIterator() {} + DownslopeNeighborIterator(NeighborTuple* n): NeighborIterator(n) {} + DownslopeNeighborIterator(Pixel p) { + pixel = p; + next(); + } + + void next() { + long xj, yj, flow; + delete m_ptr; + m_ptr = nullptr; + if (i == 8) { + 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++; + next(); + return; + } + flow = (pixel.val >> (i * 4)) & 0xF; + if (flow) { + m_ptr = new NeighborTuple(i, xj, yj, flow); + i++; + return; + } else { + i++; + next(); + } + } +}; + +class DownslopeNeighborNoSkipIterator: public NeighborIterator { +public: + DownslopeNeighborNoSkipIterator(): NeighborIterator() {} + DownslopeNeighborNoSkipIterator(NeighborTuple* n): NeighborIterator(n) {} + DownslopeNeighborNoSkipIterator(Pixel p) { + pixel = p; + next(); + } + + void next() { + long xj, yj, flow; + delete m_ptr; + m_ptr = nullptr; + if (i == 8) { + m_ptr = &endVal; + return; + } + xj = pixel.x + COL_OFFSETS[i]; + yj = pixel.y + ROW_OFFSETS[i]; + flow = (pixel.val >> (i * 4)) & 0xF; + if (flow) { + m_ptr = new NeighborTuple(i, xj, yj, flow); + i++; + return; + } else { + i++; + next(); + } + } +}; + +class UpslopeNeighborIterator: public NeighborIterator { +public: + UpslopeNeighborIterator(): NeighborIterator() {} + UpslopeNeighborIterator(NeighborTuple* n): NeighborIterator(n) {} + UpslopeNeighborIterator(Pixel p) { + pixel = p; + next(); + } + + 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; + 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++; + next(); + return; + } + flow_dir_j = pixel.raster.get(xj, yj); + flow_ji = (0xF & (flow_dir_j >> (4 * FLOW_DIR_REVERSE_DIRECTION[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++; + return; + } else { + i++; + next(); + } + } +}; + +class UpslopeNeighborNoDivideIterator: public NeighborIterator { +public: + UpslopeNeighborNoDivideIterator(): NeighborIterator() {} + UpslopeNeighborNoDivideIterator(NeighborTuple* n): NeighborIterator(n) {} + UpslopeNeighborNoDivideIterator(Pixel p) { + pixel = p; + next(); + } + + 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; + 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++; + next(); + return; + } + flow_dir_j = pixel.raster.get(xj, yj); + flow_ji = (0xF & (flow_dir_j >> (4 * FLOW_DIR_REVERSE_DIRECTION[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++; + return; + } else { + i++; + next(); + } + } +}; + +class Neighbors { +public: + Pixel pixel; + Neighbors() {} + Neighbors(Pixel pixel): pixel(pixel) {} + NeighborIterator begin() { return NeighborIterator(pixel); } + NeighborIterator end() { return NeighborIterator(&endVal); } +}; + +class DownslopeNeighbors: public Neighbors { +public: + using Neighbors::Neighbors; + DownslopeNeighborIterator begin() { return DownslopeNeighborIterator(pixel); } + DownslopeNeighborIterator end() { return DownslopeNeighborIterator(&endVal); } +}; + +class DownslopeNeighborsNoSkip: public Neighbors { +public: + using Neighbors::Neighbors; + DownslopeNeighborNoSkipIterator begin() { return DownslopeNeighborNoSkipIterator(pixel); } + DownslopeNeighborNoSkipIterator end() { return DownslopeNeighborNoSkipIterator(&endVal); } +}; + +class UpslopeNeighbors: public Neighbors { +public: + using Neighbors::Neighbors; + UpslopeNeighborIterator begin() { return UpslopeNeighborIterator(pixel); } + UpslopeNeighborIterator end() { return UpslopeNeighborIterator(&endVal); } +}; + +class UpslopeNeighborsNoDivide: public Neighbors { +public: + using Neighbors::Neighbors; + UpslopeNeighborNoDivideIterator begin() { return UpslopeNeighborNoDivideIterator(pixel); } + UpslopeNeighborNoDivideIterator end() { return UpslopeNeighborNoDivideIterator(&endVal); } +}; + +inline bool is_close(double x, double y) { + if (isnan(x) and isnan(y)) { + return true; + } + return abs(x - y) <= (pow(10, -8) + pow(10, -05) * abs(y)); +} diff --git a/src/natcap/invest/managed_raster/managed_raster.pxd b/src/natcap/invest/managed_raster/managed_raster.pxd index b3c7d0b4bb..f2d7b6d370 100644 --- a/src/natcap/invest/managed_raster/managed_raster.pxd +++ b/src/natcap/invest/managed_raster/managed_raster.pxd @@ -3,19 +3,12 @@ from libcpp.list cimport list as clist from libcpp.pair cimport pair from libcpp.set cimport set as cset +from libcpp.string cimport string from libcpp.vector cimport vector from libc.math cimport isnan -cdef struct s_neighborTuple: - int direction - int x - int y - float flow_proportion - -ctypedef s_neighborTuple NeighborTuple - # this is a least recently used cache written in C++ in an external file, -# exposing here so _ManagedRaster can use it +# exposing here so ManagedRaster can use it cdef extern from "LRUCache.h" nogil: cdef cppclass LRUCache[KEY_T, VAL_T]: LRUCache(int) @@ -25,48 +18,162 @@ cdef extern from "LRUCache.h" nogil: bint exist(KEY_T &) VAL_T get(KEY_T &) -cdef class _ManagedRaster: - cdef LRUCache[int, double*]* lru_cache - cdef cset[int] dirty_blocks - cdef int block_xsize - cdef int block_ysize - cdef int block_xmod - cdef int block_ymod - cdef int block_xbits - cdef int block_ybits - cdef long raster_x_size - cdef long raster_y_size - cdef int block_nx - cdef int block_ny - cdef int write_mode - cdef bytes raster_path - cdef int band_id - cdef int closed - - cdef inline void set(_ManagedRaster self, long xi, long yi, double value) - cdef inline double get(_ManagedRaster self, long xi, long yi) - cdef void _load_block(_ManagedRaster self, int block_index) except * - - -cdef class ManagedFlowDirRaster(_ManagedRaster): - - cdef bint is_local_high_point(ManagedFlowDirRaster self, long xi, long yi) - - cdef vector[NeighborTuple] get_upslope_neighbors(ManagedFlowDirRaster self, long xi, long yi) - - cdef vector[NeighborTuple] get_downslope_neighbors(ManagedFlowDirRaster self, long xi, long yi, bint skip_oob=*) - - -# These offsets are for the neighbor rows and columns according to the -# ordering: 3 2 1 -# 4 x 0 -# 5 6 7 -cdef int *ROW_OFFSETS -cdef int *COL_OFFSETS -cdef int *FLOW_DIR_REVERSE_DIRECTION -cdef int *INFLOW_OFFSETS - -cdef inline int is_close(double x, double y): - if isnan(x) and isnan(y): - return 1 - return abs(x - y) <= (1e-8 + 1e-05 * abs(y)) +cdef extern from "ManagedRaster.h": + cdef cppclass ManagedRaster: + LRUCache[int, double*]* lru_cache + cset[int] dirty_blocks + int block_xsize + int block_ysize + int block_xmod + int block_ymod + int block_xbits + int block_ybits + long raster_x_size + long raster_y_size + int block_nx + int block_ny + int write_mode + string raster_path + int band_id + int closed + + ManagedRaster() except + + ManagedRaster(char*, int, bool) except + + void set(long xi, long yi, double value) + double get(long xi, long yi) + void _load_block(int block_index) except * + void close() + + cdef cppclass ManagedFlowDirRaster: + LRUCache[int, double*]* lru_cache + cset[int] dirty_blocks + int block_xsize + int block_ysize + int block_xmod + int block_ymod + int block_xbits + int block_ybits + long raster_x_size + long raster_y_size + int block_nx + int block_ny + int write_mode + string raster_path + int band_id + int closed + + bint is_local_high_point(int xi, int yi) + + ManagedFlowDirRaster() except + + ManagedFlowDirRaster(char*, int, bool) except + + void set(long xi, long yi, double value) + double get(long xi, long yi) + void close() + + cdef cppclass NeighborTuple: + NeighborTuple() except + + NeighborTuple(int, int, int, float) except + + int direction, x, y + float flow_proportion + + cdef cppclass UpslopeNeighborIteratorSkip: + ManagedFlowDirRaster raster + int col + int row + int n_dir + int flow_dir + + 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: + DownslopeNeighborNoSkipIterator() + DownslopeNeighborNoSkipIterator(NeighborTuple* n) + DownslopeNeighborNoSkipIterator(Pixel) + NeighborTuple operator*() + DownslopeNeighborNoSkipIterator operator++() + bint operator==(DownslopeNeighborNoSkipIterator) + bint operator!=(DownslopeNeighborNoSkipIterator) + + cdef cppclass UpslopeNeighborIterator: + UpslopeNeighborIterator() + UpslopeNeighborIterator(NeighborTuple* n) + UpslopeNeighborIterator(Pixel) + NeighborTuple operator*() + UpslopeNeighborIterator operator++() + bint operator==(UpslopeNeighborIterator) + bint operator!=(UpslopeNeighborIterator) + + cdef cppclass UpslopeNeighborNoDivideIterator: + 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() + + bint is_close(double, double) + + int[8] INFLOW_OFFSETS + int[8] COL_OFFSETS + int[8] ROW_OFFSETS + int[8] FLOW_DIR_REVERSE_DIRECTION diff --git a/src/natcap/invest/managed_raster/managed_raster.pyx b/src/natcap/invest/managed_raster/managed_raster.pyx deleted file mode 100644 index 15f506ad3c..0000000000 --- a/src/natcap/invest/managed_raster/managed_raster.pyx +++ /dev/null @@ -1,430 +0,0 @@ -# cython: language_level=3 -# distutils: language = c++ -import os - -import numpy -import pygeoprocessing -cimport numpy -cimport cython -from osgeo import gdal - -from cpython.mem cimport PyMem_Malloc, PyMem_Free -from cython.operator cimport dereference as deref -from cython.operator cimport preincrement as inc -from libc.time cimport time as ctime -from libcpp.pair cimport pair -from libcpp.set cimport set as cset -from libcpp.list cimport list as clist -from libcpp.stack cimport stack -from libcpp.vector cimport vector - -# this ctype is used to store the block ID and the block buffer as one object -# inside Managed Raster -ctypedef pair[int, double*] BlockBufferPair -# Number of raster blocks to hold in memory at once per Managed Raster -cdef int MANAGED_RASTER_N_BLOCKS = 2**6 - -# given the pixel neighbor numbering system -# 3 2 1 -# 4 x 0 -# 5 6 7 -# These offsets are for the neighbor rows and columns -cdef int *ROW_OFFSETS = [0, -1, -1, -1, 0, 1, 1, 1] -cdef int *COL_OFFSETS = [1, 1, 0, -1, -1, -1, 0, 1] -cdef int *FLOW_DIR_REVERSE_DIRECTION = [4, 5, 6, 7, 0, 1, 2, 3] - -# if a pixel `x` has a neighbor `n` in position `i`, -# then `n`'s neighbor in position `inflow_offsets[i]` -# is the original pixel `x` -cdef int *INFLOW_OFFSETS = [4, 5, 6, 7, 0, 1, 2, 3] - -# a class to allow fast random per-pixel access to a raster for both setting -# and reading pixels. Copied from src/pygeoprocessing/routing/routing.pyx, -# revision 891288683889237cfd3a3d0a1f09483c23489fca. -cdef class _ManagedRaster: - - def __cinit__(self, raster_path, band_id, 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. - """ - raster_info = pygeoprocessing.get_raster_info(raster_path) - self.raster_x_size, self.raster_y_size = raster_info['raster_size'] - self.block_xsize, self.block_ysize = raster_info['block_size'] - self.block_xmod = self.block_xsize-1 - self.block_ymod = self.block_ysize-1 - - if not (1 <= band_id <= raster_info['n_bands']): - err_msg = ( - "Error: band ID (%s) is not a valid band number. " - "This exception is happening in Cython, so it will cause a " - "hard seg-fault, but it's otherwise meant to be a " - "ValueError." % (band_id)) - print(err_msg) - raise ValueError(err_msg) - self.band_id = band_id - - if (self.block_xsize & (self.block_xsize - 1) != 0) or ( - self.block_ysize & (self.block_ysize - 1) != 0): - # If inputs are not a power of two, this will at least print - # an error message. Unfortunately with Cython, the exception will - # present itself as a hard seg-fault, but I'm leaving the - # ValueError in here at least for readability. - err_msg = ( - "Error: Block size is not a power of two: " - "block_xsize: %d, %d, %s. This exception is happening" - "in Cython, so it will cause a hard seg-fault, but it's" - "otherwise meant to be a ValueError." % ( - self.block_xsize, self.block_ysize, raster_path)) - print(err_msg) - raise ValueError(err_msg) - - self.block_xbits = numpy.log2(self.block_xsize) - self.block_ybits = numpy.log2(self.block_ysize) - self.block_nx = ( - self.raster_x_size + (self.block_xsize) - 1) // self.block_xsize - self.block_ny = ( - self.raster_y_size + (self.block_ysize) - 1) // self.block_ysize - - self.lru_cache = new LRUCache[int, double*](MANAGED_RASTER_N_BLOCKS) - self.raster_path = raster_path - self.write_mode = write_mode - self.closed = 0 - - def __dealloc__(self): - """Deallocate _ManagedRaster. - - This operation manually frees memory from the LRUCache and writes any - dirty memory blocks back to the raster if `self.write_mode` is True. - """ - self.close() - - def close(self): - """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 self.closed: - return - self.closed = 1 - cdef int xi_copy, yi_copy - cdef numpy.ndarray[double, ndim=2] block_array = numpy.empty( - (self.block_ysize, self.block_xsize)) - cdef double *double_buffer - cdef int block_xi - cdef int block_yi - # initially the win size is the same as the block size unless - # we're at the edge of a raster - cdef int win_xsize - cdef int win_ysize - - # we need the offsets to subtract from global indexes for cached array - cdef int xoff - cdef int yoff - - cdef clist[BlockBufferPair].iterator it = self.lru_cache.begin() - cdef clist[BlockBufferPair].iterator end = self.lru_cache.end() - if not self.write_mode: - while it != end: - # write the changed value back if desired - PyMem_Free(deref(it).second) - inc(it) - return - - raster = gdal.OpenEx( - self.raster_path, gdal.GA_Update | gdal.OF_RASTER) - raster_band = raster.GetRasterBand(self.band_id) - - # if we get here, we're in write_mode - cdef cset[int].iterator dirty_itr - while it != end: - double_buffer = deref(it).second - block_index = deref(it).first - - # write to disk if block is dirty - dirty_itr = self.dirty_blocks.find(block_index) - if dirty_itr != self.dirty_blocks.end(): - self.dirty_blocks.erase(dirty_itr) - block_xi = block_index % self.block_nx - block_yi = block_index / self.block_nx - - # we need the offsets to subtract from global indexes for - # cached array - xoff = block_xi << self.block_xbits - yoff = block_yi << self.block_ybits - - win_xsize = self.block_xsize - win_ysize = self.block_ysize - - # clip window sizes if necessary - if xoff+win_xsize > self.raster_x_size: - win_xsize = win_xsize - ( - xoff+win_xsize - self.raster_x_size) - if yoff+win_ysize > self.raster_y_size: - win_ysize = win_ysize - ( - yoff+win_ysize - self.raster_y_size) - - for xi_copy in xrange(win_xsize): - for yi_copy in xrange(win_ysize): - block_array[yi_copy, xi_copy] = ( - double_buffer[ - (yi_copy << self.block_xbits) + xi_copy]) - raster_band.WriteArray( - block_array[0:win_ysize, 0:win_xsize], - xoff=xoff, yoff=yoff) - PyMem_Free(double_buffer) - inc(it) - raster_band.FlushCache() - raster_band = None - raster = None - - cdef inline void set(self, long xi, long yi, double value): - """Set the pixel at `xi,yi` to `value`.""" - cdef int block_xi = xi >> self.block_xbits - cdef int block_yi = yi >> self.block_ybits - # this is the flat index for the block - cdef int block_index = block_yi * self.block_nx + block_xi - if not self.lru_cache.exist(block_index): - self._load_block(block_index) - self.lru_cache.get( - block_index)[ - ((yi & (self.block_ymod))<> self.block_xbits - cdef int block_yi = yi >> self.block_ybits - # this is the flat index for the block - cdef int block_index = block_yi * self.block_nx + block_xi - if not self.lru_cache.exist(block_index): - self._load_block(block_index) - return self.lru_cache.get( - block_index)[ - ((yi & (self.block_ymod))< self.raster_x_size: - win_xsize = win_xsize - (xoff+win_xsize - self.raster_x_size) - if yoff+win_ysize > self.raster_y_size: - win_ysize = win_ysize - (yoff+win_ysize - self.raster_y_size) - - raster = gdal.OpenEx(self.raster_path, gdal.OF_RASTER) - raster_band = raster.GetRasterBand(self.band_id) - block_array = raster_band.ReadAsArray( - xoff=xoff, yoff=yoff, win_xsize=win_xsize, - win_ysize=win_ysize).astype( - numpy.float64) - raster_band = None - raster = None - double_buffer = PyMem_Malloc( - (sizeof(double) << self.block_xbits) * win_ysize) - for xi_copy in xrange(win_xsize): - for yi_copy in xrange(win_ysize): - double_buffer[(yi_copy<block_index, double_buffer, removed_value_list) - - if self.write_mode: - raster = gdal.OpenEx( - self.raster_path, gdal.GA_Update | gdal.OF_RASTER) - raster_band = raster.GetRasterBand(self.band_id) - - block_array = numpy.empty( - (self.block_ysize, self.block_xsize), dtype=numpy.double) - while not removed_value_list.empty(): - # write the changed value back if desired - double_buffer = removed_value_list.front().second - - if self.write_mode: - block_index = removed_value_list.front().first - - # write back the block if it's dirty - dirty_itr = self.dirty_blocks.find(block_index) - if dirty_itr != self.dirty_blocks.end(): - self.dirty_blocks.erase(dirty_itr) - - block_xi = block_index % self.block_nx - block_yi = block_index // self.block_nx - - xoff = block_xi << self.block_xbits - yoff = block_yi << self.block_ybits - - win_xsize = self.block_xsize - win_ysize = self.block_ysize - - if xoff+win_xsize > self.raster_x_size: - win_xsize = win_xsize - ( - xoff+win_xsize - self.raster_x_size) - if yoff+win_ysize > self.raster_y_size: - win_ysize = win_ysize - ( - yoff+win_ysize - self.raster_y_size) - - for xi_copy in xrange(win_xsize): - for yi_copy in xrange(win_ysize): - block_array[yi_copy, xi_copy] = double_buffer[ - (yi_copy << self.block_xbits) + xi_copy] - raster_band.WriteArray( - block_array[0:win_ysize, 0:win_xsize], - xoff=xoff, yoff=yoff) - PyMem_Free(double_buffer) - removed_value_list.pop_front() - - if self.write_mode: - raster_band = None - raster = None - - -cdef class ManagedFlowDirRaster(_ManagedRaster): - - cdef bint is_local_high_point(self, long xi, long 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. - """ - return self.get_upslope_neighbors(xi, yi).size() == 0 - - @cython.cdivision(True) - cdef vector[NeighborTuple] get_upslope_neighbors( - ManagedFlowDirRaster self, long xi, long yi): - """Return upslope neighbors of a given pixel. - - 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: - libcpp.vector of NeighborTuples. Each NeighborTuple has - the attributes ``direction`` (integer flow direction 0-7 - of the neighbor relative to the original pixel), ``x`` - and ``y`` (integer coordinates of the neighbor in pixel - space), and ``flow_proportion`` (fraction of the flow - from the neighbor that flows to the original pixel). - """ - cdef int n_dir, flow_dir_j, idx - cdef long xj, yj - cdef float flow_ji, flow_dir_j_sum - - cdef NeighborTuple n - cdef vector[NeighborTuple] upslope_neighbor_tuples - - for n_dir in range(8): - xj = xi + COL_OFFSETS[n_dir] - yj = yi + ROW_OFFSETS[n_dir] - if (xj < 0 or xj >= self.raster_x_size or - yj < 0 or yj >= self.raster_y_size): - continue - flow_dir_j = self.get(xj, yj) - flow_ji = (0xF & (flow_dir_j >> (4 * FLOW_DIR_REVERSE_DIRECTION[n_dir]))) - - if flow_ji: - flow_dir_j_sum = 0 - for idx in range(8): - flow_dir_j_sum += (flow_dir_j >> (idx * 4)) & 0xF - n.direction = n_dir - n.x = xj - n.y = yj - n.flow_proportion = flow_ji / flow_dir_j_sum - upslope_neighbor_tuples.push_back(n) - - return upslope_neighbor_tuples - - @cython.cdivision(True) - cdef vector[NeighborTuple] get_downslope_neighbors( - ManagedFlowDirRaster self, long xi, long yi, bint skip_oob=True): - """Return downslope neighbors of a given pixel. - - 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 - skip_oob (bool): if True, do not return neighbors that fall - outside the raster bounds. - - Returns: - libcpp.vector of NeighborTuples. Each NeighborTuple has - the attributes ``direction`` (integer flow direction 0-7 - of the neighbor relative to the original pixel), ``x`` - and ``y`` (integer coordinates of the neighbor in pixel - space), and ``flow_proportion`` (fraction of the flow - from the neighbor that flows to the original pixel). - """ - cdef int n_dir - cdef long xj, yj - cdef float flow_ij - - cdef NeighborTuple n - cdef vector[NeighborTuple] downslope_neighbor_tuples - - cdef int flow_dir = self.get(xi, yi) - cdef float flow_sum = 0 - - cdef int i = 0 - for n_dir in range(8): - # flows in this direction - xj = xi + COL_OFFSETS[n_dir] - yj = yi + ROW_OFFSETS[n_dir] - if skip_oob and (xj < 0 or xj >= self.raster_x_size or - yj < 0 or yj >= self.raster_y_size): - continue - flow_ij = (flow_dir >> (n_dir * 4)) & 0xF - flow_sum += flow_ij - if flow_ij: - n = NeighborTuple() - n.direction = n_dir - n.x = xj - n.y = yj - n.flow_proportion = flow_ij - downslope_neighbor_tuples.push_back(n) - i += 1 - - for j in range(i): - downslope_neighbor_tuples[j].flow_proportion = ( - downslope_neighbor_tuples[j].flow_proportion / flow_sum) - - return downslope_neighbor_tuples diff --git a/src/natcap/invest/ndr/ndr_core.pyx b/src/natcap/invest/ndr/ndr_core.pyx index f2785c3581..2f67113bed 100644 --- a/src/natcap/invest/ndr/ndr_core.pyx +++ b/src/natcap/invest/ndr/ndr_core.pyx @@ -11,10 +11,18 @@ 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 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 cdef extern from "time.h" nogil: ctypedef int time_t @@ -68,19 +76,21 @@ def ndr_eff_calculation( # 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, 1, False) - cdef _ManagedRaster crit_len_raster = _ManagedRaster( - crit_len_path, 1, False) + 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 _ManagedRaster retention_eff_lulc_raster = _ManagedRaster( - retention_eff_lulc_path, 1, False) cdef float retention_eff_nodata = pygeoprocessing.get_raster_info( retention_eff_lulc_path)['nodata'][0] - cdef _ManagedRaster effective_retention_raster = _ManagedRaster( - effective_retention_path, 1, True) - cdef ManagedFlowDirRaster mfd_flow_direction_raster = ManagedFlowDirRaster( - mfd_flow_direction_path, 1, False) # create direction raster in bytes def _mfd_to_flow_dir_op(mfd_array): @@ -97,18 +107,19 @@ def ndr_eff_calculation( [(mfd_flow_direction_path, 1)], _mfd_to_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, 1, True) + 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 + 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): @@ -126,15 +137,25 @@ def ndr_eff_calculation( should_seed = 0 # see if this pixel drains to nodata or the edge, if so it's # a drain - for neighbor in ( - mfd_flow_direction_raster.get_downslope_neighbors( - global_col, global_row, skip_oob=False)): - if (neighbor.x < 0 or neighbor.x >= n_cols or - neighbor.y < 0 or neighbor.y >= n_rows or - to_process_flow_directions_raster.get( - neighbor.x, neighbor.y) == 0): - should_seed = 1 - outflow_dirs &= ~(1 << neighbor.direction) + 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 @@ -156,7 +177,7 @@ def ndr_eff_calculation( 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. + # 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 @@ -166,14 +187,16 @@ def ndr_eff_calculation( 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 - for neighbor in mfd_flow_direction_raster.get_downslope_neighbors( - global_col, global_row, skip_oob=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: @@ -208,14 +231,16 @@ def ndr_eff_calculation( 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): - for neighbor in mfd_flow_direction_raster.get_upslope_neighbors( - global_col, global_row): + 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 = ( @@ -236,5 +261,10 @@ def ndr_eff_calculation( # 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() - os.remove(to_process_flow_directions_path) diff --git a/src/natcap/invest/scenic_quality/viewshed.pyx b/src/natcap/invest/scenic_quality/viewshed.pyx index 0d7404e7aa..523700913c 100644 --- a/src/natcap/invest/scenic_quality/viewshed.pyx +++ b/src/natcap/invest/scenic_quality/viewshed.pyx @@ -39,7 +39,7 @@ from libcpp.queue cimport queue from libc cimport math cimport numpy cimport cython -from ..managed_raster.managed_raster cimport _ManagedRaster +from ..managed_raster.managed_raster cimport ManagedRaster LOGGER = logging.getLogger(__name__) LOGGER.addHandler(logging.NullHandler()) @@ -315,8 +315,8 @@ def viewshed(dem_raster_path_band, (pixel_xsize, pixel_ysize)) # Verify that the block sizes are powers of 2 and are square. - # This is needed for the _ManagedRaster classes. If this is not asserted - # here, the _ManagedRaster classes will crash with a segfault. + # This is needed for the ManagedRaster classes. If this is not asserted + # here, the ManagedRaster classes will crash with a segfault. block_xsize, block_ysize = dem_raster_info['block_size'] if (block_xsize & (block_xsize - 1) != 0 or ( block_ysize & (block_ysize - 1) != 0)) or ( @@ -350,12 +350,12 @@ def viewshed(dem_raster_path_band, raster_driver_creation_tuple=BYTE_GTIFF_CREATION_OPTIONS) # LRU-cached rasters for easier access to individual pixels. - cdef _ManagedRaster dem_managed_raster = ( - _ManagedRaster(dem_raster_path_band[0], dem_raster_path_band[1], 0)) - cdef _ManagedRaster aux_managed_raster = ( - _ManagedRaster(aux_filepath, 1, 1)) - cdef _ManagedRaster visibility_managed_raster = ( - _ManagedRaster(visibility_filepath, 1, 1)) + cdef ManagedRaster dem_managed_raster = ( + ManagedRaster(dem_raster_path_band[0].encode('utf-8'), dem_raster_path_band[1], 0)) + cdef ManagedRaster aux_managed_raster = ( + ManagedRaster(aux_filepath.encode('utf-8'), 1, 1)) + cdef ManagedRaster visibility_managed_raster = ( + ManagedRaster(visibility_filepath.encode('utf-8'), 1, 1)) # get the pixel size in terms of meters. dem_srs = osr.SpatialReference() diff --git a/src/natcap/invest/sdr/sdr_core.pyx b/src/natcap/invest/sdr/sdr_core.pyx index afe1925194..4cfb961a7b 100644 --- a/src/natcap/invest/sdr/sdr_core.pyx +++ b/src/natcap/invest/sdr/sdr_core.pyx @@ -1,17 +1,19 @@ import logging import os -import numpy import pygeoprocessing -cimport numpy 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 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 cdef extern from "time.h" nogil: @@ -84,13 +86,13 @@ def calculate_sediment_deposition( gdal.GDT_Float32, [target_nodata]) cdef ManagedFlowDirRaster mfd_flow_direction_raster = ManagedFlowDirRaster( - mfd_flow_direction_path, 1, False) - cdef _ManagedRaster e_prime_raster = _ManagedRaster( - e_prime_path, 1, False) - cdef _ManagedRaster sdr_raster = _ManagedRaster(sdr_path, 1, False) - cdef _ManagedRaster f_raster = _ManagedRaster(f_path, 1, True) - cdef _ManagedRaster sediment_deposition_raster = _ManagedRaster( - target_sediment_deposition_path, 1, True) + 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) @@ -103,7 +105,7 @@ def calculate_sediment_deposition( 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 long xs, ys + cdef int xs, ys cdef long flat_index cdef long seed_col = 0 cdef long seed_row = 0 @@ -116,6 +118,8 @@ def calculate_sediment_deposition( 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): @@ -134,6 +138,10 @@ def calculate_sediment_deposition( 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 @@ -154,14 +162,12 @@ def calculate_sediment_deposition( # the weighted sum of flux flowing onto this pixel from # all neighbors f_j_weighted_sum = 0 - for neighbor in ( - mfd_flow_direction_raster.get_upslope_neighbors( - global_col, global_row)): - + 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 @@ -171,9 +177,12 @@ def calculate_sediment_deposition( # neighbor # (sum over k ∈ K of SDR_k * p(i,k) in the equation above) downslope_sdr_weighted_sum = 0 - for neighbor in ( - mfd_flow_direction_raster.get_downslope_neighbors( - global_col, global_row)): + 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 @@ -185,7 +194,6 @@ def calculate_sediment_deposition( 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 @@ -194,18 +202,14 @@ def calculate_sediment_deposition( # completed upslope_neighbors_processed = 1 # iterate over each neighbor-of-neighbor - for neighbor_of_neighbor in ( - mfd_flow_direction_raster.get_upslope_neighbors( - neighbor.x, neighbor.y)): - # no need to push the one we're currently - # calculating back onto the stack - if (INFLOW_OFFSETS[neighbor_of_neighbor.direction] == - neighbor.direction): + 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): + 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 @@ -222,6 +226,9 @@ def calculate_sediment_deposition( 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 @@ -257,5 +264,10 @@ def calculate_sediment_deposition( f_raster.set(global_col, global_row, f_i) n_pixels_processed += win_xsize * win_ysize - LOGGER.info('Sediment deposition 100% complete') 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') 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 da93418534..3e37c25e5f 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 @@ -15,9 +15,15 @@ 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 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: @@ -77,7 +83,7 @@ cpdef calculate_local_recharge( 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, p_ij_base + 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 @@ -86,7 +92,7 @@ cpdef calculate_local_recharge( 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 ManagedRaster et0_m_raster, qf_m_raster, kc_m_raster cdef numpy.ndarray[numpy.npy_float32, ndim=1] alpha_month_array = ( numpy.array( @@ -101,42 +107,44 @@ cpdef calculate_local_recharge( 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, 1, 0) + 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 - et0_m_raster_list = [] + cdef vector[ManagedRaster] et0_m_rasters et0_m_nodata_list = [] for et0_path in et0_path_list: - et0_m_raster_list.append(_ManagedRaster(et0_path, 1, 0)) + 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) - precip_m_raster_list = [] + cdef vector[ManagedRaster] precip_m_rasters precip_m_nodata_list = [] for precip_m_path in precip_path_list: - precip_m_raster_list.append(_ManagedRaster(precip_m_path, 1, 0)) + 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) - qf_m_raster_list = [] + cdef vector[ManagedRaster] qf_m_rasters qf_m_nodata_list = [] for qf_m_path in qf_m_path_list: - qf_m_raster_list.append(_ManagedRaster(qf_m_path, 1, 0)) + 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]) - kc_m_raster_list = [] + cdef vector[ManagedRaster] kc_m_rasters kc_m_nodata_list = [] for kc_m_path in kc_path_list: - kc_m_raster_list.append(_ManagedRaster(kc_m_path, 1, 0)) + 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]) @@ -144,33 +152,32 @@ cpdef calculate_local_recharge( 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, 1, 1) + 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, 1, 1) + 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, 1, 1) + 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, 1, 1) + 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, 1, 1) - + 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): @@ -194,8 +201,7 @@ cpdef calculate_local_recharge( 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)) + work_queue.push(pair[long, long](xs_root, ys_root)) while work_queue.size() > 0: xi = work_queue.front().first @@ -212,8 +218,10 @@ cpdef calculate_local_recharge( upslope_defined = 1 # initialize to 0 so we indicate we haven't tracked any # mfd values yet - l_sum_avail_i = 0 - for neighbor in flow_raster.get_upslope_neighbors(xi, yi): + 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) @@ -225,11 +233,14 @@ cpdef calculate_local_recharge( # 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 @@ -241,13 +252,13 @@ cpdef calculate_local_recharge( for m_index in range(12): precip_m_raster = ( - <_ManagedRaster?>precip_m_raster_list[m_index]) + precip_m_rasters[m_index]) qf_m_raster = ( - <_ManagedRaster?>qf_m_raster_list[m_index]) + qf_m_rasters[m_index]) et0_m_raster = ( - <_ManagedRaster?>et0_m_raster_list[m_index]) + et0_m_rasters[m_index]) kc_m_raster = ( - <_ManagedRaster?>kc_m_raster_list[m_index]) + kc_m_rasters[m_index]) et0_nodata = et0_m_nodata_list[m_index] precip_nodata = precip_m_nodata_list[m_index] @@ -289,9 +300,25 @@ cpdef calculate_local_recharge( target_li_raster.set(xi, yi, l_i) target_li_avail_raster.set(xi, yi, l_avail_i) - for neighbor in flow_raster.get_downslope_neighbors(xi, yi): + 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() + def route_baseflow_sum( flow_dir_mfd_path, l_path, l_avail_path, l_sum_path, @@ -323,8 +350,7 @@ def route_baseflow_sum( 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 float p_ij - cdef int flow_dir_i, p_ij_base + 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 @@ -345,16 +371,16 @@ def route_baseflow_sum( flow_dir_mfd_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, 1, 1) - cdef _ManagedRaster target_b_raster = _ManagedRaster( - target_b_path, 1, 1) - cdef _ManagedRaster l_raster = _ManagedRaster(l_path, 1, 0) - cdef _ManagedRaster l_avail_raster = _ManagedRaster(l_avail_path, 1, 0) - cdef _ManagedRaster l_sum_raster = _ManagedRaster(l_sum_path, 1, 0) + 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, 1, 0) - cdef _ManagedRaster stream_raster = _ManagedRaster(stream_path, 1, 0) + 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( @@ -376,8 +402,9 @@ def route_baseflow_sum( # search for a pixel that has no downslope neighbors, # or whose downslope neighbors all have nodata in the stream raster (?) outlet = 1 - for neighbor in flow_dir_mfd_raster.get_downslope_neighbors( - xs_root, ys_root): + 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 @@ -403,7 +430,8 @@ def route_baseflow_sum( b_sum_i = 0.0 downslope_defined = 1 - for neighbor in flow_dir_mfd_raster.get_downslope_neighbors(xi, yi): + 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 @@ -417,7 +445,7 @@ def route_baseflow_sum( 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 += p_ij * ( + b_sum_i += neighbor.flow_proportion * ( (1 - l_avail_j / l_sum_j) * ( b_sum_j / (l_sum_j - l_j))) else: @@ -439,5 +467,14 @@ def route_baseflow_sum( target_b_sum_raster.set(xi, yi, b_sum_i) current_pixel += 1 - for neighbor in flow_dir_mfd_raster.get_upslope_neighbors(xi, yi): + 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()