Commit 9d36f12e authored by Anthony Schrapffer's avatar Anthony Schrapffer
Browse files

Correction of the sub_area and output to check area

parent 82f98503
......@@ -23,12 +23,15 @@ OutGraphFile=config.get("OverAll","GraphFile")
lonint=np.array(config.get("OverAll", "DiagLon").split(","),dtype=float)
latint=np.array(config.get("OverAll", "DiagLat").split(","),dtype=float)
#
EarthRadius=config.getfloat("OverAll", "EarthRadius")
#
import ModelGrid as MG
import Partition as PA
import RPPtools as RPP
import HydroGrid as HG
import diagplots as DP
import Interface as IF
from spherical_geometry import vector
#
# Logging in MPI : https://groups.google.com/forum/#!topic/mpi4py/SaNzc8bdj6U
#
......@@ -67,7 +70,7 @@ modelgrid=MG.ModelGrid(part.ihstart+gg.igstart,part.nih,part.jhstart+gg.jgstart,
INFO("Longitude interval on proc "+str(rank)+" = "+str(modelgrid.box_land[0][0])+" : "+str(modelgrid.box_land[0][1]))
INFO("Latitude interval on proc "+str(rank)+" = "+str(modelgrid.box_land[1][0])+" : "+str(modelgrid.box_land[1][1]))
print("modelgrid.box_land",modelgrid.box_land)
hydrogrid=HG.HydroGrid(modelgrid.box_land)
icell = 0
......@@ -78,38 +81,50 @@ sub_lat = []
#
# Only compute the weights if the file does not exist.
#
if not os.path.exists(wdir+"/"+wfile) :
INFO("Computing weights")
for cell, area in zip(modelgrid.polylist, modelgrid.area) :
overlap=np.zeros((hydrogrid.jjm,hydrogrid.iim), dtype=float)
llinside=vector.vector_to_lonlat(cell.polygons[0].inside[0],cell.polygons[0].inside[1],cell.polygons[0].inside[2])
area_in=np.zeros((hydrogrid.jjm,hydrogrid.iim), dtype=float)
for hydrocell,index in zip(hydrogrid.polylist, hydrogrid.index) :
if cell.intersects_poly(hydrocell) :
overlap[index[0],index[1]] = cell.overlap(hydrocell)
#
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)
##############
# 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.
#
if np.abs(1-np.sum(overlap)) > 0.05 :
print("ERROR on overlap area of ", icell, " : ", np.sum(overlap))
print("overlap largest values : ",np.sort(overlap.flatten(),)[-6:])
jmax,imax = np.unravel_index(np.argmax(overlap), (hydrogrid.jjm,hydrogrid.iim))
"""
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(overlap) ================")
print(overlap[jmax-1:jmax+2,imax-1:imax+2])
overlap[overlap > 10] = 0.0
#
overlap[:,:]=overlap[:,:]/np.sum(overlap)
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(overlap))+
" Nb of Hydro grid overlaping :"+str(np.count_nonzero(overlap)))
sub_index.append(np.where(overlap > 0.0))
sub_area.append(np.extract(overlap > 0.0, overlap*area))
sub_lon.append(np.extract(overlap > 0.0, hydrogrid.lon))
sub_lat.append(np.extract(overlap > 0.0, hydrogrid.lat))
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))
sub_lat.append(np.extract(area_in > 0.0, hydrogrid.lat))
icell = icell + 1
#
# Save the weights to a dedicated file for future use
......@@ -152,7 +167,7 @@ IF.initatmgrid(rank, nbcore, nbpt, modelgrid)
INFO("hoverlap")
hoverlap = IF.HydroOverlap(nbpt, nbvmax, sub_pts, sub_index, sub_area, sub_lon, sub_lat, modelgrid, hydrodata)
hoverlap = IF.HydroOverlap(nbpt, nbvmax, sub_pts, sub_index, sub_area, sub_lon, sub_lat, part, modelgrid, hydrodata)
#
# Do some memory management and synchronize procs.
......
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