Merging Development

.gitignore

@@ -1,6 +1,8 @@
 # python
 __pycache__
 *.py[cod]
+*$py.class
+
 # binary blobs
 *.nc
 *.nc4
@@ -13,9 +15,6 @@ __pycache__
 # mypy
 .mypy_cache/
 
-# vscode
-.vscode
-
 # os
 .directory
 .DS_Store
@@ -35,6 +34,14 @@ project/
 .ipynb_checkpoints/
 
 #build files 
+eggs/
+.eggs/
+*.egg
 *.egg-info/
+wheels/
 build/
 dist/
+
+#ide
+.idea/
+.vscode

readme.md

@@ -3,7 +3,7 @@ Python implementation of Tank Hydrologic Model, a conceptual rainfall-runoff mod
 
 
 <p align="center">
-<img align="center" height="600px"  src="https://raw.githubusercontent.com/nzahasan/tank-model/master/assets/tank-model-schamatic.svg" >
+<img align="center" height="500px"  src="https://raw.githubusercontent.com/nzahasan/tank-model/master/assets/tank-model-schamatic.svg" >
 </p>
 
 ### Installation
@@ -14,14 +14,41 @@ Tank-Model can be installed as a python package using the following commands
 $ pip install https://github.com/nzahasan/tank-model/zipball/master
 ```
 
+after successful installation `tank_cmd.py` should be available which can be used for setting up new project, optimizing the project and computation.
+
+```bash
+# get help text command line utility
+$ tank_cmd.py --help
+
+# get help text of subcommand
+$ tank_cmd.py new-project --help
+```
 
 ### Setting up a new model:
 
-New project can be created using the following command. This creates a folder with a project defination inside it.
+New project can be created using the following command. This command creates a folder in working directory with a json formatted project definition inside it.
 ```bash
 $ tank_cmd.py new_project project_name
 ```
 
+A sample project definition looks like this
+```json
+{
+  "interval": 24.0,
+  "basin": "sample_project.basin.json",
+  "precipitation": "sample_project.pr.csv",
+  "evapotranspiration": "sample_project.et.csv",
+  "discharge": "sample_project.q.csv",
+  "result": "sample_project.result.csv",
+  "statistics": "sample_project.stats.json"
+}
+```
+here `interval` is the time step of simulation in hours. The other attributes are file locations; `precipitation`, `evapotranspiration` and `discharge` are CSV files containing time-series data. These files should be formatted according to the file format mentioned here <a href="file-format-spec.md">file-format-spec.md</a>
+
+`precipitation` & `evapotranspiration` serve as input data for the model simulation and resulting output is stored in the `result` file following the time-series CSV format mentioned earlier. Data in the `discharge` is used for model calibration. And performance matrices are stored in the `statistics` file.
+
+
+
 ### Converting HEC-HMS basin to tank model basin
 
 ```bash

scripts/gefs_download.py

@@ -1,12 +1,12 @@
 #!/usr/bin/env python3
 
-import requests as req 
-import click
-from multiprocessing.pool import ThreadPool
-from pathlib import Path
 import time
+import click
 import shutil
 import subprocess
+import requests as req 
+from pathlib import Path
+from multiprocessing.pool import ThreadPool
 
 SLEEP_TIME = 2*60 # in seconds
 
@@ -50,7 +50,7 @@ def download_file(parameters:tuple, logging=True) -> None:
         time.sleep(SLEEP_TIME)
         download_file(parameters)
 
-def merge_grib_convert_nc(output_dir):
+def merge_grib_convert_nc(output_dir, cycle):
 
     for ens in range(31):
 
@@ -59,7 +59,7 @@ def merge_grib_convert_nc(output_dir):
         file_paths = list()
 
         for hr in range(3,243,3):
-            _fname = f'{_pre}{ens:02d}.t00z.pgrb2s.0p25.f{hr:03d}'
+            _fname = f'{_pre}{ens:02d}.t{cycle}z.pgrb2s.0p25.f{hr:03d}'
 
             _full_path = output_dir / _fname 
 
@@ -68,7 +68,6 @@ def merge_grib_convert_nc(output_dir):
         file_checks = [fpath.exists() for fpath in file_paths]
 
 
-
         if False in file_checks:
             print(f'file check# all files are not available')
             print(f'file check# skipping for ens - {ens:02}')
