Included the special case where the intersection polygone is of length zero.

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))