From 6ab622bb917f14127f72a0ded04bb8820b0cda48 Mon Sep 17 00:00:00 2001
From: Peter Burek <burek@iiasa.ac.at>
Date: Fri, 8 Oct 2021 14:39:54 +0200
Subject: [PATCH] fix: some small fixes to read tif files changes: made CWaTM
 to run as library inside calibration tool

---
 cwatm/hydrological_modules/readmeteo.py   |  2 +-
 cwatm/management_modules/data_handling.py | 20 ++++++++++++--------
 2 files changed, 13 insertions(+), 9 deletions(-)

diff --git a/cwatm/hydrological_modules/readmeteo.py b/cwatm/hydrological_modules/readmeteo.py
index 7b732b70..4bef65e6 100644
--- a/cwatm/hydrological_modules/readmeteo.py
+++ b/cwatm/hydrological_modules/readmeteo.py
@@ -5,7 +5,7 @@
 # Author:      PB
 #
 # Created:     13/07/2016
-# Copyright:   (c) PB 2016
+# Copyright:   (c) PB 201664
 # -------------------------------------------------------------------------
 
 from cwatm.management_modules.data_handling import *
diff --git a/cwatm/management_modules/data_handling.py b/cwatm/management_modules/data_handling.py
index 5e9365ee..ba9a1656 100644
--- a/cwatm/management_modules/data_handling.py
+++ b/cwatm/management_modules/data_handling.py
@@ -106,11 +106,12 @@ def setmaskmapAttr(x,y,col,row,cell):
     if (x-int(x)) != 0.:
         if abs(x - int(x)) > 1e9:
             x = 1/round(1/(x-int(x)),4) + int(x)
-        else: x = round(x,4)
+        else: x = round(x,6)
     if (y - int(y)) != 0.:
         if abs(y - int(y)) > 1e9:
             y = 1 / round(1 / (y - int(y)), 4) + int(y)
-        else: y = round(y,4)
+        else: y = round(y,6)
+    # This is still not ok! Some rounding issues still appear sometimes
 
     maskmapAttr['x'] = x
     maskmapAttr['y'] = y
@@ -694,8 +695,11 @@ def mapattrNetCDF(name, check=True):
 
     xx = maskmapAttr['x']
     yy = maskmapAttr['y']
-    cut0 = int(0.0001 + np.abs(xx - lon) * invcell)  # argmin() ??
-    cut2 = int(0.0001 + np.abs(yy - lat) * invcell)
+    #cut0 = int(0.0001 + np.abs(xx - lon) * invcell) # argmin() ??
+    #cut2 = int(0.0001 + np.abs(yy - lat) * invcell) #
+    cut0 = int(np.abs(xx + maskmapAttr['cell']/2 - lon) * invcell) # argmin() ??
+    cut2 = int(np.abs(yy - maskmapAttr['cell']/2 - lat) * invcell) #
+
 
     cut1 = cut0 + maskmapAttr['col']
     cut3 = cut2 + maskmapAttr['row']
@@ -779,7 +783,7 @@ def mapattrTiff(nf2):
     cellSize = geotransform[1]
 
     #invcell = round(1/cellSize,0)
-    if cellSize > 0:
+    if cellSize > 1:
         invcell = 1 / cellSize
     else:
         invcell = round(1/cellSize,0)
@@ -846,14 +850,14 @@ def multinetdf(meteomaps, startcheck = 'dateBegin'):
             except:
                 datediv = 1
 
-            datestart = num2date(nctime[0] ,units=nctime.units,calendar=nctime.calendar)
+            datestart = num2date(nctime[:][0] ,units=nctime.units,calendar=nctime.calendar)
 
             # sometime daily records have a strange hour to start with -> it is changed to 0:00 to haqve the same record
             datestart = datestart.replace(hour=0, minute=0)
-            dateend = num2date(nctime[-1], units=nctime.units, calendar=nctime.calendar)
+            dateend = num2date(nctime[:][-1], units=nctime.units, calendar=nctime.calendar)
 
             datestartint = int(nctime[0]) // datediv
-            dateendint = int(nctime[-1]) // datediv
+            dateendint = int(nctime[:][-1]) // datediv
 
             dateend = dateend.replace(hour=0, minute=0)
             #if dateVar['leapYear'] > 0: