From 3acee20846bf7e01733d180d5f4f3c025e1d4788 Mon Sep 17 00:00:00 2001
From: POLCHER Jan <jan.polcher@lmd.jussieu.fr>
Date: Wed, 15 Apr 2020 18:48:24 +0200
Subject: [PATCH] Move code out of the main code to make it simpler.
---
RPPtools.py | 91 ++++++++++++++++++++++++++++++++++++++++++++++-
RoutingPreProc.py | 73 ++-----------------------------------
2 files changed, 93 insertions(+), 71 deletions(-)
diff --git a/RPPtools.py b/RPPtools.py
index 9fc4229..cf0cf23 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 849f7f1..fa75507 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,16 +15,13 @@ 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"})
+config=configparser.ConfigParser({"DiagLon":"0.0, 0.0", "DiagLat":"0.0, 0.0"})
config.read("run.def")
-saveweights=config.get("OverAll", "SaveWeights")
nbasmax=int(config.getfloat("OverAll", "nbasmax"))
OutGraphFile=config.get("OverAll","GraphFile")
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
@@ -33,7 +29,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
#
@@ -74,72 +69,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)
#
#
#
--
GitLab