Merge branch 'master' of https://gitlab.in2p3.fr/jan.polcher/routingpp

Interface.py

@@ -6,6 +6,7 @@ import pickle
 from netCDF4 import Dataset
 import RPPtools as RPP
 from mpi4py import MPI
+import gc
 #
 import sys
 from inspect import currentframe, getframeinfo

RPPtools.py

 #
 import numpy as np
 from numba import jit
+import os, sys
 #
 FillValue=np.nan
 IntFillValue=999999
 #
+import pickle
+from spherical_geometry import vector
+#
+# Configuration
+#
+import configparser
+config=configparser.ConfigParser({'SaveWeights':'true'})
+config.read("run.def")
+saveweights=config.get("OverAll", "SaveWeights")
+EarthRadius=config.getfloat("OverAll", "EarthRadius")
+#
+# Logging in MPI : https://groups.google.com/forum/#!topic/mpi4py/SaNzc8bdj6U
+#
+import getargs
+log_master, log_world = getargs.getLogger()
+INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
+INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
+#
 # Overlap with the the lat/lon box. The function only tests if there is a potential overlap.
 # The default answer is that there is an overlap, unless we find it impossible.
 #
@@ -91,5 +110,75 @@ def dumpfield(x, lon, lat, filename, varname) :
     #
     return
 #
+# Compute the weights of the overlap of modelgrif and hydrogrid
+# Only if there is no file which already contains the information.
 #
-#
+def compweights(wfile, modelgrid, hydrogrid) :
+    icell = 0
+    sub_index = []
+    sub_area = []
+    sub_lon = []
+    sub_lat = []
+    #
+    # Only compute the weights if the file does not exist.
+    #
+    if not os.path.exists(wfile) :
+        INFO("Computing weights")
+        for cell, center, radius, area in zip(modelgrid.polylist, modelgrid.centers, modelgrid.radius, modelgrid.area) :
+            llinside=vector.vector_to_lonlat(cell.polygons[0].inside[0],cell.polygons[0].inside[1],cell.polygons[0].inside[2])
+            area_in=np.zeros((hydrogrid.jjm,hydrogrid.iim), dtype=float)
+
+            selected = hydrogrid.select(center, radius)
+            INFO(str(icell)+" Number of HydroCells selected : "+str(len(selected))+ " out of "+str(len(hydrogrid.index)))
+            for isel in selected :
+                hydrocell = hydrogrid.polylist[isel]
+                index = hydrogrid.index[isel]
+                if cell.intersects_poly(hydrocell):
+                    inter = cell.intersection(hydrocell)
+                    # Another strange behaviour of SphericalHarmonics. There should be an overlap but the intersection polygone is empty.
+                    if len(inter.polygons) > 0 :
+                        inside = inter.polygons[0]._inside
+                        if hydrocell.contains_point(inside) and cell.contains_point(inside):
+                            area_in[index[0],index[1]] = inter.area()*(EarthRadius**2)
+                        else:
+                            area_in[index[0],index[1]] = (4*np.pi-inter.area())*(EarthRadius**2)
+        
+            ##############
+            # Output of Overlap areas for each grid point
+            #
+            #ncfile="Area_"+str("%.2f"%llinside[0])+"_"+str("%.2f"%llinside[1])+"_"+str(rank+1)+"of"+str(nbcore)+".nc"
+            #RPP.dumpfield(area_in, hydrogrid.lon, hydrogrid.lat, ncfile, "Intersect area")
+            #
+
+            INFO(str(icell)+" Sum of overlap "+str(np.sum(area_in)/area)+
+                 " Nb of Hydro grid overlaping : "+str(np.count_nonzero(area_in)))
+            sub_index.append(np.where(area_in > 0.0))
+            sub_area.append(np.extract(area_in > 0.0, area_in))
+            sub_lon.append(np.extract(area_in > 0.0, hydrogrid.lon))
+            sub_lat.append(np.extract(area_in > 0.0, hydrogrid.lat))
+            icell = icell + 1
+        #
+        # Save the weights to a dedicated file for future use
+        #
+        if ( saveweights.lower() == 'true' ) :
+            with open(wfile, 'wb') as fp:
+                pickle.dump(sub_index, fp)
+                pickle.dump(sub_area, fp)
+                pickle.dump(sub_lon, fp)
+                pickle.dump(sub_lat, fp)
+            fp.close()
+    else :
+        #
+        # Extract the weights from the file if it exists
+        #
+        INFO("Reading weights")
+        with open (wfile, 'rb') as fp:
+            sub_index = pickle.load(fp)
+            sub_area = pickle.load(fp)
+            sub_lon = pickle.load(fp)
+            sub_lat = pickle.load(fp)
+        fp.close()
+    #
+    #
+    #
+    return sub_lon, sub_lat, sub_index, sub_area

RoutingPreProc.py

@@ -3,7 +3,6 @@
 #
 import numpy as np
 import os, sys
-import pickle
 from mpi4py import MPI
 import gc
 #
@@ -16,17 +15,15 @@ import time
 # Gert the information from the configuration file.
 #
 import configparser
-config=configparser.ConfigParser({'SaveWeights':'true', "DiagLon":"0.0, 0.0", "DiagLat":"0.0, 0.0", "numop" : 100})
+config=configparser.ConfigParser({"DiagLon":"0.0, 0.0", "DiagLat":"0.0, 0.0", "numop":'100'})
 config.read("run.def")
