Commit 3acee208 authored by POLCHER Jan's avatar POLCHER Jan 🚴🏾
Browse files

Move code out of the main code to make it simpler.

parent 70d96a98
#
import numpy as np
from numba import jit
import os, sys
#
FillValue=np.nan
IntFillValue=999999
#
import pickle
from spherical_geometry import vector
#
# Configuration
#
import configparser
config=configparser.ConfigParser({'SaveWeights':'true'})
config.read("run.def")
saveweights=config.get("OverAll", "SaveWeights")
EarthRadius=config.getfloat("OverAll", "EarthRadius")
#
# Logging in MPI : https://groups.google.com/forum/#!topic/mpi4py/SaNzc8bdj6U
#
import getargs
log_master, log_world = getargs.getLogger()
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
#
# Overlap with the the lat/lon box. The function only tests if there is a potential overlap.
# The default answer is that there is an overlap, unless we find it impossible.
#
......@@ -91,5 +110,75 @@ def dumpfield(x, lon, lat, filename, varname) :
#
return
#
# Compute the weights of the overlap of modelgrif and hydrogrid
# Only if there is no file which already contains the information.
#
#
def compweights(wfile, modelgrid, hydrogrid) :
icell = 0
sub_index = []
sub_area = []
sub_lon = []
sub_lat = []
#
# Only compute the weights if the file does not exist.
#
if not os.path.exists(wfile) :
INFO("Computing weights")
for cell, center, radius, area in zip(modelgrid.polylist, modelgrid.centers, modelgrid.radius, modelgrid.area) :
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)
selected = hydrogrid.select(center, radius)
INFO(str(icell)+" Number of HydroCells selected : "+str(len(selected))+ " out of "+str(len(hydrogrid.index)))
for isel in selected :
hydrocell = hydrogrid.polylist[isel]
index = hydrogrid.index[isel]
if cell.intersects_poly(hydrocell):
inter = cell.intersection(hydrocell)
# 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
#
#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))
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
#
if ( saveweights.lower() == 'true' ) :
with open(wfile, 'wb') as fp:
pickle.dump(sub_index, fp)
pickle.dump(sub_area, fp)
pickle.dump(sub_lon, fp)
pickle.dump(sub_lat, fp)
fp.close()
else :
#
# Extract the weights from the file if it exists
#
INFO("Reading weights")
with open (wfile, 'rb') as fp:
sub_index = pickle.load(fp)
sub_area = pickle.load(fp)
sub_lon = pickle.load(fp)
sub_lat = pickle.load(fp)
fp.close()
#
#
#
return sub_lon, sub_lat, sub_index, sub_area
......@@ -3,7 +3,6 @@
#
import numpy as np
import os, sys
import pickle
from mpi4py import MPI
import gc
#
......@@ -16,16 +15,13 @@ import time
# Gert the information from the configuration file.
#
import configparser
config=configparser.ConfigParser({'SaveWeights':'true', "DiagLon":"0.0, 0.0", "DiagLat":"0.0, 0.0"})
config=configparser.ConfigParser({"DiagLon":"0.0, 0.0", "DiagLat":"0.0, 0.0"})
config.read("run.def")
saveweights=config.get("OverAll", "SaveWeights")
nbasmax=int(config.getfloat("OverAll", "nbasmax"))
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
......@@ -33,7 +29,6 @@ import HydroGrid as HG
import diagplots as DP
import Interface as IF
from Truncate import trunc as TR
from spherical_geometry import vector
#
# Logging in MPI : https://groups.google.com/forum/#!topic/mpi4py/SaNzc8bdj6U
#
......@@ -74,72 +69,10 @@ INFO("Latitude interval on proc "+str(rank)+" = "+str(modelgrid.box_land[1][0])+
print("modelgrid.box_land",modelgrid.box_land)
hydrogrid=HG.HydroGrid(modelgrid.box_land)
icell = 0
sub_index = []
sub_area = []
sub_lon = []
sub_lat = []
#
# Only compute the weights if the file does not exist.
# Computes weights of overlap of modelgrid and hydrogrid
#
if not os.path.exists(wdir+"/"+wfile) :
INFO("Computing weights")
for cell, center, radius, area in zip(modelgrid.polylist, modelgrid.centers, modelgrid.radius, modelgrid.area) :
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)
selected = hydrogrid.select(center, radius)
INFO(str(icell)+" Number of HydroCells selected : "+str(len(selected))+ " out of "+str(len(hydrogrid.index)))
for isel in selected :
hydrocell = hydrogrid.polylist[isel]
index = hydrogrid.index[isel]
if cell.intersects_poly(hydrocell):
inter = cell.intersection(hydrocell)
# 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
#
#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))
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
#
if ( saveweights.lower() == 'true' ) :
with open(wdir+"/"+wfile, 'wb') as fp:
pickle.dump(sub_index, fp)
pickle.dump(sub_area, fp)
pickle.dump(sub_lon, fp)
pickle.dump(sub_lat, fp)
fp.close()
else :
#
# Extract the weights from the file if it exists
#
INFO("Reading weights")
with open (wdir+"/"+wfile, 'rb') as fp:
sub_index = pickle.load(fp)
sub_area = pickle.load(fp)
sub_lon = pickle.load(fp)
sub_lat = pickle.load(fp)
fp.close()
sub_lon, sub_lat, sub_index, sub_area = RPP.compweights(wdir+"/"+wfile, modelgrid, hydrogrid)
#
#
#
......
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