Commit d26e421f authored by Anthony's avatar Anthony
Browse files

Clean up in RoutingPreProc

parent 09fbfb17
...@@ -97,28 +97,14 @@ if not os.path.exists(wdir+"/"+wfile) : ...@@ -97,28 +97,14 @@ if not os.path.exists(wdir+"/"+wfile) :
area_in[index[0],index[1]] = inter.area()*(EarthRadius**2) area_in[index[0],index[1]] = inter.area()*(EarthRadius**2)
else: else:
area_in[index[0],index[1]] = (4*np.pi-inter.area())*(EarthRadius**2) 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 # Output of Overlap areas for each grid point
# 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.
# #
""" #ncfile="Area_"+str("%.2f"%llinside[0])+"_"+str("%.2f"%llinside[1])+"_"+str(rank+1)+"of"+str(nbcore)+".nc"
if np.abs((area-np.sum(area_in))/area) > 0.05 : #RPP.dumpfield(area_in, hydrogrid.lon, hydrogrid.lat, ncfile, "Intersect area")
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")
# #
INFO(str(icell)+" Sum of overlap"+str(np.sum(area_in)/area)+ INFO(str(icell)+" Sum of overlap"+str(np.sum(area_in)/area)+
" Nb of Hydro grid overlaping :"+str(np.count_nonzero(area_in))) " Nb of Hydro grid overlaping :"+str(np.count_nonzero(area_in)))
sub_index.append(np.where(area_in > 0.0)) sub_index.append(np.where(area_in > 0.0))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment