Skip to content

GitLab

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

A version which works halfways. Further tests are needed.

parent be249805
Hide whitespace changes
Inline Side-by-side
Showing with 188 additions and 55 deletions
GraphHydro.py
View file @ b9ab8b9b
......@@ -97,8 +97,8 @@ class HydroGraph :
#
def position(self, modelgrid, locations) :
npt, nbas = self.routing_fetch.shape
if len(locations.id) > 0 :
for i,id in enumerate(locations.id) :
if len(locations.lid) > 0 :
for i,lid in enumerate(locations.lid) :
dist = RPP.distances(np.array(modelgrid.coordll)[:,0], np.array(modelgrid.coordll)[:,1], \
locations.lon[i], locations.lat[i])
# Calculation in km
......@@ -114,9 +114,10 @@ class HydroGraph :
#
#
#
def add_variable(self, outnf, procgrid, NCFillValue, part, coord, name, title, units, data, vtyp, orig_type = "float"):
var = procgrid.landscatter(data.astype(vtyp), order='F')
def add_variable(self, outnf, procgrid, NCFillValue, part, coord, name, title, units, data, vtyp, orig_type = "float", arrayorder = 'F'):
var = procgrid.landscatter(data.astype(vtyp), order=arrayorder)
if orig_type == "float":
var[np.isnan(var)] = NCFillValue
elif orig_type == "int":
......@@ -128,18 +129,20 @@ class HydroGraph :
ncvar.title = title
ncvar.units = units
else :
ncvar = np.zeros((1,1,1))
ncvar[:] = part.gather(var)
ncvar = np.zeros((1,)*len(coord))
ncvar[:] = part.gather(var, default=NCFillValue)
#
#
#
def dumpnetcdf(self, filename, globalgrid, procgrid, part, tcst) :
def dumpnetcdf(self, filename, globalgrid, procgrid, part, tcst, locations) :
#
NCFillValue=1.0e20
vtyp=np.float64
cornerind=[0,2,4,6]
nbcorners = len(cornerind)
#
nlmax = locations.maxlocations(self.nbasmax, part)
#
if part.rank == 0 :
outnf=Dataset(filename, 'w', format='NETCDF4_CLASSIC')
# Dimensions
......@@ -149,6 +152,7 @@ class HydroGraph :
outnf.createDimension('z', self.nbasmax)
outnf.createDimension('bnd', nbcorners)
outnf.createDimension('inflow', self.max_inflow)
outnf.createDimension('stations', nlmax)
else :
outnf = None
#
......@@ -170,63 +174,71 @@ class HydroGraph :
#
#
# The field route_togrid is with indices on the local grid. That needs to be converted to the global grid.
self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, "int")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, orig_type="int")
# route to basin
self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "routetobasin", "Basin in to which the water goes", "-", self.route_tobasin[:,:], vtyp, "int")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "routetobasin", "Basin in to which the water goes", "-", self.route_tobasin[:,:], vtyp, orig_type="int")
# basin id
self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basinid", "ID of basin", "-", self.global_basinid[:,:], vtyp, "int")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basinid", "ID of basin", "-", self.global_basinid[:,:], vtyp, orig_type="int")
#
# basin area
self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_area", "area of basin", "m^2", self.routing_area[:,:], vtyp, "float")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basin_area", "area of basin", "m^2", self.routing_area[:,:], vtyp)
self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_orog_mean", "Mean orography", "m", self.routing_orog_mean[:,:], vtyp, "float")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basin_orog_mean", "Mean orography", "m", self.routing_orog_mean[:,:], vtyp)
self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_orog_min", "Min orography", "m", self.routing_orog_min[:,:], vtyp, "float")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basin_orog_min", "Min orography", "m", self.routing_orog_min[:,:], vtyp)
self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_orog_max", "Max orography", "m", self.routing_orog_max[:,:], vtyp, "float")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basin_orog_max", "Max orography", "m", self.routing_orog_max[:,:], vtyp)
self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basin_floodp", "Fraction of floodplains", "-", self.routing_floodp[:,:], vtyp, "float")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basin_floodp", "Fraction of floodplains", "-", self.routing_floodp[:,:], vtyp)
# route number into basin
self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "routenbintobas", "Number of basin into current one", "-", self.route_nbintobas[:], vtyp, "int")
self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "routenbintobas", "Number of basin into current one", "-", self.route_nbintobas[:], vtyp, orig_type="int")
#
# original number into basin
self.add_variable(outnf, procgrid, NCFillValue, part, ( 'y','x'), "originnbintobas", "Number of sub-grid basin into current one before truncation", "-", self.origin_nbintobas[:], vtyp, "int")
self.add_variable(outnf, procgrid, NCFillValue, part, ( 'y','x'), "originnbintobas", "Number of sub-grid basin into current one before truncation", "-", \
self.origin_nbintobas[:], vtyp, orig_type="int")
#
# latitude of outlet
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp, "float")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp)
# longitude of outlet
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outletlon", "Longitude of Outlet", "degrees east", self.route_outlet[:,:,1], vtyp, "float")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outletlon", "Longitude of Outlet", "degrees east", self.route_outlet[:,:,1], vtyp)
# type of outlet
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp, "float")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp)
#
# topographic index
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "topoindex", "Topographic index of the retention time", "m", self.topo_resid[:,:], vtyp, "float")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "topoindex", "Topographic index of the retention time", "m", self.topo_resid[:,:], vtyp)
#
# Inflow number
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, "int")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, orig_type="int")
#
# Inflow grid
#gingrid = part.l2glandindex(self.inflow_grid[:,:,:inflow_size])
self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'z','y','x'), "route_ingrid", "Grid from which the water flows", "-", self.route_ingrid[:,:,:], vtyp, "int")
self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow','z','y','x'), "route_ingrid", "Grid from which the water flows", "-", \
self.route_ingrid[:,:,:], vtyp, orig_type="int")
#
# Inflow basin
self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'z','y','x'), "route_inbasin", "Basin from which the water flows", "-", self.route_inbasin[:,:,:], vtyp, "int")
self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow','z','y','x'), "route_inbasin", "Basin from which the water flows", "-", \
self.route_inbasin[:,:,:], vtyp, orig_type="int")
#
# Save centre of gravity of HTU
#
# Check if it works
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "CG_lon", "Longitude of centre of gravity of HTU", "degrees east", self.routing_cg[:,:,1], vtyp, "float")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "CG_lon", "Longitude of centre of gravity of HTU", "degrees east", self.routing_cg[:,:,1], vtyp)
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "CG_lat", "Latitude of centre of gravity of HTU", "degrees north", self.routing_cg[:,:,0], vtyp, "float")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "CG_lat", "Latitude of centre of gravity of HTU", "degrees north", self.routing_cg[:,:,0], vtyp)
#
# Save the fetch of each basin
#
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp, "float")
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp)
#
# 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)
#
# Close the file
#
if part.rank == 0:
outnf.close()
#
......
Locations.py
View file @ b9ab8b9b
......@@ -2,56 +2,174 @@
#
import numpy as np
from netCDF4 import Dataset
import sys
import getargs
config = getargs.SetupConfig()
#
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
#
class Locations :
def __init__(self, bbox) :
self.id = []
self.lid = []
self.lon = []
self.lat = []
self.upstream = []
self.type = []
self.ltype = []
self.gnbloc = 0
structure_file = config.get("OverAll","StructuresFile",fallback=None)
if structure_file != None :
self.addstructures(structure_file, bbox)
self.gnbloc += self.addstructures(structure_file, bbox)
grdc_file = config.get("OverAll","GRDCFile",fallback=None)
if grdc_file != None :
self.addgrdcstations(grdc_file, bbox)
self.gnbloc += self.addgrdcstations(grdc_file, bbox)
#
self.ij = [-1]*len(self.id)
self.ib = [-1]*len(self.id)
self.cost = [np.nan]*len(self.id)
self.ij = [-1]*len(self.lid)
self.ib = [-1]*len(self.lid)
self.cost = [np.nan]*len(self.lid)
return
#
# Get the locations of infrastructures from an ASCII file
#
def addstructures(self, f, bbox) :
nb = 0
for line in open(f, 'r') :
if line.find("#") != 0 :
nb += 1
col = line.split()
lon = float(col[1]); lat = float(col[2])
if min(bbox[0]) < lon < max(bbox[0]) and min(bbox[1]) < lat < max(bbox[1]) :
self.id.append(int(col[0]))
self.lid.append(int(col[0]))
self.lon.append(lon)
self.lat.append(lat)
self.upstream.append(float(col[3]))
self.type.append(" ".join(col[4:]).lower())
return
self.ltype.append(" ".join(col[4:]).lower())
return nb
#
# Get all the locations from the GRDC database
#
def addgrdcstations(self, f, bbox) :
g = Dataset(f, 'r')
nb = len(g.dimensions['stations'])
for i, id in enumerate(g.variables["number"][:]) :
if min(bbox[0]) < g.variables["lon"][i] < max(bbox[0]) and \
min(bbox[1]) < g.variables["lat"][i] < max(bbox[1]) :
self.id.append(id)
self.lid.append(id)
self.lon.append(g.variables["lon"][i])
self.lat.append(g.variables["lat"][i])
self.upstream.append(g.variables["area"][i])
self.type.append("gauging station")
self.ltype.append("gauging station")
g.close()
return
return nb
#
# Adds position of location in the HTU graph
#
def addhtupos(self, i, ij, ib, cost) :
self.ij[i] = ij
self.ib[i] = ib
self.cost[i] = cost
return
#
# Function to keep only locations which could be placed in the HTU graph
#
def cleanup(self) :
if len(self.ij) > 0 :
self.lid = np.array(self.lid)[np.array(self.ij) >= 0]
self.lon = np.array(self.lon)[np.array(self.ij) >= 0]
self.lat = np.array(self.lat)[np.array(self.ij) >= 0]
self.upstream = np.array(self.upstream)[np.array(self.ij) >= 0]
self.ltype = np.array(self.ltype)[np.array(self.ij) >= 0]
self.ij = np.array(self.ij)[np.array(self.ij) >= 0]
self.ib = np.array(self.ib)[np.array(self.ij) >= 0]
self.cost = np.array(self.cost)[np.array(self.ij) >= 0]
else :
self.lid = np.array(self.lid)
self.lon = np.array(self.lon)
self.lat = np.array(self.lat)
self.upstream = np.array(self.upstream)
self.ltype = np.array(self.ltype)
self.ij = np.array(self.ij)
self.ib = np.array(self.ib)
self.cost = np.array(self.cost)
return
#
# Get all the unique locations identified in the HTU graph
#
def unique(self, nbasmax, part) :
# Number of stations positioned
glid = part.gatherbypoint(self.lid.astype(float))
if part.rank == 0 :
nblocated = np.unique(glid[~np.isnan(glid)]).shape[0]
else :
nblocated = None
nblocated = part.bcast(nblocated)
# Number of stations per HTU
f = 10**int(np.log10(nbasmax)+1)
htuloc = self.ij + self.ib/f
ghl = part.gatherbypoint(htuloc)
if part.rank == 0 :
u,n = np.unique(ghl[~np.isnan(ghl)], return_counts=True)
if len(u) > 0 :
nlmax = max(n)
else :
nlmax = 1
else :
nlmax = None
nlmax = part.bcast(nlmax)
return nblocated, nlmax
#
# 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)
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
#
# 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) :
#
name = "locations"
title = "Locations of stations or structures"
units = "-"
#
locmap = np.zeros((nbpt,nbasmax,self.nlmax), dtype=vtyp, order='F')
locmap[:,:,:] = np.nan
costmap = np.zeros((nbpt,nbasmax,self.nlmax), dtype=vtyp, order='F')
costmap[:,:,:] = np.nan
#
for i,lid in enumerate(self.lid) :
ip = 0
while ~np.isnan(locmap[self.ij[i],self.ib[i],ip]) :
ip += 1
locmap[self.ij[i],self.ib[i],ip] = lid
costmap[self.ij[i],self.ib[i],ip] = self.cost[i]
#
gl = part.gather(procgrid.landscatter(locmap, order='F'))
gc = part.gather(procgrid.landscatter(costmap, order='F'))
#
# Eliminate locations with lower cost
#
if part.rank == 0:
u,n = np.unique(gl[~np.isnan(gl)], return_counts=True)
# Go through all IDs which are not unique
for lid in u[np.where(n > 1)[0]] :
maxcost = np.max(gc[gl == lid])
gl[(gl == lid) & (gc < maxcost)] = np.nan
gc[(gl == lid) & (gc < maxcost)] = np.nan
#
# 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
return
#
Partition.py
View file @ b9ab8b9b
......@@ -538,18 +538,18 @@ class partition :
#
# Gather all fields partitioned in the 2D domain onto the root proc
#
def gather(self, x) :
def gather(self, x, default=np.nan) :
#
indim = x.shape
outdim = indim[:-2]+(self.njg,self.nig)
if self.rank == 0 :
xout = np.zeros(outdim, dtype=x.dtype)
if len(outdim) == 2 :
xout[:,:] = np.nan
xout[:,:] = default
elif len(outdim) == 3 :
xout[:,:,:] = np.nan
xout[:,:,:] = default
elif len(outdim) == 4 :
xout[:,:,:,:] = np.nan
xout[:,:,:,:] = default
else :
ERROR("Unforessen rank of field to be gathered")
sys.exit()
......@@ -709,6 +709,11 @@ class partition :
def domainmin(self, x) :
return min(self.comm.allgather(x))
#
# Simple access to bcast function
#
def bcast(self, x) :
return self.comm.bcast(x, root=0)
#
# Function to gather vectors of points on different procs on root
#
def gatherbypoint(self, x) :
......
RoutingPreProc.py
View file @ b9ab8b9b
......@@ -130,15 +130,7 @@ INFO("=================== Locate gauging stations and other structures =========
loc = LO.Locations(modelgrid.box_land)
hgraph.position(modelgrid, loc)
for i, id in enumerate(loc.id) :
if loc.ij[i] >= 0 :
print(rank," --> Locations to be placed : ", id)
print("Locations : lon, lat, upstream : ", loc.lon[i], loc.lat[i], loc.upstream[i])
print("Locations : ij, ib, cost :", loc.ij[i], loc.ib[i], loc.cost[i])
print("Locations in HTU :", hgraph.routing_cg[loc.ij[i],loc.ib[i],:], hgraph.routing_fetch[loc.ij[i],loc.ib[i]]/1.e+6)
print("Locations : type : ", loc.type[i])
hgraph.dumpnetcdf(OutGraphFile, gg, modelgrid, part, hydrotcst)
hgraph.dumpnetcdf(OutGraphFile, gg, modelgrid, part, hydrotcst, loc)
IF.closeinterface(comm)
tests/Mallorca/Structures.txt
View file @ b9ab8b9b
# ID Lon(decimal) lat upstream area (km^2) Nature
1 2.7918 39.7810 763.0 Reservoir
1 2.7918 39.7810 43.0 Reservoir
2 3.1861 39.5861 4.0 Run of river
tests/Mallorca/run_EuroSW.def
View file @ b9ab8b9b
......@@ -24,3 +24,9 @@ nbasmax = 9
#
GraphFile = EuroSW_test_graph.nc
#
# File containing infrastructures to be placed.
# Maximum error in the distance of the station in km^2
MaxDistErr = 25.0
# Maximum error in the upstream area in %
MaxUpstrErr = 25.0
StructuresFile = Structures.txt
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