Changes in Hydrogrid for land points outside Hydrosheds (>60°S) + some clean up

HydroGrid.py

@@ -6,6 +6,7 @@ from netCDF4 import Dataset
 # Doc : https://spacetelescope.github.io/spherical_geometry/api/spherical_geometry.polygon.SphericalPolygon.html
 from spherical_geometry import polygon
 import RPPtools as RPP
+import os 
 #
 import configparser
 config=configparser.ConfigParser()
@@ -31,35 +32,27 @@ def getbox(ncdf, corners) :
 #
 def corners(lon, lat) :
     jjm,iim = lon.shape
-    cornerspoly = []
-    cornersll = []
-    centersll=[]
-    radiusll = []
-    index = []
+    #
     hdlon=np.mean(np.abs(np.diff(lon[0,:])))
     hdlat=np.mean(np.abs(np.diff(lat[:,0])))
     #
-    for i in range(iim) :
-        for j in range(jjm) :
-            #
-            boxll = RPP.boxit([lon[j,i], lat[j,i]], hdlon, hdlat, 2)
-            lons = [p[0] for p in boxll]
-            lats = [p[1] for p in boxll]
-            #
-            cornerspoly.append(polygon.SphericalPolygon.from_lonlat(lons, lats, center=[lon[j,i], lat[j,i]]))
-            #
-            centersll.append([lon[j,i], lat[j,i]])
-            radiusll.append(RPP.maxradius([lon[j,i], lat[j,i]], lons, lats))
-            #
-            index.append([j,i])
-            cornersll.append(boxll)
-
+    cornersll = [RPP.boxit([lon[j,i], lat[j,i]], hdlon, hdlat, 2) for i in range(iim) for j in range(jjm)]
+    Llons = [[p[0] for p in boxll] for boxll in cornersll]
+    Llats = [[p[1] for p in boxll] for boxll in cornersll]
+    index = [[j,i] for i in range(iim) for j in range(jjm)]
+    centersll = [[lon[j,i], lat[j,i]] for j,i in index]
+    cornerspoly = [ polygon.SphericalPolygon.from_lonlat(lons, lats, center=cent) for lons, lats, cent in zip(Llons, Llats, centersll)]
+    radiusll = [RPP.maxradius(cent, lons, lats) for cent, lons, lats in zip(centersll,Llons,Llats)]
+    #
     return cornersll, cornerspoly, centersll, radiusll, index
 #
-def gather(x, index) :
+def gather(x, index, default = 0) :
     y=[]
     for ia in index :
-        y.append(list(x[ia[0,i],ia[1,i]] for i in range(ia.shape[1]) ))
+        if (ia[:,0] == [-1, -1]).all():
+            y.append([default])
+        else:
+            y.append(list(x[ia[0,i],ia[1,i]] for i in range(ia.shape[1]) ))
     return y
 #
 def getattrcontaining(nc, varname, substr) :
@@ -71,7 +64,7 @@ def getattrcontaining(nc, varname, substr) :
     return att
 #
 class HydroGrid :
-    def __init__(self, lolacorners) :
+    def __init__(self, lolacorners, wfile) :
         #
         self.source=config.get("OverAll", "HydroFile")
         INFO("Opening in HydroGrid : "+self.source)
@@ -83,31 +76,32 @@ class HydroGrid :
         self.jjm,self.iim = self.lon.shape
         DEBUG("# Range Lon :"+str(np.min(self.lon))+" -- "+str(np.max(self.lon)))
         DEBUG("# Range Lat :"+str(np.min(self.lat))+" -- "+str(np.max(self.lat)))
-        self.polyll, self.polylist, self.centers, self.radius, self.index = corners(self.lon, self.lat)
+        #
+        if not os.path.exists(wfile): 
+            self.polyll, self.polylist, self.centers, self.radius, self.index = corners(self.lon, self.lat)
 
     def select(self, c, r) :
         indices=[]
