From 6f440c1a8b676101739e8b94ecbd1cb54257ecef Mon Sep 17 00:00:00 2001
From: POLCHER Jan <jan.polcher@lmd.jussieu.fr>
Date: Fri, 20 Mar 2020 10:36:03 +0100
Subject: [PATCH] Included the special case where the intersection polygone is
of length zero.
---
RoutingPreProc.py | 17 +++++++++--------
1 file changed, 9 insertions(+), 8 deletions(-)
diff --git a/RoutingPreProc.py b/RoutingPreProc.py
index c1489d7..896f431 100644
--- a/RoutingPreProc.py
+++ b/RoutingPreProc.py
@@ -96,12 +96,13 @@ if not os.path.exists(wdir+"/"+wfile) :
index = hydrogrid.index[isel]
if cell.intersects_poly(hydrocell):
inter = cell.intersection(hydrocell)
- 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)
+ # 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
@@ -110,8 +111,8 @@ if not os.path.exists(wdir+"/"+wfile) :
#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)))
+ 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))
--
GitLab