diff --git a/README.md b/README.md index 2e01fd0..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. @@ -152,13 +153,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/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 0aa1d39..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, @@ -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/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/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 9c85336..8565565 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 @@ -145,15 +148,21 @@ 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) + out._fobj = None + 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 -def sameExtend(fobjs): +def sameExtent(fobjs): bbox = commonBbox(fobjs) for fobj in fobjs: if fobj.bbox != bbox: @@ -161,20 +170,32 @@ def sameExtend(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 - 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("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)) +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)) + 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): - 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) @@ -186,6 +207,13 @@ def maskData(data, mask): def main(config, gauges): + flowacc = None + flowdir = None + filedict_main = {} + gaugedict_main = {} + gaugefile_main = None + updated_gauge = None + for gauge in gauges: logging.info("processing gauge: %s", gauge.id) @@ -194,35 +222,44 @@ 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"]) - 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") - gauge = matchFlowacc(gauge, flowacc, **config["matching"]) + updated_gauge, error = matchFlowacc(gauge, flowacc, **config["matching"]) - if not 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") - 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") fitem = GridFile( - fname=config["gauge"].get("fname", "idgauges.asc"), + fname=fname, outpath=config["gauge"].get("outpath"), ) - gaugefile = gaugeGrid(flowacc, 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("reding gauge file") + logging.debug("reading gauge file") mask = gridBasinMask(gauge) for fdict in config.get("gridfiles", []): @@ -253,19 +290,23 @@ 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) writeFiles(bpath, filedict) logging.debug("writing report") - writeReport(bpath, mask, config["matching"]["scaling_factor"], gauge) + 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():