diff --git a/RoutingPreProc.py b/RoutingPreProc.py index 033ed5e7afa726db85a776f636b89505b45ca551..3fded15df050180f817028bf16c42af5b11a1d0c 100644 --- a/RoutingPreProc.py +++ b/RoutingPreProc.py @@ -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.