add multi-gauge facilitation

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,15 +37,15 @@ 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)
 ```
 
 #### MacOS
 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

setup.cfg

@@ -51,6 +51,7 @@ install_requires =
     pyyaml
     gdal
     netcdf4<1.6
+    pandas
 python_requires = >=3.6
 zip_safe = False
 

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

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))
 

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
     ------------------

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,36 +148,54 @@ 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:
             return False
     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():