diff --git a/RoutingPreProc.py b/RoutingPreProc.py index e358847eabf431ec330eb79c7fa16dc36434e647..849f7f11dde61e1aaf0d7857ad7db7d68d6037eb 100644 --- a/RoutingPreProc.py +++ b/RoutingPreProc.py @@ -97,12 +97,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 @@ -111,8 +112,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))