diff --git a/RoutingPreProc.py b/RoutingPreProc.py index c1489d7fea1dd160cdc88a923f42b77f29b3b61e..896f4311175a180497482381b73a99915318c553 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))