From d26e421f5bd18f13dd5f19c567354b4e26de8a60 Mon Sep 17 00:00:00 2001
From: Anthony <anthony.schrapffer@polytechnique.fr>
Date: Thu, 27 Feb 2020 13:52:40 +0100
Subject: [PATCH] Clean up in RoutingPreProc
---
RoutingPreProc.py | 24 +++++-------------------
1 file changed, 5 insertions(+), 19 deletions(-)
diff --git a/RoutingPreProc.py b/RoutingPreProc.py
index 3fded15..9513a1d 100644
--- a/RoutingPreProc.py
+++ b/RoutingPreProc.py
@@ -97,28 +97,14 @@ if not os.path.exists(wdir+"/"+wfile) :
area_in[index[0],index[1]] = inter.area()*(EarthRadius**2)
else:
area_in[index[0],index[1]] = (4*np.pi-inter.area())*(EarthRadius**2)
-
+
##############
- # Because of an issue in spherical geometry we need to exclude points where
- # the area of the intersecting polygone is not correctly computes.
- # This is temporary while the spherical geometru code is corrected.
- # We found that it only occurs on hydrogrids with very small overlap.
+ # Output of Overlap areas for each grid point
#
- """
- if np.abs((area-np.sum(area_in))/area) > 0.05 :
- print("ERROR on overlap area of ", llinside, " : ", np.sum(area_in)/area)
- print("overlap area largest values : ",np.sort(area_in.flatten(),)[-6:])
- jmax,imax = np.unravel_index(np.argmax(area_in), (hydrogrid.jjm,hydrogrid.iim))
- print("Coordinates of Max : ", hydrogrid.lon[jmax,imax], hydrogrid.lat[jmax,imax])
- print("============ Values around max(area_in) ================")
- print(area_in[jmax-1:jmax+2,imax-1:imax+2])
- area_in[area_in > 10*(EarthRadius**2)] = 0.0
- """
- ##############
-
- 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")
+ #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))
--
GitLab