diff --git a/RoutingPreProc.py b/RoutingPreProc.py index 3fded15df050180f817028bf16c42af5b11a1d0c..9513a1d2009aa0f54e19054df9305d2a5b71f5c8 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))