-        for i in range(len(self.centers)) :
-            if RPP.loladist(c,self.centers[i]) <= r+self.radius[i] :
-                indices.append(i)
+        indices = [i for i in range(len(self.centers)) if (RPP.loladist(c,self.centers[i]) <= r+self.radius[i])]
         return indices
 
 class HydroData :
     def __init__(self, nf, box, index) :
         istr, iend, jstr, jend = box[:]
-        self.trip=gather(nf.variables["trip"][jstr:jend,istr:iend].astype(np.float32), index)
+        self.trip=gather(nf.variables["trip"][jstr:jend,istr:iend].astype(np.float32), index, 97)
         self.tripdesc=nf.variables["trip"].long_name
         #
-        self.basins=gather(nf.variables["basins"][jstr:jend,istr:iend], index)
+        self.basins=gather(nf.variables["basins"][jstr:jend,istr:iend], index, 999)
         self.basinsdesc=nf.variables["basins"].long_name
         att = getattrcontaining(nf, "basins", "max")
         if len(att) > 0 :
             self.basinsmax=att[0]
         else :
             INFO("We need to scan full file to find maximum number of basins")
-            self.basinsmax=np.max(np.where(nf.variables["basins"][:,:] <  1.e10))
+            # This variable seems not to be used further
+            self.basinsmax = part.domainmax(np.ma.max(ma.masked_where(self.basins < 1.e10, self.basins)))
         #
-        self.topoind=gather(nf.variables["topoind"][jstr:jend,istr:iend].astype(np.float32), index)
+        self.topoind=gather(nf.variables["topoind"][jstr:jend,istr:iend].astype(np.float32), index, 10)
         self.topoinddesc=nf.variables["topoind"].long_name
         att = getattrcontaining(nf, "topoind", "min")
         if len(att) > 0 :
@@ -117,7 +111,7 @@ class HydroData :
             self.topoindmin=np.min(np.where(nf.variables["topoind"][:,:] <  1.e15))
         #
         #
-        self.topoindh=gather(nf.variables["topoind_h"][jstr:jend,istr:iend].astype(np.float32), index)
+        self.topoindh=gather(nf.variables["topoind_h"][jstr:jend,istr:iend].astype(np.float32), index, 10)
         self.topoindhdesc=nf.variables["topoind_h"].long_name
         att = getattrcontaining(nf, "topoind_h", "min")
         if len(att) > 0 :
@@ -126,20 +120,20 @@ class HydroData :
             INFO("We need to scan full file to find minimum topoind_h over domain")
             self.topoindhmin=np.min(np.where(nf.variables["topoind_h"][:,:] <  1.e15))
         #
-        self.disto=gather(nf.variables["disto"][jstr:jend,istr:iend].astype(np.float32), index)
+        self.disto=gather(nf.variables["disto"][jstr:jend,istr:iend].astype(np.float32), index, 0)
         self.distodesc=nf.variables["disto"].long_name
         #
-        self.fac=gather(nf.variables["fac"][jstr:jend,istr:iend].astype(np.float32), index)
+        self.fac=gather(nf.variables["fac"][jstr:jend,istr:iend].astype(np.float32), index, 0)
         self.facdesc=nf.variables["fac"].long_name
         #
         if "orog" in nf.variables.keys():
-            self.orog = gather(nf.variables["orog"][jstr:jend,istr:iend].astype(np.float32), index)
+            self.orog = gather(nf.variables["orog"][jstr:jend,istr:iend].astype(np.float32), index, 0)
             self.orogdesc=nf.variables["orog"].long_name
         else:
             self.orog = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)
         #
         if "floodplains" in nf.variables.keys():
-            self.floodplains = gather(nf.variables["floodplains"][jstr:jend,istr:iend].astype(np.float32), index)
+            self.floodplains = gather(nf.variables["floodplains"][jstr:jend,istr:iend].astype(np.float32), index, 0)
             self.floodplainsdesc=nf.variables["floodplains"].long_name
         else:
             self.floodplains = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)