diff --git a/Interface.py b/Interface.py
index 71c0e66d1ec9527dd88e0210166871fdc1fda0ac..086a88a22399bfb0ef99b88c84628500d0293213 100644
--- a/Interface.py
+++ b/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
diff --git a/RPPtools.py b/RPPtools.py
index 9fc422916467c3df8c18df5f927415ad03053b60..cf0cf23d3d800f50349da09eaf4832821bb66d7a 100644
--- a/RPPtools.py
+++ b/RPPtools.py
@@ -1,10 +1,29 @@
#
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
diff --git a/RoutingPreProc.py b/RoutingPreProc.py
index dd8a0dff14da41176c00c3130a12830768004e01..dc99f8fb208163ca70bdf86bcd75e5c07d286bd1 100644
--- a/RoutingPreProc.py
+++ b/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)