From 0fabe7c99694910aae9801e8ed3dd0990f000413 Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Mon, 26 Sep 2022 17:26:07 +0200 Subject: [PATCH 1/7] - added optional read of flwdir and flwacc fields --- src/basinex/main.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/basinex/main.py b/src/basinex/main.py index 9c85336..a9191b4 100644 --- a/src/basinex/main.py +++ b/src/basinex/main.py @@ -186,6 +186,9 @@ def maskData(data, mask): def main(config, gauges): + flowacc = None + flowdir = None + for gauge in gauges: logging.info("processing gauge: %s", gauge.id) @@ -194,11 +197,13 @@ def main(config, gauges): if not gauge.path: # create mask if not given - logging.debug("reading flow accumulation") - flowacc = ga.fromfile(config["flowacc"]).astype(np.int32) + if flowacc is None: + logging.debug("reading flow accumulation") + flowacc = ga.fromfile(config["flowacc"]).astype(np.int32) - logging.debug("reading flow direction") - flowdir = ga.fromfile(config["flowdir"]).astype(np.int32) + if flowdir is None: + logging.debug("reading flow direction") + flowdir = ga.fromfile(config["flowdir"]).astype(np.int32) if gauge.size: logging.debug("moving gauge to streamflow") From 0d62d88cd70940a5cf54e19e548e552fcc8c197d Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Tue, 27 Sep 2022 01:27:20 +0200 Subject: [PATCH 2/7] - fixed typos --- README.md | 4 ++-- src/basinex/geoarray/wrapper.py | 2 +- src/basinex/main.py | 10 +++++----- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 2e01fd0..fcb821c 100644 --- a/README.md +++ b/README.md @@ -152,13 +152,13 @@ ncfiles: - `path`: path to the mask file - `varname`: name of the mask variable (optional, only needed if the mask is stored in a netcdf file) - `latitude-size-correction: False` - **Optional**: - perform a latitude correction for the given basin size (default: False) + perform a latitude correction for the basin size of a given gauge (default: False) - `AREA = N_cells * res_x * ( cos(LAT) * res_y ) * scaling factor^2` - `matching:` - **Required**: gauge matching parameters - **Note**: The gauge matching is based on the flowaccumulation data. The value for any given cell in the flowaccumulation grid is interpreted as the size - [in cells] of a river basin drainig into the respective cell. + [in cells] of a river basin draining into the respective cell. During gauge matching the flowaccumulation grid is searched for a cell with a corresponding basin size close to the given gauge basin size. The search radius will be increased succesively and can be limited to a diff --git a/src/basinex/geoarray/wrapper.py b/src/basinex/geoarray/wrapper.py index 9dc0524..e39e702 100644 --- a/src/basinex/geoarray/wrapper.py +++ b/src/basinex/geoarray/wrapper.py @@ -190,7 +190,7 @@ def full(shape, value, dtype=np.float64, *args, **kwargs): Arguments --------- shape : tuple # shape of the returned grid - fill_value : scalar # fille value + fill_value : scalar # fill value Optional Arguments ------------------ diff --git a/src/basinex/main.py b/src/basinex/main.py index a9191b4..5053f22 100644 --- a/src/basinex/main.py +++ b/src/basinex/main.py @@ -153,7 +153,7 @@ def gaugeGrid(grid_template, gauge): return out -def sameExtend(fobjs): +def sameExtent(fobjs): bbox = commonBbox(fobjs) for fobj in fobjs: if fobj.bbox != bbox: @@ -227,7 +227,7 @@ def main(config, gauges): filedict[fitem] = maskData(gaugefile, mask) else: - logging.debug("reding gauge file") + logging.debug("reading gauge file") mask = gridBasinMask(gauge) for fdict in config.get("gridfiles", []): @@ -258,13 +258,13 @@ def main(config, gauges): filedict[fitem] = mask if filedict: - logging.debug("finding common extend") + logging.debug("finding common extent") bbox = commonBbox(tuple(filedict.values())) - logging.debug("enlarging data to common extend") + logging.debug("enlarging data to common extent") filedict = enlargeFiles(filedict, bbox) - if not sameExtend(tuple(filedict.values())): + if not sameExtent(tuple(filedict.values())): raise RuntimeError("incompatible cellsizes") bpath = os.path.join(config["outpath"], gauge.id) From 633c073e2a65254b5d25c70d5b4a61f60cc15cfe Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Tue, 27 Sep 2022 01:29:09 +0200 Subject: [PATCH 3/7] - presumably fixed bug in !2, bbox comparison was imprecise due to rounding errors --- src/basinex/geoarray/spatial.py | 19 +++++++++++-------- src/basinex/main.py | 2 +- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/src/basinex/geoarray/spatial.py b/src/basinex/geoarray/spatial.py index 8f4b086..0a32644 100644 --- a/src/basinex/geoarray/spatial.py +++ b/src/basinex/geoarray/spatial.py @@ -4,6 +4,9 @@ import numpy as np +# precision for rounding to avoid numerical instabilities +PRECISION = 10 + class SpatialMixin(object): def trim(self): @@ -83,10 +86,10 @@ def shrink(self, ymin=None, ymax=None, xmin=None, xmax=None): } cellsize = [float(abs(cs)) for cs in self.cellsize] - top = floor((self.bbox["ymax"] - bbox["ymax"]) / cellsize[0]) - left = floor((bbox["xmin"] - self.bbox["xmin"]) / cellsize[1]) - bottom = floor((bbox["ymin"] - self.bbox["ymin"]) / cellsize[0]) - right = floor((self.bbox["xmax"] - bbox["xmax"]) / cellsize[1]) + top = floor(round((self.bbox["ymax"] - bbox["ymax"]) / cellsize[0], PRECISION)) + left = floor(round((bbox["xmin"] - self.bbox["xmin"]) / cellsize[1], PRECISION)) + bottom = floor(round((bbox["ymin"] - self.bbox["ymin"]) / cellsize[0], PRECISION)) + right = floor(round((self.bbox["xmax"] - bbox["xmax"]) / cellsize[1], PRECISION)) return self.removeCells( max(top, 0), max(left, 0), max(bottom, 0), max(right, 0) @@ -175,10 +178,10 @@ def enlarge(self, ymin=None, ymax=None, xmin=None, xmax=None): cellsize = [float(abs(cs)) for cs in self.cellsize] - top = ceil((bbox["ymax"] - self.bbox["ymax"]) / cellsize[0]) - left = ceil((self.bbox["xmin"] - bbox["xmin"]) / cellsize[1]) - bottom = ceil((self.bbox["ymin"] - bbox["ymin"]) / cellsize[0]) - right = ceil((bbox["xmax"] - self.bbox["xmax"]) / cellsize[1]) + top = ceil(round((bbox["ymax"] - self.bbox["ymax"]) / cellsize[0], PRECISION)) + left = ceil(round((self.bbox["xmin"] - bbox["xmin"]) / cellsize[1], PRECISION)) + bottom = ceil(round((self.bbox["ymin"] - bbox["ymin"]) / cellsize[0], PRECISION)) + right = ceil(round((bbox["xmax"] - self.bbox["xmax"]) / cellsize[1], PRECISION)) return self.addCells(max(top, 0), max(left, 0), max(bottom, 0), max(right, 0)) diff --git a/src/basinex/main.py b/src/basinex/main.py index 5053f22..af4c6f9 100644 --- a/src/basinex/main.py +++ b/src/basinex/main.py @@ -174,7 +174,7 @@ def writeReport(bpath, mask, scaling_factor, gauge): def maskData(data, mask): - if all(x == y for x, y in zip(mask.cellsize, data.cellsize)): + if all(np.isclose(x, y) for x, y in zip(mask.cellsize, data.cellsize)): return data.setMask(mask.mask) enlarged_mask = mask.enlarge(**data.bbox).astype(float) From 915e3c3598da745d4e715a97aaef1047479b55eb Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Tue, 27 Sep 2022 01:33:25 +0200 Subject: [PATCH 4/7] - refactored writeReport to be more precise - implemented first version of batchprocessing of idgauges.asc file (still buggy, facc values are written although full_like is used) --- src/basinex/gauges.py | 4 +++- src/basinex/main.py | 52 ++++++++++++++++++++++++++++++------------- 2 files changed, 39 insertions(+), 17 deletions(-) diff --git a/src/basinex/gauges.py b/src/basinex/gauges.py index 0aa1d39..50f220f 100755 --- a/src/basinex/gauges.py +++ b/src/basinex/gauges.py @@ -99,4 +99,6 @@ def matchFlowacc(gauge, facc, max_distance, max_error, scaling_factor=1): # the cell cordinates y, x = grid.coordinatesOf(river_cells_y[nn], river_cells_x[nn]) - return Gauge(id=gauge.id, y=y, x=x, size=size) + return Gauge(id=gauge.id, y=y, x=x, size=grid.data[river_cells_y[nn], river_cells_x[nn]]), error + else: + return None, error diff --git a/src/basinex/main.py b/src/basinex/main.py index af4c6f9..822e688 100644 --- a/src/basinex/main.py +++ b/src/basinex/main.py @@ -145,9 +145,14 @@ def commonBbox(fobjs): return bbox -def gaugeGrid(grid_template, gauge): - out = ga.full_like(grid_template, grid_template.fill_value) +def gaugeGrid(grid_template, gauge, out_in=None): + if out_in is None: + out = ga.full_like(grid_template, grid_template.fill_value) + else: + out = out_in idx = out.indexOf(gauge.y, gauge.x) + if out_in is not None and not out.mask[idx]: + warnings.warn(f"There is already a gauge at {idx}, with id '{out.data[idx]}', replacing it by '{gauge.id}'.") out.data[idx] = gauge.id out.mask[idx] = False return out @@ -161,16 +166,17 @@ def sameExtent(fobjs): return True -def writeReport(bpath, mask, scaling_factor, gauge): - size = (np.sum(~mask.mask) * np.prod(np.abs(mask.cellsize))) * scaling_factor**2 - error_size = (size - gauge.size) / gauge.size * 100 +def writeReport(bpath, updated_gauge, gauge, error): Path(bpath).mkdir(exist_ok=True, parents=True) with open(os.path.join(bpath, "report.out"), "w") as f: - f.write("calculated_catchment_size: {:}\n".format(size)) + f.write("error_catchment_size (%) : {:}\n".format(error)) + if updated_gauge is not None: + f.write("new_catchment_size : {:}\n".format(updated_gauge.size)) + f.write("new_y : {:}\n".format(updated_gauge.y)) + f.write("new_x : {:}\n".format(updated_gauge.x)) f.write("input_catchment_size : {:}\n".format(gauge.size)) - f.write("error_catchment_size (%) : {:}\n".format(error_size)) - f.write("adjusted_y : {:}\n".format(gauge.y)) - f.write("adjusted_x : {:}\n".format(gauge.x)) + f.write("input_y : {:}\n".format(gauge.y)) + f.write("input_x : {:}\n".format(gauge.x)) def maskData(data, mask): @@ -188,6 +194,9 @@ def main(config, gauges): flowacc = None flowdir = None + filedict_main = {} + gaugefile_main = None + updated_gauge = None for gauge in gauges: logging.info("processing gauge: %s", gauge.id) @@ -199,7 +208,7 @@ def main(config, gauges): # create mask if not given if flowacc is None: logging.debug("reading flow accumulation") - flowacc = ga.fromfile(config["flowacc"]).astype(np.int32) + flowacc = ga.fromfile(config["flowacc"]) if flowdir is None: logging.debug("reading flow direction") @@ -207,23 +216,32 @@ def main(config, gauges): if gauge.size: logging.debug("moving gauge to streamflow") - gauge = matchFlowacc(gauge, flowacc, **config["matching"]) + updated_gauge, error = matchFlowacc(gauge, flowacc, **config["matching"]) - if not gauge: + if not updated_gauge: warnings.warn("Failed to match the gauge to the flow accumulation grid") continue logging.debug("generating basin mask") - mask = gaugeBasinMask(flowdir, gauge) + mask = gaugeBasinMask(flowdir, updated_gauge) # write gauge grid if desired if "gauge" in config: logging.debug("writing gauge file") + fname = config["gauge"].get("fname", "idgauges.asc") + if len(gauges) > 1: + if gaugefile_main is None: + fitem_main = GridFile( + fname=fname, + outpath=config["gauge"].get("outpath"), + ) + gaugefile_main = gaugeGrid(flowacc, updated_gauge, gaugefile_main) + filedict_main[fitem_main] = gaugefile_main fitem = GridFile( - fname=config["gauge"].get("fname", "idgauges.asc"), + fname=fname, outpath=config["gauge"].get("outpath"), ) - gaugefile = gaugeGrid(flowacc, gauge).shrink(**mask.bbox) + gaugefile = gaugeGrid(flowacc, updated_gauge).shrink(**mask.bbox) filedict[fitem] = maskData(gaugefile, mask) else: @@ -270,7 +288,9 @@ def main(config, gauges): bpath = os.path.join(config["outpath"], gauge.id) writeFiles(bpath, filedict) logging.debug("writing report") - writeReport(bpath, mask, config["matching"]["scaling_factor"], gauge) + writeReport(bpath, updated_gauge, gauge, error) + + writeFiles(config["outpath"], filedict_main) def initArgparser(): From 15268675f5dd31b906c9aa650957b111c151f1c7 Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Wed, 28 Sep 2022 17:31:47 +0200 Subject: [PATCH 5/7] - fixed small bug preventing correct writing of idgauges.asc --- src/basinex/main.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/basinex/main.py b/src/basinex/main.py index 822e688..a8e2bdb 100644 --- a/src/basinex/main.py +++ b/src/basinex/main.py @@ -148,6 +148,7 @@ def commonBbox(fobjs): def gaugeGrid(grid_template, gauge, out_in=None): if out_in is None: out = ga.full_like(grid_template, grid_template.fill_value) + out._fobj = None else: out = out_in idx = out.indexOf(gauge.y, gauge.x) From 7cee75f5a494087174ace772c541c3493f3fb539 Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Fri, 11 Nov 2022 00:01:04 +0100 Subject: [PATCH 6/7] added complete summary as table of moving gauge to streamflow, added pandas dependency, fixed bug when writing basin mask, refactored types of gauge attributes --- README.md | 7 +++-- setup.cfg | 1 + src/basinex/gauges.py | 12 ++++---- src/basinex/main.py | 69 ++++++++++++++++++++++++++----------------- 4 files changed, 53 insertions(+), 36 deletions(-) diff --git a/README.md b/README.md index fcb821c..a3494c2 100644 --- a/README.md +++ b/README.md @@ -16,6 +16,7 @@ The [mHM](https://mhm-ufz.org/) basin extractor. Extract basins for given gaugin ## Dependencies - numpy v1.14.5 or later +- pandas - netCDF4 - GDAL - pyyaml @@ -36,7 +37,7 @@ To get a recent version of GDAL, you can use the ppa of [ubuntugis](https://laun sudo add-apt-repository ppa:ubuntugis/ppa sudo apt-get update sudo apt install gdal-bin libgdal-dev -pip install wheel numpy +pip install wheel numpy pandas pip install GDAL==$(gdal-config --version) ``` @@ -44,7 +45,7 @@ pip install GDAL==$(gdal-config --version) GDAL can be installed with [homebrew](https://formulae.brew.sh/formula/gdal): ``` brew install gdal -pip install wheel numpy +pip install wheel numpy pandas pip install GDAL==$(gdal-config --version) ``` @@ -60,7 +61,7 @@ pipwin install gdal It is best to use basinex with conda to have gdal and NetCDF installed properly. To use the development version of basinex, download this repository and do the following in your conda environment: - conda install -y gdal netcdf4 pyyaml cxx-compiler + conda install -y gdal netcdf4 pyyaml cxx-compiler pandas pip install . Then you can execute `basinex` in that conda environment. diff --git a/setup.cfg b/setup.cfg index 7f31553..8440d3a 100644 --- a/setup.cfg +++ b/setup.cfg @@ -51,6 +51,7 @@ install_requires = pyyaml gdal netcdf4<1.6 + pandas python_requires = >=3.6 zip_safe = False diff --git a/src/basinex/gauges.py b/src/basinex/gauges.py index 50f220f..a915ca7 100755 --- a/src/basinex/gauges.py +++ b/src/basinex/gauges.py @@ -10,9 +10,9 @@ def __init__( self, id, y=None, x=None, size=None, path=None, varname=None, lat_fix=False ): self.id = id - self.y = y - self.x = x - self.size = float(size) / np.cos(np.deg2rad(float(y))) if lat_fix else size + self.y = float(y) + self.x = float(x) + self.size = float(size) / np.cos(np.deg2rad(self.y)) if lat_fix else float(size) self.path = path self.varname = varname @@ -66,9 +66,9 @@ def matchFlowacc(gauge, facc, max_distance, max_error, scaling_factor=1): to convert from map units to km^2 of the catchment area """ # print "gauge before:", gauge.y, gauge.x - y = float(gauge.y) - x = float(gauge.x) - size = float(gauge.size) + y = gauge.y + x = gauge.x + size = gauge.size bbox = { "ymin": y - max_distance, "ymax": y + max_distance, diff --git a/src/basinex/main.py b/src/basinex/main.py index a8e2bdb..2155851 100644 --- a/src/basinex/main.py +++ b/src/basinex/main.py @@ -7,6 +7,7 @@ from pathlib import Path import numpy as np +import pandas as pd import yaml from . import __version__ @@ -62,10 +63,10 @@ def gaugeBasinMask(flowdir, gauge): extract(np.array(flowdir, dtype=np.int32, copy=True), *gauge_idx), dtype=np.int32, ) - mask[mask == 0] = flowdir.fill_value - out = ga.array(mask, **flowdir.header) - return out.trim() + out = ga.array(mask, **flowdir.header).trim() + out._fobj = None + return out def gridBasinMask(gauge): @@ -89,6 +90,7 @@ def gridBasinMask(gauge): fill_value=var.fill_value, cellsize=nc.cellsize, ) + out._fobj = None return out.setMask(out <= 0) @@ -114,6 +116,7 @@ def openGridFiles(flist): flist = flist or () for fdict in flist: out[GridFile(**fdict)] = ga.fromfile(fdict["fname"]) + out._fobj = None return out @@ -167,17 +170,28 @@ def sameExtent(fobjs): return True -def writeReport(bpath, updated_gauge, gauge, error): - Path(bpath).mkdir(exist_ok=True, parents=True) - with open(os.path.join(bpath, "report.out"), "w") as f: - f.write("error_catchment_size (%) : {:}\n".format(error)) - if updated_gauge is not None: - f.write("new_catchment_size : {:}\n".format(updated_gauge.size)) - f.write("new_y : {:}\n".format(updated_gauge.y)) - f.write("new_x : {:}\n".format(updated_gauge.x)) - f.write("input_catchment_size : {:}\n".format(gauge.size)) - f.write("input_y : {:}\n".format(gauge.y)) - f.write("input_x : {:}\n".format(gauge.x)) +def writeReport(bpath, updated_gauge, gauge, error, do_write=False): + if do_write: + Path(bpath).mkdir(exist_ok=True, parents=True) + with open(os.path.join(bpath, "report.out"), "w") as f: + f.write("error_catchment_size (%) : {:}\n".format(error)) + if updated_gauge is not None: + f.write("new_catchment_size : {:}\n".format(updated_gauge.size)) + f.write("new_y : {:}\n".format(updated_gauge.y)) + f.write("new_x : {:}\n".format(updated_gauge.x)) + f.write("input_catchment_size : {:}\n".format(gauge.size)) + f.write("input_y : {:}\n".format(gauge.y)) + f.write("input_x : {:}\n".format(gauge.x)) + return {gauge.id: { + 'gauge_size': gauge.size, + 'new_size': updated_gauge.size, + 'error_size': error, + 'gauge_x': gauge.x, + 'new_x': updated_gauge.x, + 'gauge_y': gauge.y, + 'new_y': updated_gauge.y, + 'error_dist': ((gauge.x -updated_gauge.x)**2 + (gauge.y -updated_gauge.y)**2)**(0.5), + }} def maskData(data, mask): @@ -196,6 +210,7 @@ def main(config, gauges): flowacc = None flowdir = None filedict_main = {} + gaugedict_main = {} gaugefile_main = None updated_gauge = None @@ -219,8 +234,8 @@ def main(config, gauges): logging.debug("moving gauge to streamflow") updated_gauge, error = matchFlowacc(gauge, flowacc, **config["matching"]) - if not updated_gauge: - warnings.warn("Failed to match the gauge to the flow accumulation grid") + if updated_gauge is None: + warnings.warn(f"Failed to match the gauge {gauge.id} to the flow accumulation grid") continue logging.debug("generating basin mask") @@ -230,20 +245,18 @@ def main(config, gauges): if "gauge" in config: logging.debug("writing gauge file") fname = config["gauge"].get("fname", "idgauges.asc") - if len(gauges) > 1: - if gaugefile_main is None: - fitem_main = GridFile( - fname=fname, - outpath=config["gauge"].get("outpath"), - ) - gaugefile_main = gaugeGrid(flowacc, updated_gauge, gaugefile_main) - filedict_main[fitem_main] = gaugefile_main fitem = GridFile( fname=fname, outpath=config["gauge"].get("outpath"), ) - gaugefile = gaugeGrid(flowacc, updated_gauge).shrink(**mask.bbox) - filedict[fitem] = maskData(gaugefile, mask) + if len(gauges) > 1: + if gaugefile_main is None: + fitem_main = fitem + gaugefile_main = gaugeGrid(flowacc, updated_gauge, gaugefile_main) + filedict_main[fitem_main] = gaugefile_main + else: + gaugefile = gaugeGrid(flowacc, updated_gauge).shrink(**mask.bbox) + filedict[fitem] = maskData(gaugefile, mask) else: logging.debug("reading gauge file") @@ -289,9 +302,11 @@ def main(config, gauges): bpath = os.path.join(config["outpath"], gauge.id) writeFiles(bpath, filedict) logging.debug("writing report") - writeReport(bpath, updated_gauge, gauge, error) + gaugedict_main.update(writeReport(bpath, updated_gauge, gauge, error)) writeFiles(config["outpath"], filedict_main) + if gaugedict_main: + pd.DataFrame.from_dict(gaugedict_main, orient='index').to_csv(f'{config["outpath"]}/summary.csv') def initArgparser(): From 6573834aab772527d931afb96ef202c2bf44eb88 Mon Sep 17 00:00:00 2001 From: Robert Schweppe Date: Fri, 11 Nov 2022 00:09:46 +0100 Subject: [PATCH 7/7] made individual report dependent on verbose flag --- src/basinex/main.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/basinex/main.py b/src/basinex/main.py index 2155851..8565565 100644 --- a/src/basinex/main.py +++ b/src/basinex/main.py @@ -170,8 +170,8 @@ def sameExtent(fobjs): return True -def writeReport(bpath, updated_gauge, gauge, error, do_write=False): - if do_write: +def writeReport(bpath, updated_gauge, gauge, error): + if logging.DEBUG >= logging.root.level: Path(bpath).mkdir(exist_ok=True, parents=True) with open(os.path.join(bpath, "report.out"), "w") as f: f.write("error_catchment_size (%) : {:}\n".format(error))