-saveweights=config.get("OverAll", "SaveWeights")
-nbasmax=int(config.getfloat("OverAll", "nbasmax"))
-numop=int(config.getfloat("OverAll", "numop"))
+nbasmax=config.getint("OverAll", "nbasmax")
+numop=config.getint("OverAll", "numop")
 OutGraphFile=config.get("OverAll","GraphFile")
+DumpHydroSuper=config.getboolean("OverAll","DumpHydroSuper",fallback=False)
 lonint=np.array(config.get("OverAll", "DiagLon").split(","),dtype=float)
 latint=np.array(config.get("OverAll", "DiagLat").split(","),dtype=float)
 #
-EarthRadius=config.getfloat("OverAll", "EarthRadius")
-#
 import ModelGrid as MG
 import Partition as PA
 import RPPtools as RPP
@@ -34,7 +31,6 @@ import HydroGrid as HG
 import diagplots as DP
 import Interface as IF
 from Truncate import trunc as TR
-from spherical_geometry import vector
 #
 # Logging in MPI : https://groups.google.com/forum/#!topic/mpi4py/SaNzc8bdj6U
 #
@@ -75,72 +71,10 @@ INFO("Latitude interval on proc "+str(rank)+" = "+str(modelgrid.box_land[1][0])+
 
 print("modelgrid.box_land",modelgrid.box_land)
 hydrogrid=HG.HydroGrid(modelgrid.box_land)
-
-icell = 0
-sub_index = []
-sub_area = []
-sub_lon = []
-sub_lat = []
 #
-# Only compute the weights if the file does not exist.
+# Computes weights of overlap of modelgrid and hydrogrid
 #
-
-if not os.path.exists(wdir+"/"+wfile) :
-    INFO("Computing weights")
-    for cell, center, radius, area in zip(modelgrid.polylist, modelgrid.centers, modelgrid.radius, modelgrid.area) :
-        llinside=vector.vector_to_lonlat(cell.polygons[0].inside[0],cell.polygons[0].inside[1],cell.polygons[0].inside[2])
-        area_in=np.zeros((hydrogrid.jjm,hydrogrid.iim), dtype=float)
-
-        selected = hydrogrid.select(center, radius)
-        INFO(str(icell)+" Number of HydroCells selected : "+str(len(selected))+ " out of "+str(len(hydrogrid.index)))
-        for isel in selected :
-            hydrocell = hydrogrid.polylist[isel]
-            index = hydrogrid.index[isel]
-            if cell.intersects_poly(hydrocell):
-                inter = cell.intersection(hydrocell)
-                # Another strange behaviour of SphericalHarmonics. There should be an overlap but the intersection polygone is empty.
-                if len(inter.polygons) > 0 :
-                    inside = inter.polygons[0]._inside
-                    if hydrocell.contains_point(inside) and cell.contains_point(inside):
-                        area_in[index[0],index[1]] = inter.area()*(EarthRadius**2)
-                    else:
-                        area_in[index[0],index[1]] = (4*np.pi-inter.area())*(EarthRadius**2)
-        
-        ##############
-        # Output of Overlap areas for each grid point
-        #
-        #ncfile="Area_"+str("%.2f"%llinside[0])+"_"+str("%.2f"%llinside[1])+"_"+str(rank+1)+"of"+str(nbcore)+".nc"
-        #RPP.dumpfield(area_in, hydrogrid.lon, hydrogrid.lat, ncfile, "Intersect area")
-        #
-
-        INFO(str(icell)+" Sum of overlap "+str(np.sum(area_in)/area)+
-             " Nb of Hydro grid overlaping : "+str(np.count_nonzero(area_in)))
-        sub_index.append(np.where(area_in > 0.0))
-        sub_area.append(np.extract(area_in > 0.0, area_in))
-        sub_lon.append(np.extract(area_in > 0.0, hydrogrid.lon))
-        sub_lat.append(np.extract(area_in > 0.0, hydrogrid.lat))
-        icell = icell + 1
-    #
-    # Save the weights to a dedicated file for future use
-    #
-    if ( saveweights.lower() == 'true' ) :
-        with open(wdir+"/"+wfile, 'wb') as fp:
-            pickle.dump(sub_index, fp)
-            pickle.dump(sub_area, fp)
-            pickle.dump(sub_lon, fp)
-            pickle.dump(sub_lat, fp)
-        fp.close()
-else :
-    #
-    # Extract the weights from the file if it exists
-    #
-    INFO("Reading weights")
-    with open (wdir+"/"+wfile, 'rb') as fp:
-        sub_index = pickle.load(fp)
-        sub_area = pickle.load(fp)
-        sub_lon = pickle.load(fp)
-        sub_lat = pickle.load(fp)
-    fp.close()
+sub_lon, sub_lat, sub_index, sub_area = RPP.compweights(wdir+"/"+wfile, modelgrid, hydrogrid)
 #
 #
 #
@@ -205,7 +139,9 @@ t1 = time.time()
 print("Time for fetch: {:0.2f} s.".format(t1-t)) 
 comm.Barrier()
 
-hsuper.dumpnetcdf(OutGraphFile.replace(".nc","_HydroSuper.nc"), gg, modelgrid, part)
+if DumpHydroSuper :
+    print("Dumping HydroSuper")
+    hsuper.dumpnetcdf(OutGraphFile.replace(".nc","_HydroSuper.nc"), gg, modelgrid, part)
 
 print("Rank : {0} - Basin_count Before Truncate : ".format(part.rank), hsuper.basin_count)
 hs = TR(hsuper, part, comm, modelgrid, numop = numop)