From 561618ed821bd488da9553dbf1054b03989068bb Mon Sep 17 00:00:00 2001 From: Richard Barnes Date: Fri, 13 May 2022 13:32:53 -0500 Subject: [PATCH] Add Depression Hierarchy to pyrichdem wrapper --- .gitignore | 1 + .../depressions/depression_hierarchy.hpp | 35 ++---- include/richdem/terrain_generation.hpp | 2 +- src/terrain_generation/terrain_generation.cpp | 17 ++- wrappers/pyrichdem/richdem/__init__.py | 108 ++++++++++++------ wrappers/pyrichdem/setup.py | 4 +- wrappers/pyrichdem/src/pywrapper.cpp | 36 +++++- wrappers/pyrichdem/src/pywrapper.hpp | 2 + wrappers/pyrichdem/tests/tests.py | 34 ++++++ 9 files changed, 168 insertions(+), 71 deletions(-) create mode 100644 wrappers/pyrichdem/tests/tests.py diff --git a/.gitignore b/.gitignore index 549084e5..7c41ae27 100644 --- a/.gitignore +++ b/.gitignore @@ -10,3 +10,4 @@ travis.log build/ .vscode/ .log +*.so diff --git a/include/richdem/depressions/depression_hierarchy.hpp b/include/richdem/depressions/depression_hierarchy.hpp index e202e55f..2d322a5f 100644 --- a/include/richdem/depressions/depression_hierarchy.hpp +++ b/include/richdem/depressions/depression_hierarchy.hpp @@ -31,14 +31,14 @@ typedef uint32_t dh_label_t; typedef uint32_t flat_c_idx; //Some special valuess -const dh_label_t NO_PARENT = std::numeric_limits::max(); -const dh_label_t NO_VALUE = std::numeric_limits::max(); +constexpr dh_label_t NO_PARENT = std::numeric_limits::max(); +constexpr dh_label_t NO_VALUE = std::numeric_limits::max(); //This class holds information about a depression. Its pit cell and outlet cell -//(in flat-index form) as well as the elevations of these cells. It also notes //what is flat-index form? +//(in flat-index form) as well as the elevations of these cells. It also notes //the depression's parent. The parent of the depression is the outlet through //which it must flow in order to reach the ocean. If a depression has more than -//one outlet at the same level one of them is arbitrarily chosen; hopefully this //so, everything should have a parent except for the ocean, right? +//one outlet at the same level one of them is arbitrarily chosen; hopefully this //happens only rarely in natural environments. template struct Depression { @@ -141,7 +141,7 @@ struct Outlet { //determine if we've already found an outlet for a depression. We'll look at //outlets from lowest to highest, so if an outlet already exists for a //depression, it is that depression's lowest outlet. - bool operator==(const Outlet &o) const { //so beyond just checking, is this somehow preventing it from being recorded if one already exists? How does this work? + bool operator==(const Outlet &o) const { //Outlets are the same if they link two depressions, regardless of the //depressions' labels storage order within this class. return depa==o.depa && depb==o.depb; @@ -164,25 +164,16 @@ struct OutletHash { }; - -//The regular mod function allows negative numbers to stay negative. This mod -//function wraps negative numbers around. For instance, if a=-1 and n=100, then -//the result is 99. -int ModFloor(int a, int n) { - return ((a % n) + n) % n; -} - - template using PriorityQueue = radix_heap::pair_radix_heap; // using PriorityQueue = GridCellZk_high_pq pq; //Cell is not part of a depression -const dh_label_t NO_DEP = std::numeric_limits::max(); +constexpr dh_label_t NO_DEP = std::numeric_limits::max(); //Cell is part of the ocean and a place from which we begin searching for //depressions. -const dh_label_t OCEAN = 0; +constexpr dh_label_t OCEAN = 0; template using DepressionHierarchy = std::vector>; @@ -252,7 +243,7 @@ std::ostream& operator<<(std::ostream &out, const DepressionHierarchy &d // //@param dem - 2D array of elevations. May be in any data format. // -//@return label - A label indiciate which depression the cell belongs to. +//@return label - A label indicating which depression the cell belongs to. // The indicated label is always the leaf of the depression // hierarchy, or the OCEAN. // @@ -370,7 +361,6 @@ DepressionHierarchy GetDepressionHierarchy( bool has_lower = false; //Pretend we have no lower neighbours for(int n=1;n<=neighbours;n++){ //Check out our neighbours //Use offset to get neighbour x coordinate, wrapping as needed - // const int nx = ModFloor(x+dx[n],dem.width()); const int nx = x+dx[n]; //Use offset to get neighbour y coordinate const int ny = y+dy[n]; @@ -476,10 +466,9 @@ DepressionHierarchy GetDepressionHierarchy( auto &newdep = depressions.emplace_back(); //Add the next flat (increases size by 1) newdep.pit_cell = dem.xyToI(cx,cy); //Make a note of the pit cell's location newdep.pit_elev = celev; //Make a note of the pit cell's elevation - newdep.dep_label = clabel; //I am storing the label in the object so that I can find it later and call up the number of cells and volume (better way of doing this?) -- I have since realised I can use the index in the depressions array. So perhaps the label is no longer needed? + newdep.dep_label = clabel; //Make a note of the depression's label //TODO: It might be possible to remove this variable entirely label(ci) = clabel; //Update cell with new label } else { - //Cell has already been assigned to a depression. In this case, one of two //things is true. (1) This cell is on the frontier of our search, in which //case the cell has neighbours which have not yet been seen. (2) This cell @@ -489,12 +478,8 @@ DepressionHierarchy GetDepressionHierarchy( //However, it is harmless to check on them again. } - //TODO: Update the appropriate depression's cell_count and dep_vol variables I did this in the else if above, and then in the if and the else below. I add a cell to the count whenever it is added to the depression and add its elevation to the total elevations. - //here. Then I calculate the total volume only when we find an outlet (Good way to test this? Print values of volumes only of those that make the outlet queue? I get some negative values sometimes so I may have done something wrong, but what if it's an 'outlet' at the highest point of the depression?) - //Consider the cell's neighbours for(int n=1;n<=neighbours;n++){ - // const int nx = ModFloor(cx+dx[n],dem.width()); //Get neighbour's x-coordinate using an offset and wrapping const int nx = cx + dx[n]; //Get neighbour's y-coordinate using an offset const int ny = cy + dy[n]; //Get neighbour's y-coordinate using an offset if(!dem.inGrid(nx,ny)) //Is this cell in the grid? @@ -720,7 +705,7 @@ DepressionHierarchy GetDepressionHierarchy( //resize! const auto depa_pitcell_temp = depa.pit_cell; - auto &newdep = depressions.emplace_back(); //is it right to create a new depression for the metadepression like this? + auto &newdep = depressions.emplace_back(); newdep.lchild = depa_set; newdep.rchild = depb_set; newdep.dep_label = newlabel; diff --git a/include/richdem/terrain_generation.hpp b/include/richdem/terrain_generation.hpp index 2ee3ce6a..1377560a 100644 --- a/include/richdem/terrain_generation.hpp +++ b/include/richdem/terrain_generation.hpp @@ -4,6 +4,6 @@ namespace richdem { -Array2D perlin(const int size, const uint32_t seed); +void generate_perlin_terrain(Array2D& arr, const uint32_t seed); } \ No newline at end of file diff --git a/src/terrain_generation/terrain_generation.cpp b/src/terrain_generation/terrain_generation.cpp index 59c6dcb3..43e03d4f 100644 --- a/src/terrain_generation/terrain_generation.cpp +++ b/src/terrain_generation/terrain_generation.cpp @@ -1,21 +1,26 @@ +#include "PerlinNoise.h" + #include #include #include -#include "PerlinNoise.h" +#include namespace richdem { -Array2D perlin(const int size, const uint32_t seed){ +void generate_perlin_terrain(Array2D& arr, const uint32_t seed){ PerlinNoise pn(seed); - Array2D terrain(size,size); + if(arr.width() != arr.height()){ + throw std::runtime_error("Perlin noise array must be square!"); + } + + const auto size = arr.width(); + for(int y=0;y(size),10*y/static_cast(size),0.8); + arr(x,y) = pn.noise(10*x/static_cast(size),10*y/static_cast(size),0.8); } - - return terrain; } } \ No newline at end of file diff --git a/wrappers/pyrichdem/richdem/__init__.py b/wrappers/pyrichdem/richdem/__init__.py index 379d02aa..fe016b2d 100644 --- a/wrappers/pyrichdem/richdem/__init__.py +++ b/wrappers/pyrichdem/richdem/__init__.py @@ -2,9 +2,10 @@ import datetime from optparse import Option import pkg_resources -from typing import Optional +from typing import Any, Final, List, Optional, Tuple, Union import numpy as np +import sys try: import _richdem @@ -12,6 +13,8 @@ print("COULD NOT LOAD RichDEM ENGINE! NOTHING WILL WORK!") raise e +from _richdem import depression_hierarchy + try: from osgeo import gdal @@ -20,6 +23,7 @@ except: GDAL_AVAILABLE = False +STANDARD_GEOTRANSFORM: Final[np.ndarray] = np.array([0, 1, 0, 0, 0, -1]) def _RichDEMVersion() -> str: return "RichDEM (Python {pyver}) (hash={hash}, hashdate={compdate})".format( @@ -29,8 +33,8 @@ def _RichDEMVersion() -> str: ) -def _AddAnalysis(rda, analysis): - if type(rda) not in [rdarray, rd3array]: +def _AddAnalysis(rda: "rdarray", analysis: str) -> None: + if type(rda) not in {rdarray, rd3array}: raise Exception("An rdarray or rd3array is required!") metastr = "\n{nowdate} | {verstr} | {analysis}".format( @@ -47,24 +51,24 @@ def _AddAnalysis(rda, analysis): def rdShow( - rda, + rda: "rdarray", ignore_colours=[], - show=True, - axes=True, - cmap="gray", - vmin=None, - vmax=None, - xmin=None, - xmax=None, - ymin=None, - ymax=None, - zxmin=None, - zxmax=None, - zymin=None, - zymax=None, - figsize=(4, 4), - zcolor="red", - zloc=1, + show: bool = True, + axes: bool = True, + cmap: str = "gray", + vmin: Optional[int] = None, + vmax: Optional[int] = None, + xmin: Optional[int] = None, + xmax: Optional[int] = None, + ymin: Optional[int] = None, + ymax: Optional[int] = None, + zxmin: Optional[int] = None, + zxmax: Optional[int] = None, + zymin: Optional[int] = None, + zymax: Optional[int] = None, + figsize: Tuple[int, int] = (4, 4), + zcolor: str = "red", + zloc: int = 1, ): if type(rda) is np.ndarray: rda = rdarray(rda) @@ -148,8 +152,8 @@ def rdShow( class rdarray(np.ndarray): def __new__( - cls, array, meta_obj=None, no_data=None, dtype=None, order=None, **kwargs - ): + cls, array, meta_obj=None, no_data: Optional[Union[float, int]]=None, dtype=None, order=None, **kwargs: Any + ) -> "rdarray": obj = np.asarray(array, dtype=dtype, order=order).view(cls) if meta_obj is not None: @@ -189,7 +193,7 @@ def wrap(self): } dtype = str(self.dtype) if dtype not in richdem_arrs: - raise Exception("No equivalent RichDEM datatype.") + raise Exception(f"No equivalent RichDEM datatype to '{dtype}'.") rda = richdem_arrs[dtype](self) @@ -199,7 +203,7 @@ def wrap(self): else: rda.setNoData(self.no_data) - if self.geotransform: + if self.geotransform is not None: rda.geotransform = np.array(self.geotransform, dtype="float64") else: print( @@ -216,7 +220,7 @@ def copyFromWrapped(self, wrapped): class rd3array(np.ndarray): - def __new__(cls, array, meta_obj=None, no_data=None, order=None, **kwargs): + def __new__(cls, array, meta_obj=None, no_data=None, order=None, **kwargs: Any) -> "rd3array": obj = np.asarray(array, dtype=np.float32, order=order).view(cls) if meta_obj is not None: @@ -329,7 +333,7 @@ def LoadGDAL(filename: str, no_data: Optional[float] = None) -> rdarray: srcdata.metadata[k] = v _AddAnalysis( - srcdata, "LoadGDAL(filename={0}, no_data={1})".format(filename, no_data) + srcdata, f"LoadGDAL(filename={filename}, no_data={no_data})" ) return srcdata @@ -393,7 +397,7 @@ def FillDepressions(dem: rdarray, epsilon: bool = False, in_place: bool = False, if not in_place: dem = dem.copy() - _AddAnalysis(dem, "FillDepressions(dem, epsilon={0})".format(epsilon)) + _AddAnalysis(dem, f"FillDepressions(dem, epsilon={epsilon})") demw = dem.wrap() @@ -467,9 +471,7 @@ def ResolveFlats(dem: rdarray, in_place: bool = False) -> Optional[rdarray]: if not in_place: dem = dem.copy() - _AddAnalysis( - dem, "ResolveFlats(dem, in_place={in_place})".format(in_place=in_place) - ) + _AddAnalysis(dem, f"ResolveFlats(dem, in_place={in_place})") demw = dem.wrap() @@ -702,10 +704,7 @@ def FlowProportions(dem: rdarray, method: Optional[str] = None, exponent: Option _AddAnalysis( fprops, - "FlowProportions(dem, method={method}, exponent={exponent})".format( - method=method, - exponent=exponent, - ), + f"FlowProportions(dem, method={method}, exponent={exponent})", ) if method in fprop_methods: @@ -781,7 +780,8 @@ def TerrainAttribute(dem: rdarray, attrib: str, zscale: float = 1.0) -> rdarray: resultw = result.wrap() _AddAnalysis( - result, "TerrainAttribute(dem, attrib={0}, zscale={1})".format(attrib, zscale) + result, + f"TerrainAttribute(dem, attrib={attrib}, zscale={zscale})" ) terrain_attribs[attrib](dem.wrap(), resultw, zscale) @@ -800,5 +800,41 @@ def GeneratePerlinTerrain(size: int, seed: int) -> rdarray: Returns: A DEM """ + dem = rdarray(np.zeros(shape=(size, size), dtype="float64"), no_data=-9999) + dem.geotransform = STANDARD_GEOTRANSFORM + demw = dem.wrap() + _richdem.generate_perlin_terrain(demw, seed) + _AddAnalysis(dem, f"GeneratePerlinTerrain(size={size}, seed={seed})") + dem.copyFromWrapped(demw) + return dem + +def get_depression_hierarchy(dem: rdarray, labels: rdarray) -> Tuple[List[depression_hierarchy.Depression], rdarray]: + """Fills all depressions in a DEM. + + Args: + dem: An elevation model + labels: Should have `OCEAN` for cells representing the "ocean" (the + place to which depressions drain) and `NO_DEP` for all other + cells. + + Returns: + label: Modified in-place to contain a label indicating which depression + the cell belongs to. The indicated label is always the leaf of + the depression hierarchy, or the OCEAN. + + flowdirs: A value [0,7] indicated which direction water from the cell + flows in order to go "downhill". All cells have a flow + direction (even flats) except for pit cells. + """ + if type(dem) is not rdarray: + raise Exception("A richdem.rdarray or numpy.ndarray is required!") + + demw = dem.wrap() + labelsw = labels.wrap() + + flowdirs = rdarray(_richdem.NO_FLOW * np.ones(dem.shape, dtype=np.int8), no_data=-9999, geotransform=STANDARD_GEOTRANSFORM) + flowdirsw = flowdirs.wrap() + + dhret = depression_hierarchy.get_depression_hierarchy(demw, labelsw, flowdirsw) - return _richdem.generate_perlin_terrain(size, seed) + return dhret, flowdirs diff --git a/wrappers/pyrichdem/setup.py b/wrappers/pyrichdem/setup.py index 0d5c01b5..410bfd40 100644 --- a/wrappers/pyrichdem/setup.py +++ b/wrappers/pyrichdem/setup.py @@ -90,7 +90,7 @@ def build_extensions(self): ("DOCTEST_CONFIG_DISABLE", None), ("RICHDEM_COMPILE_TIME", f'"\\"{richdem_compile_time}\\""'), ("RICHDEM_GIT_HASH", f'"\\"{richdem_git_hash}\\""'), - ("RICHDEM_LOGGING", None), + # ("RICHDEM_LOGGING", None), ( "_USE_MATH_DEFINES", None, @@ -143,7 +143,7 @@ def build_extensions(self): # ':python_version < "3.4"': [ # 'numpy>=1.7,<1.12' # ] - # }, + # }, #python_requires = ' >= 2.6, !=3.0.*, !=3.1.*, !=3.2.*, <4', #TODO: https://pypi.python.org/pypi?%3Aaction=list_classifiers diff --git a/wrappers/pyrichdem/src/pywrapper.cpp b/wrappers/pyrichdem/src/pywrapper.cpp index be562334..d543ba03 100644 --- a/wrappers/pyrichdem/src/pywrapper.cpp +++ b/wrappers/pyrichdem/src/pywrapper.cpp @@ -1,5 +1,7 @@ #include "pywrapper.hpp" +#include + #include #include #include @@ -14,6 +16,8 @@ using namespace richdem; PYBIND11_MODULE(_richdem, m) { m.doc() = "Internal library used by pyRichDEM for calculations"; + m.attr("NO_FLOW") = &richdem::NO_FLOW; + //py::bind_vector>(m, "VecDouble"); py::bind_map>(m, "MapStringString"); @@ -129,5 +133,35 @@ PYBIND11_MODULE(_richdem, m) { } ); - m.def("generate_perlin_terrain", &perlin, "Generate random terrain using perlin noise", py::arg("size"), py::arg("seed")); + m.def("generate_perlin_terrain", &richdem::generate_perlin_terrain, "Generate random terrain using perlin noise", py::arg("array"), py::arg("seed")); + + py::module_ dephier_module = m.def_submodule("depression_hierarchy", "Depression Hierarchies"); + + dephier_module.attr("NO_PARENT") = &dephier::NO_PARENT; + dephier_module.attr("NO_VALUE") = &dephier::NO_VALUE; + dephier_module.attr("NO_DEP") = &dephier::NO_DEP; + dephier_module.attr("OCEAN") = &dephier::OCEAN; + + // Depression Hierarchy + py::class_>(dephier_module, "Depression") + .def(py::init<>()) + .def_readwrite("pit_cell", &dephier::Depression::pit_cell, "Flat index of the pit cell, the lowest cell in the depression. If more than one cell shares this lowest elevation, then one is arbitrarily chosen.") + .def_readwrite("out_cell", &dephier::Depression::out_cell, "Flat index of the outlet cell. If there is more than one outlet cell at this cell's elevation, then one is arbitrarily chosen.") + .def_readwrite("parent", &dephier::Depression::parent, "Parent depression. If both this depression and its neighbour fill up, this parent depression is the one which will contain the overflow.") + .def_readwrite("odep", &dephier::Depression::odep, "Outlet depression. The metadepression into which this one overflows. Usually its neighbour depression, but sometimes the ocean.") + .def_readwrite("geolink", &dephier::Depression::geolink, "When a metadepression overflows it does so into the metadepression indicated by `odep`. However, odep must flood from the bottom up. Therefore, we keep track of the `geolink`, which indicates what leaf depression the overflow is initially routed into.") + .def_readwrite("pit_elev", &dephier::Depression::pit_elev, "Elevation of the pit cell. Since the pit cell has the lowest elevation of any cell in the depression, we initialize this to infinity.") + .def_readwrite("out_elev", &dephier::Depression::out_elev, "Elevation of the outlet cell. Since the outlet cell has the lowest elevation of any path leading from a depression, we initialize this to infinity.") + .def_readwrite("lchild", &dephier::Depression::lchild, "The depressions form a binary tree. Each depression has two child depressions: one left and one right.") + .def_readwrite("rchild", &dephier::Depression::rchild, "The depressions form a binary tree. Each depression has two child depressions: one left and one right.") + .def_readwrite("ocean_parent", &dephier::Depression::ocean_parent, "Indicates whether the parent link is to either the ocean or a depression that links to the ocean.") + .def_readwrite("ocean_linked", &dephier::Depression::ocean_linked, "Indicates depressions which link to the ocean through this depression, but are not subdepressions. That is, these ocean-linked depressions may be at the top of high cliffs and spilling into this depression.") + .def_readwrite("dep_label", &dephier::Depression::dep_label, "The label of the depression, for calling it up again.") + .def_readwrite("cell_count", &dephier::Depression::cell_count, "Number of cells contained within the depression and its children.") + .def_readwrite("dep_vol", &dephier::Depression::dep_vol, "Volume of the depression and its children. Used in the Water Level Equation (see below).") + .def_readwrite("water_vol", &dephier::Depression::water_vol, "Water currently contained within the depression. Used in the Water Level Equation (see below).") + .def_readwrite("total_elevation", &dephier::Depression::total_elevation, "Total elevation of cells contained with the depression and its children.") + ; //Ends the class definition above + + dephier_module.def("get_depression_hierarchy", &dephier::GetDepressionHierarchy, "Calculate the hierarchy of depressions. Takes as input a digital elevation model and a set of labels. The labels should have `OCEAN` for cells"); } diff --git a/wrappers/pyrichdem/src/pywrapper.hpp b/wrappers/pyrichdem/src/pywrapper.hpp index 8bb6fdef..c6419ef5 100644 --- a/wrappers/pyrichdem/src/pywrapper.hpp +++ b/wrappers/pyrichdem/src/pywrapper.hpp @@ -1,3 +1,5 @@ +#pragma once + #include #include #include diff --git a/wrappers/pyrichdem/tests/tests.py b/wrappers/pyrichdem/tests/tests.py new file mode 100644 index 00000000..85b50f6a --- /dev/null +++ b/wrappers/pyrichdem/tests/tests.py @@ -0,0 +1,34 @@ +import unittest + +import numpy as np +import richdem as rd +from richdem import depression_hierarchy as dephier + +class py_richdem_tests(unittest.TestCase): + def test_depression_defaults(self) -> None: + d = dephier.Depression() + self.assertEqual(d.out_cell, dephier.NO_VALUE) + self.assertEqual(d.pit_cell, dephier.NO_VALUE) + self.assertEqual(d.parent, dephier.NO_PARENT) + self.assertEqual(d.odep, dephier.NO_VALUE) + self.assertEqual(d.geolink, dephier.NO_VALUE) + self.assertEqual(d.lchild, dephier.NO_VALUE) + self.assertEqual(d.rchild, dephier.NO_VALUE) + self.assertEqual(d.ocean_parent, False) + self.assertEqual(d.dep_label, 0) + self.assertEqual(d.cell_count, 0) + self.assertEqual(d.dep_vol, 0) + self.assertEqual(d.water_vol, 0) + self.assertEqual(d.total_elevation, 0) + + def test_random_terrain(self) -> None: + dem = rd.GeneratePerlinTerrain(20, 20) + # A labels array where all the edge cells are in the ocean and all the + # interior cells are not yet assigned to a depression + labels = rd.rdarray( + dephier.OCEAN*np.ones(dem.shape, dtype = np.uint32), + no_data=-9999, + geotransform=rd.STANDARD_GEOTRANSFORM + ) + labels[1:-1,1:-1] = dephier.NO_DEP + dh, flowdirs = rd.get_depression_hierarchy(dem, labels) \ No newline at end of file