Skip to content

GitLab

Commit c2b62560 authored by Anthony's avatar Anthony
Browse files

Merge branch 'master' of https://gitlab.in2p3.fr/ipsl/lmd/intro/routingpp

parents fe9dce1f ff23aa60
Hide whitespace changes
Inline Side-by-side
Showing with 90 additions and 26 deletions
GraphHydro.py
View file @ c2b62560
......@@ -15,7 +15,7 @@ MaxDistErr = config.getfloat("OverAll", "MaxDistErr", fallback=25.0)
# Maximum error in the upstream area in %
MaxUpstrErr = config.getfloat("OverAll", "MaxUpstrErr", fallback=10.0)
#
log_master, log_world = getargs.getLogger()
log_master, log_world = getargs.getLogger(__name__)
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
#
......@@ -116,8 +116,25 @@ class HydroGraph :
ij, ib = np.unravel_index(np.nanargmax(cost),cost.shape)
locations.addhtupos(i, ij, ib, np.nanmax(cost))
else :
INFO("Could not locate station number : "+str(lid)+" of type "+locations.ltype[i]+" at "+locations.lname[i])
#
# Only report a missed station is there was a geographical solution on the domain.
if np.count_nonzero(~np.isnan(disterr)) > 0 :
if lid < 9999 :
typename="Structure"
else :
typename="GRDC"
#
INFO("Could not locate "+typename+" number : "+str(lid)+" of type "+\
locations.ltype[i]+" at "+locations.lname[i])
ij = np.nanargmax(disterr)
INFO(typename+" number : "+str(lid)+" best distance match [km] : "+str(dist[ij]/1000.))
if np.count_nonzero(~np.isnan(upstrerr)) > 0 :
ib, ij = np.unravel_index(np.nanargmax(upstrerr),upstrerr.shape)
INFO(typename+" number : "+str(lid)+" smallest upstream error [%] : "+\
str(upstrerr[ib,ij]*MaxUpstrErr+MaxUpstrErr))
else :
INFO(typename+" number : "+str(lid)+" No appropriate HTU with an upstream "+\
"error less than "+str(MaxUpstrErr)+" % ")
#
#
#
def add_variable(self, outnf, procgrid, NCFillValue, part, coord, name, title, units, data, vtyp, orig_type = "float", arrayorder = 'F'):
......@@ -147,7 +164,7 @@ class HydroGraph :
cornerind=[0,2,4,6]
nbcorners = len(cornerind)
#
nlmax = locations.maxlocations(self.nbasmax, part)
nlmax, nblocated = locations.maxlocations(self.nbasmax, part)
#
if part.rank == 0 :
outnf=Dataset(filename, 'w', format='NETCDF4_CLASSIC')
......@@ -158,7 +175,8 @@ class HydroGraph :
outnf.createDimension('z', self.nbasmax)
outnf.createDimension('bnd', nbcorners)
outnf.createDimension('inflow', self.max_inflow)
outnf.createDimension('stations', nlmax)
outnf.createDimension('stnperhtu', nlmax)
outnf.createDimension('locations', nblocated)
else :
outnf = None
#
......@@ -243,7 +261,7 @@ class HydroGraph :
#
# Build and save a map of all the structures located on the graph
#
locations.addtonetcdf(self.nbpt, self.nbasmax, outnf, procgrid, NCFillValue, part, ('stations','z','y','x'), vtyp)
locations.addtonetcdf(self.nbpt, self.nbasmax, outnf, procgrid, part, ('y','x'), ('locations',), vtyp)
#
# Close the file
#
......
Locations.py
View file @ c2b62560
......@@ -8,7 +8,7 @@ import codecs
import getargs
config = getargs.SetupConfig()
#
log_master, log_world = getargs.getLogger()
log_master, log_world = getargs.getLogger(__name__)
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
#
......@@ -122,29 +122,40 @@ class Locations :
u,n = np.unique(ghl[~np.isnan(ghl)], return_counts=True)
if len(u) > 0 :
nlmax = max(n)
# Maximum number of HTU to monitor per grid
iu = [int(ii) for ii in u]
uhtu,nhtu = np.unique(iu, return_counts=True)
htumonmax=max(nhtu)
else :
nlmax = 1
htumonmax = 1
else :
nlmax = None
htumonmax = None
nlmax = part.bcast(nlmax)
return nblocated, nlmax
htumonmax = part.bcast(htumonmax)
return nblocated, nlmax, htumonmax
#
# Function to get the maximum number of locations per HTU
#
def maxlocations(self, nbasmax, part) :
self.cleanup()
self.nblocated, self.nlmax = self.unique(nbasmax, part)
self.nblocated, self.nlmax, self.htumonmax = self.unique(nbasmax, part)
INFO("Out of a total of "+str(self.gnbloc)+" stations "+str(self.nblocated)+" could be located on the HTU graph.")
INFO("Maximum number of stations per HTU : "+str(self.nlmax))
return self.nlmax
INFO("Maximum number of HTU to monitor per grid : "+str(self.htumonmax))
return self.nlmax, self.nblocated
#
# Function to place all IDs of structures onto a map and remove duplicated locations
#
def addtonetcdf(self, nbpt, nbasmax, outnf, procgrid, NCFillValue, part, coord, vtyp) :
def addtonetcdf(self, nbpt, nbasmax, outnf, procgrid, part, spacecoord, loccoord, vtyp) :
#
name = "locations"
title = "Locations of stations or structures"
units = "-"
locname = "locations"
loctitle = "Basic information of located stations or structures"
locunits = "-"
monname = "HTUmonitor"
montitle = "Index of HTU to be monitored in the model"
monunits = "-"
#
locmap = np.zeros((nbpt,nbasmax,self.nlmax), dtype=vtyp, order='F')
locmap[:,:,:] = np.nan
......@@ -161,6 +172,11 @@ class Locations :
gl = part.gather(procgrid.landscatter(locmap, order='F'))
gc = part.gather(procgrid.landscatter(costmap, order='F'))
#
htumon = np.zeros((self.htumonmax,part.njg,part.nig), dtype=int)
htumon[:,:,:] = -1
monnb = np.zeros((part.njg,part.nig), dtype=int)
locationlist = np.zeros((5,self.nblocated), dtype=int)
#
# Eliminate locations with lower cost
#
if part.rank == 0:
......@@ -170,16 +186,46 @@ class Locations :
maxcost = np.max(gc[gl == lid])
gl[(gl == lid) & (gc < maxcost)] = np.nan
gc[(gl == lid) & (gc < maxcost)] = np.nan
#
# Build the list of all IDs which are now unique
#
u = np.unique(gl[~np.isnan(gl)])
for i,uid in enumerate(u) :
f = np.where(gl == uid)
if np.any(htumon[:,f[2][0],f[3][0]] == f[1][0]+1) :
# HTU already monitored
imom = np.where(htumon[:,f[2][0],f[3][0]] == f[1][0]+1)[0][0]
else :
# Add HTU to those monitored
imon = monnb[f[2][0],f[3][0]]
htumon[imon,f[2][0],f[3][0]] = f[1][0]+1
monnb[f[2][0],f[3][0]] += 1
# Information saved.
locationlist[0,i] = uid
locationlist[1,i] = f[3][0]+1
locationlist[2,i] = f[2][0]+1
locationlist[3,i] = f[1][0]+1
locationlist[4,i] = imon+1
#
# Add to NetCDF file
#
if part.rank == 0:
ncvar = outnf.createVariable(name, vtyp, coord, fill_value=NCFillValue)
ncvar.title = title
ncvar.units = units
gl[np.isnan(gl)] = NCFillValue
else :
ncvar = np.zeros((1,)*len(coord))
ncvar[:,:,:,:] = gl
# Create new dimension
monmax = np.max(monnb)
moncoord = ('htumon',)
outnf.createDimension(moncoord[0], monmax)
infocoord = ('locinfo',)
outnf.createDimension(infocoord[0], 5)
# Create variable for the list of stations
locvar = outnf.createVariable(locname, vtyp, infocoord+loccoord)
locvar.title = loctitle
locvar.units = locunits
locvar.description = "ID, i-index, j-index, HTU-index, monitoring-index"
locvar[:,:] = locationlist
# Create variable for the indices of the HTU to monitor
monvar = outnf.createVariable(monname, vtyp, moncoord+spacecoord)
monvar.title = montitle
monvar.units= monunits
monvar[:,:,:] = htumon[0:monmax,:,:]
return
#
RPPtools.py
View file @ c2b62560
......@@ -24,7 +24,7 @@ EarthRadius=config.getfloat("OverAll", "EarthRadius", fallback=6370000.0)
# Logging in MPI : https://groups.google.com/forum/#!topic/mpi4py/SaNzc8bdj6U
#
import getargs
log_master, log_world = getargs.getLogger()
log_master, log_world = getargs.getLogger(__name__)
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
#
......
Truncate.py
View file @ c2b62560
......@@ -17,7 +17,7 @@ localdir=os.path.dirname(getframeinfo(currentframe()).filename)
sys.path.append(localdir+'/F90subroutines')
import routing_interface
import getargs
log_master, log_world = getargs.getLogger()
log_master, log_world = getargs.getLogger(__name__)
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
######################
......
WeightsOnly.py
View file @ c2b62560
......@@ -31,7 +31,7 @@ import HydroGrid as HG
# Logging in MPI : https://groups.google.com/forum/#!topic/mpi4py/SaNzc8bdj6U
#
import getargs
log_master, log_world = getargs.getLogger()
log_master, log_world = getargs.getLogger(__name__)
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
#
......
tests/Iberia/run_regular.def
View file @ c2b62560
......@@ -6,7 +6,7 @@ EarthRadius = 6370000.
#
ModelGridFile = /bdd/MEDI/workspaces/polcher/NewRouting/EM_WFDEI_CRU_2000.nc
WEST_EAST = -9.75, 5.25
SOUTH_NORTH = 35.5, 44.5
SOUTH_NORTH = 35.5, 43.5
HydroFile = /bdd/ORCHIDEE_Forcing/Routing/Hydro4ORCH/HydroSHEDS_MED.nc
#
# FORTRAN interface parameters
......@@ -15,7 +15,7 @@ Documentation = true
#
# Configuration for the graph to be generated
#
nbasmax = 35
nbasmax = 32
#
# Number of operation of simplification performed together
#
......
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