@@ -124,7 +123,7 @@ def main(date:str, cycle:str, left:float, right:float, bottom:float, top:float,
 
         for hr in range(3,243,3):
 
-            _fname = f'{_pre}{ens:02d}.t00z.pgrb2s.0p25.f{hr:03d}'
+            _fname = f'{_pre}{ens:02d}.t{cycle}z.pgrb2s.0p25.f{hr:03d}'
 
             full_url  = f"https://nomads.ncep.noaa.gov/cgi-bin/filter_gefs_atmos_0p25s.pl?"
             full_url += f"dir=/gefs.{date}/{cycle}/atmos/pgrb2sp25"
@@ -141,7 +140,7 @@ def main(date:str, cycle:str, left:float, right:float, bottom:float, top:float,
     pool.close()
 
     # merge grib files and convert to netcdf4
-    merge_grib_convert_nc(output_dir)
+    merge_grib_convert_nc(output_dir, cycle)
 
 
 

scripts/tank_cmd.py

@@ -184,7 +184,8 @@ def optimize(project_file):
     '''Automatically optimizes tank basin parameters for a given projects'''
     project_dir = os.path.dirname(os.path.abspath(project_file))
 
-    project = ioh.read_project_file(project_file)
+    # project must have discharge as this is mandatory 
+    project = ioh.read_project_file(project_file, check_discharge_file=True)
 
     basin_file = os.path.join(project_dir, project['basin'])
     precipitation_file = os.path.join(project_dir, project['precipitation'])

tank_core/__init__.py

@@ -6,5 +6,5 @@
 Tank-Model is a conceptual rainfall-runoff model 
 proposed by Sugawara and Funiyuki (1956).
 '''
-__version__ = '0.2.5'
+__version__ = '0.2.6'
 __author__ = 'Nazmul Ahasan <nzahasan@gmail.com>'

tank_core/arima.py

@@ -3,7 +3,7 @@
 Time series error Correction Module for model
 ---------------------------------------------
 ARIMA(p,d,q)
-    p - autoregressive part
+    p - auto-regressive part
     d - integrated part / differentiation part
     q - moving average part
 

tank_core/channel_routing.py

@@ -19,15 +19,15 @@ def muskingum(in_flow:np.ndarray, del_t:float, k:float, x:float) -> np.ndarray:
     x        - hr
     '''
 
-    # calculate timestep
+    # calculate time-step
     n_step:int = in_flow.shape[0]
 
     # create a zero array of out_flow
     out_flow:np.ndarray = np.zeros(n_step, dtype=np.float64)
 
     C0:float = (-k*x+0.5*del_t) / (k*(1-x)+0.5*del_t)
-    C1:float = ( k*x+0.5*del_t) / (k*(1-x)+0.5*del_t)
-    C2:float = (k*(1-x)-0.5*del_t) / (k*(1-x)+0.5*del_t)
+    C1:float = (k*x+0.5 *del_t) / (k*(1-x)+0.5*del_t)
+    C2:float = (k*(1-x) - 0.5*del_t) / (k*(1-x)+0.5*del_t)
 
     # constraints check
     if (C0+C1+C2) > 1 or x >0.5 or (del_t/k + x) > 1:

tank_core/computation_helpers.py

@@ -6,17 +6,10 @@
 import pandas as pd
 from queue import Queue
 from scipy.optimize import minimize
+from . import utils
 from .tank_basin import tank_discharge
 from .channel_routing import muskingum
-from . import utils
-from .cost_functions import (
-    PBIAS, 
-    R2, 
-    RMSE, 
-    NSE, 
-    MSE, 
-    PBIAS
-)
+from .cost_functions import (PBIAS, R2, RMSE, NSE,)
 from . import global_config as gc
 
 
@@ -85,6 +78,7 @@ def compute_project(
     n_step = len(precipitation.index)
 
     computation_result = pd.DataFrame(index=precipitation.index)
+
     model_states = dict(
         time = precipitation.index.to_numpy()
     )

tank_core/cost_functions.py

@@ -26,9 +26,14 @@
         Lower is better 
         +ve values indicate underestimation
         -ve values indicate model overestimation
+    KGE:
+        Kling-Gupta efficiency
+        KGE = 1 indicates perfect agreement
+        KGE -ve is bad model
 '''
 
 import numpy as np
+from scipy.stats import pearsonr
 from .utils import shape_alike
 
 def R2(x:np.ndarray, y:np.ndarray)->float:
@@ -85,4 +90,14 @@ def PBIAS(obs:np.ndarray, sim:np.ndarray)->float:
     if not shape_alike(sim,obs):
         raise Exception('shape mismatch between x and y')
 
-    return (obs-sim).sum() * 100 / obs.sum()
+    return (obs-sim).sum() * 100 / obs.sum()
+
+def KGE(obs:np.ndarray, sim:np.ndarray)->float:
+    """
+    Kling-Gupta efficiency
+    """
+    eMean = (np.mean(sim) / np.mean(obs)) - 1 
+    eVar = (np.std(sim) / np.std(obs)) - 1 
+    eCor = pearsonr(sim, obs).statistic - 1 
+
+    return 1 - np.sqrt(eMean**2 + eVar**2 + eCor**2)

tank_core/evapotranspiration.py

@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 '''
-    Reference evapotranspiration calculatin models
+    Reference evapotranspiration calculation models
 
     References: 
     1. Evaluation of alternative methods for estimating reference evapotranspiration

tank_core/global_config.py

@@ -120,7 +120,7 @@
     tank_param_bounds['t3_soc']['max'],
 ])
 
-# Chanel - MUSKINGUM
+# Channel - MUSKINGUM
 
 muskingum_param_bound:dict = {
     "k" : {"min": 0, "max": 5},

tank_core/io_helpers.py

@@ -2,13 +2,14 @@
 '''
 Necessary file i/o helper functions 
 - for reading files
-- checking time consistancy of input files
+- checking time consistency of input files
 '''
 import numpy as np
 import pandas as pd 
 import json
 import os
 from datetime import datetime as dt
+from pathlib import Path
 
 from . import global_config as gc
 
@@ -43,6 +44,9 @@ def write_ts_file(df:pd.DataFrame,file_path:str)->None:
     writes model input/output timeseries files (precip, et, discharge, result etc.)
     returns 
     '''
+    # abort if index name is not Time
+    if df.index.name != 'Time':
+        raise Exception('Error: invalid time-series data no Time index found')
 
     status  = df.to_csv(
         file_path,
@@ -54,7 +58,34 @@ def write_ts_file(df:pd.DataFrame,file_path:str)->None:
 
     return status
 
-def read_project_file(project_file:str)->dict:
+def check_project(project:dict, project_dir:Path, check_discharge_file:bool)->tuple:
+
+    input_keys = ["basin", "precipitation", "evapotranspiration" ]
+    if check_discharge_file: input_keys.append('discharge')
+
+    # check if mandatory keys are present in the project definition
+    mandatory_keys = ["interval" , *input_keys]
+
+    for k in project.keys():
+        if k not in mandatory_keys:
+            return (False, f'Missing mandatory field {k} in project file')
+
+    # check if time interval is okay
+    interval_checks = [
+        project['interval'] in [0.25, 0.5], 
+        project['interval'].is_integer() 
+    ]
+    if True not in interval_checks:
+        return (False, f'invalid interval {k} hr')
+
+    # check if files of input_keys are present
+    for k in input_keys:
+        if not os.path.exists(project_dir /  project[k]):
+            (False, f'no file for {k} found in project directory')
+
+    return (True, 'All checks passed')
+
+def read_project_file(project_file:str, check_discharge_file=False)->dict:
 
     if not os.path.exists(project_file):
         raise Exception('provided project file doesn\'t exists')
@@ -64,14 +95,15 @@ def read_project_file(project_file:str)->dict:
         project = json.load(pfrb)
 
         # check if project file is okay
+        project_dir = Path(project_file).resolve().parent
+        check_ok, msg = check_project(project, project_dir, check_discharge_file)
 
-        # check if basin file is  okay
-
+        if check_ok == False:
+            raise Exception(msg)
 
         return project 
 
-    return None
-
+
 def read_basin_file(basin_file:str)->dict:
 
     if not os.path.exists(basin_file):
@@ -81,11 +113,11 @@ def read_basin_file(basin_file:str)->dict:
 
         basin = json.load(basin_file_rd_buffer)
 
-        # check if basin file is  okay [will work on it later]
+        # check if basin file is  okay [will work on it later, 
+        # basically check for missing link
 
         return basin 
 
-