Docker-in-Docker (DinD) capabilities of public runners deactivated. More info

GraphHydro.py 21.7 KB
Newer Older
1 2 3 4
#
import sys
import numpy as np
from netCDF4 import Dataset
5
import datetime
6
#
7 8 9 10 11 12 13
import getargs
config = getargs.SetupConfig()
# Maximum error in the distance of the station in km
MaxDistErr = config.getfloat("OverAll", "MaxDistErr", fallback=25.0)
# Maximum error in the upstream area in %
MaxUpstrErr = config.getfloat("OverAll", "MaxUpstrErr", fallback=10.0)
#
14
log_master, log_world = getargs.getLogger(__name__)
15 16 17
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
#
18 19 20 21
import F90tools as F90T
#
F90T.adddirtopath("F90subroutines")
F90T.compilef90("routing_interface")
22 23 24 25
import routing_interface
import NcHelp as NH
import RPPtools as RPP
#
26 27
NCFillValue=1.0e20
#
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
# Function to compute quantiles with weights
#
def weighted_quantile(values, quantiles, sample_weight=None, 
                      values_sorted=False, old_style=False):
    """ Very close to numpy.percentile, but supports weights.
    NOTE: quantiles should be in [0, 1]!
    :param values: numpy.array with data
    :param quantiles: array-like with many quantiles needed
    :param sample_weight: array-like of the same length as `array`
    :param values_sorted: bool, if True, then will avoid sorting of
        initial array
    :param old_style: if True, will correct output to be consistent
        with numpy.percentile.
    :return: numpy.array with computed quantiles.
    """
    values = np.array(values)
    quantiles = np.array(quantiles)
    if sample_weight is None:
        sample_weight = np.ones(len(values))
    sample_weight = np.array(sample_weight)
    assert np.all(quantiles >= 0) and np.all(quantiles <= 1), \
        'quantiles should be in [0, 1]'

    if not values_sorted:
        sorter = np.argsort(values)
        values = values[sorter]
        sample_weight = sample_weight[sorter]

    weighted_quantiles = np.cumsum(sample_weight) - 0.5 * sample_weight
    if old_style:
        # To be convenient with numpy.percentile
        weighted_quantiles -= weighted_quantiles[0]
        weighted_quantiles /= weighted_quantiles[-1]
    else:
        weighted_quantiles /= np.sum(sample_weight)
    return np.interp(quantiles, weighted_quantiles, values)
#
65 66
# Add time constants and possibly other prameters to the graph file
#
67
def addparameters(outnf, version, part, tcst, vtyp, NCFillValue) :
68
    if part.rank == 0 :
69 70 71 72
        # Add global attributes to the netCDF file
        outnf.setncattr("RoutingPPVersion", version)
        outnf.setncattr("GenerationDate", datetime.datetime.utcnow().isoformat(sep=" "))
        # Add routing parameters
73 74 75
        outnf.createDimension('dimpara', 1)
        stream = outnf.createVariable("StreamTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        stream.title = "Time constant for the stream reservoir"
76
        stream.units = "s/km"
77 78 79
        stream[:] = tcst.stream_tcst
        fast = outnf.createVariable("FastTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        fast.title = "Time constant for the fast reservoir"
80
        fast.units = "s/km"
81 82 83
        fast[:] = tcst.fast_tcst
        slow = outnf.createVariable("SlowTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        slow.title = "Time constant for the slow reservoir"
84
        slow.units = "s/km"
85 86 87
        slow[:] = tcst.slow_tcst
        flood = outnf.createVariable("FloodTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        flood.title = "Time constant for the flood reservoir"
88
        flood.units = "s/km"
89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107
        flood[:] = tcst.flood_tcst
        swamp = outnf.createVariable("SwampCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        swamp.title = "Fraction of the river transport that flows to the swamps"
        swamp.units = "-"
        swamp[:] = tcst.swamp_cst
    return
#
#
#
class HydroGraph :
    def __init__(self, nbasmax, hydrosuper, part, modelgrid) :
        #
        self.nbasmax = nbasmax
        self.nbpt = hydrosuper.basin_count.shape[0]
        nwbas = hydrosuper.basin_topoind.shape[1]
        nbxmax_in = hydrosuper.inflow_grid.shape[2]
        #
        #
        self.routing_area, self.routing_orog_mean, self.routing_orog_min,self.routing_orog_max, \
108
            self.routing_floodp,self.routing_beta, self.routing_cg, self.topo_resid, self.topo_resid_stream, self.topo_rlen, \
109 110
            self.topo_rdz, self.route_nbbasin, self.route_togrid, self.route_tobasin, self.route_nbintobas, \
            self.global_basinid, self.route_outlet, self.route_type, self.origin_nbintobas, self.routing_fetch, self.floodcri = \
111 112 113 114 115 116 117
                                    routing_interface.finish_truncate(nbpt = self.nbpt, inflowmax = nbxmax_in, nbasmax = nbasmax, nwbas = nwbas, \
                                                                      num_largest = hydrosuper.num_largest, gridarea = modelgrid.area, \
                                                                      cfrac = modelgrid.contfrac, gridcenters = np.array(modelgrid.centers), \
                                                                      basin_count = hydrosuper.basin_count, \
                                                                      basin_notrun = hydrosuper.basin_notrun, basin_area = hydrosuper.basin_area, \
                                                                      basin_orog_mean = hydrosuper.basin_orog_mean, basin_orog_min = hydrosuper.basin_orog_min,\
                                                                      basin_orog_max = hydrosuper.basin_orog_max, basin_floodp = hydrosuper.basin_floodp, \
118
                                                                      basin_cg = hydrosuper.basin_cg, basin_topoind = hydrosuper.basin_topoind, \
119
                                                                      basin_topoindex_stream = hydrosuper.basin_topoindex_stream, \
120
                                                                      basin_rlen = hydrosuper.basin_rlen, basin_rdz = hydrosuper.basin_rdz, \
121
                                                                      basin_beta_fp = hydrosuper.basin_beta_fp , fetch_basin = hydrosuper.fetch_basin, \
122 123 124 125 126
                                                                      basin_id = hydrosuper.basin_id, \
                                                                      basin_coor = hydrosuper.basin_outcoor, basin_type = hydrosuper.basin_type, \
                                                                      basin_flowdir = hydrosuper.basin_flowdir, \
                                                                      outflow_grid = hydrosuper.outflow_grid, outflow_basin = hydrosuper.outflow_basin, \
                                                                      inflow_number = hydrosuper.inflow_number, inflow_grid = hydrosuper.inflow_grid, \
127 128
                                                                      inflow_basin = hydrosuper.inflow_basin, floodcri = hydrosuper.floodcri, \
                                                                      rfillval = NCFillValue, ifillval = RPP.IntFillValue)
129 130 131 132 133 134 135 136 137 138 139
        #
        self.num_largest = routing_interface.finalrivclass(part.landcorelist, self.route_togrid, self.route_tobasin, self.routing_fetch, \
                      hydrosuper.largest_rivarea)
        #
        # Inflows
        self.max_inflow = part.domainmax(np.max(hydrosuper.inflow_number))
        gingrid = part.l2glandindex( hydrosuper.inflow_grid[:,:,:self.max_inflow])
        self.route_innum, self.route_ingrid, self.route_inbasin = \
            routing_interface.finish_inflows(nbpt = self.nbpt, nwbas = nwbas, nbasmax = nbasmax, inf_max = self.max_inflow, \
                                             basin_count = hydrosuper.basin_count, inflow_number = hydrosuper.inflow_number, \
                                             inflow_grid = gingrid, inflow_basin = hydrosuper.inflow_basin[:,:,:self.max_inflow])
140
        #
141 142 143 144
        return
    #
    #
    #
145 146
    def position(self, modelgrid, locations) :
        npt, nbas = self.routing_fetch.shape
147
        if len(locations.lid) > 0 :
148 149 150
          for i,lid in enumerate(locations.lid) :
            # Some stations have a negative upstream area (because undefined)
            if locations.upstream[i] > 0:
151 152 153 154 155 156 157 158 159 160 161
                dist = RPP.distances(np.array(modelgrid.coordll)[:,0], np.array(modelgrid.coordll)[:,1], \
                                    locations.lon[i], locations.lat[i])
                # Calculation in km
                disterr = (MaxDistErr-dist/1000.)/MaxDistErr
                disterr[disterr < 0] = np.nan
                # Calculation in km^2 and result in %
                upstrerr = (MaxUpstrErr-np.abs(self.routing_fetch/1.e+6-locations.upstream[i])/locations.upstream[i]*100)/MaxUpstrErr
                upstrerr[upstrerr<0] = np.nan
                cost = np.tile(disterr,nbas).reshape((nbas,npt)).transpose()+upstrerr
                if np.count_nonzero(~np.isnan(cost)) > 0 :
                    ij, ib = np.unravel_index(np.nanargmax(cost),cost.shape)
162 163
                    userr = (self.routing_fetch[ij,ib]/1.e+6-locations.upstream[i])/locations.upstream[i]*100
                    locations.addhtupos(i, ij, ib, np.nanmax(cost), userr)
164
                else :
165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182
                    # 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)+" % ")
183
    #
184 185
    #
    #
186 187 188 189
    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)
            
190 191 192 193 194 195 196 197 198 199 200
        if orig_type == "float":
            var[np.isnan(var)] = NCFillValue
        elif orig_type == "int":
            var[np.isnan(var)] = RPP.IntFillValue
            var[var==RPP.IntFillValue] = NCFillValue

        if part.rank == 0:
            ncvar = outnf.createVariable(name, vtyp, coord, fill_value=NCFillValue)
            ncvar.title = title
            ncvar.units = units
        else :
201 202
            ncvar = np.zeros((1,)*len(coord))
        ncvar[:] = part.gather(var, default=NCFillValue)
203 204 205 206
        return
    #
    #
    #
207
    def add_tstepdistrib(self, outnf, procgrid, NCFillValue, part, data, area, vtyp, nbbins, tcst, arrayorder = 'F') :
208 209 210 211

        
        var = procgrid.landscatter(data.astype(vtyp), order=arrayorder)
        gvar = part.gather(var, default=NCFillValue)
212 213
        larea = procgrid.landscatter(area.astype(vtyp), order=arrayorder)
        garea = part.gather(larea, default=NCFillValue)
214 215 216
        
        if part.rank == 0:
            gvar[gvar >= NCFillValue] = np.nan
217
            #
POLCHER Jan's avatar
POLCHER Jan committed
218
            hist, bin_edges = np.histogram(gvar[~np.isnan(gvar)], bins=nbbins, range=(0,min(10000, np.nanmax(gvar))))
219 220 221
            whist, bin_edges = np.histogram(gvar[~np.isnan(gvar)], weights=garea[~np.isnan(gvar)], \
                                             bins=nbbins, range=(0,min(10000, np.nanmax(gvar))))
            #
222 223 224 225 226 227 228 229 230
            topobins = outnf.createVariable("topobins", vtyp, ('nbbnds','nbbins'))
            topobins.title = "Topographic index"
            topobins.units = "km"
            topobins[0,:] = bin_edges[:-1]
            topobins[1,:] = bin_edges[1:]
            topohist = outnf.createVariable("topocount", vtyp, ("nbbins",), fill_value=NCFillValue)
            topohist.title = "Counts of HTU by topoindex"
            topohist.units = "counts"
            topohist[:] = hist[:]
231 232 233 234
            wtopohist = outnf.createVariable("wtopocount", vtyp, ("nbbins",), fill_value=NCFillValue)
            wtopohist.title = "Area weighted counts of HTU by topoindex"
            wtopohist.units = "counts"
            wtopohist[:] = whist[:]
235
            #
236
            thist, tbin_edges = np.histogram(gvar[~np.isnan(gvar)]*tcst.stream_tcst, bins=nbbins, range=(0,RPP.OneDay/8))
237 238 239 240 241 242
            tbins = outnf.createVariable("tstepbins", vtyp, ('nbbnds','nbbins'))
            tbins.title = "Time bins"
            tbins.units = "s"
            tbins[0,:] = tbin_edges[:-1]
            tbins[1,:] = tbin_edges[1:]
            ncvar = outnf.createVariable("stabtstepcount", vtyp, ("nbbins",), fill_value=NCFillValue)
243 244
            ncvar.title = "Stable timestep distribution"
            ncvar.units = "count"
245
            ncvar[:] = thist[:]
246
            #
247
            maxtstep = weighted_quantile(gvar[~np.isnan(gvar)]*tcst.stream_tcst, 0.50, sample_weight=garea[~np.isnan(gvar)])
248 249 250 251
            mtstep = outnf.createVariable("MaxTimeStep", vtyp, ('dimpara'), fill_value=NCFillValue)
            mtstep.title = "The maximum time-step possible given the distribution of HTU properties"
            mtstep.units = "s"
            mtstep[:] = maxtstep
252
        return
253 254 255
    #
    #
    #
256
    def dumpnetcdf(self, filename, version, globalgrid, procgrid, part, tcst, locations) :
257 258 259 260
        #
        vtyp=np.float64
        cornerind=[0,2,4,6]
        nbcorners = len(cornerind)
261
        nbbins = 2000
262
        #
263
        nlmax, nblocated = locations.maxlocations(self.nbasmax, part)
264
        #
265 266 267 268 269 270 271 272 273
        if part.rank == 0 :
            outnf=Dataset(filename, 'w', format='NETCDF4_CLASSIC')
            # Dimensions
            outnf.createDimension('x', globalgrid.ni)
            outnf.createDimension('y', globalgrid.nj)
            outnf.createDimension('land', len(procgrid.area))
            outnf.createDimension('z', self.nbasmax)
            outnf.createDimension('bnd', nbcorners)
            outnf.createDimension('inflow', self.max_inflow)
274
            outnf.createDimension('stnperhtu', nlmax)
275
            outnf.createDimension('nbbins', nbbins)
276
            outnf.createDimension('nbbnds',2)
277
            outnf.createDimension('locations', nblocated)
278 279 280 281 282
        else :
            outnf = None
        #
        NH.addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
        NH.addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt)
283
        addparameters(outnf, version, part, tcst, vtyp, NCFillValue)
284 285 286 287 288 289
        #
        # land grid index -> to facilitate the analyses of the routing
        #
        nbpt_loc = np.zeros((self.nbpt,1)).astype(np.int32)
        nbpt_loc[:,0] = np.arange(1, self.nbpt+1)
        nbpt_glo = part.l2glandindex(nbpt_loc)
290 291
        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), \
                          "nbpt_glo", "Land point indices on global grid", "-", nbpt_glo[:,0], vtyp)
292 293 294 295 296 297 298 299
        #
        ################
        #
        # Conversion of grid indices for route_togrid
        grgrid = part.l2glandindex(self.route_togrid[:,:])
        #
        #
        # The field route_togrid is with indices on the local grid. That needs to be converted to the global grid.
300 301
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, orig_type="int")
302
        # route to basin
303 304
        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")
305
        # basin id
306 307
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "basinid", "ID of basin", "-", self.global_basinid[:,:], vtyp, orig_type="int")
308 309
        #
        # basin area
310 311
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "basin_area", "area of basin", "m^2", self.routing_area[:,:], vtyp)
312

313 314
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "basin_orog_mean", "Mean orography", "m", self.routing_orog_mean[:,:], vtyp)
315

316 317
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "basin_orog_min", "Min orography", "m", self.routing_orog_min[:,:], vtyp)
318

319 320
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "basin_orog_max", "Max orography", "m", self.routing_orog_max[:,:], vtyp)
321

322 323
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "basin_floodp", "Fraction of floodplains", "-", self.routing_floodp[:,:], vtyp)
324

325 326
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "floodcri", "Height at which all basin is flooded", "mm", self.floodcri[:,:], vtyp)
327

328 329
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "basin_beta_fp", "Shape of the floodplins HTUs (concave / convex)", "-", self.routing_beta[:,:], vtyp)
330

331
        # route number into basin
332 333
        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), \
                          "routenbintobas", "Number of basin into current one", "-", self.route_nbintobas[:], vtyp, orig_type="int")
334 335
        #
        # original number into basin
336 337
        self.add_variable(outnf, procgrid, NCFillValue, part, ( 'y','x'), \
                          "originnbintobas", "Number of sub-grid basin into current one before truncation", "-", \
338
                          self.origin_nbintobas[:], vtyp, orig_type="int")
339 340
        #
        # latitude of outlet
341 342
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp)
343
        # longitude of outlet
344 345
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "outletlon", "Longitude of Outlet", "degrees east", self.route_outlet[:,:,1], vtyp)
346
        # type of outlet
347 348
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp)
349 350
        #
        # topographic index
351 352
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "topoindex", "Topographic index of the retention time", "km", self.topo_resid[:,:], vtyp)
353 354
        self.add_tstepdistrib(outnf, procgrid, NCFillValue, part, self.topo_resid[:,:], self.routing_area[:,:], \
                              vtyp, nbbins, tcst)
355 356 357 358 359 360 361 362
        #
        # topoindex stream
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "topoindex_stream", "Topographic index of the retention time for stream res.", "km", self.topo_resid_stream[:,:], vtyp)
        # Add another distribution ? Or replace the other ?
        #self.add_tstepdistrib(outnf, procgrid, NCFillValue, part, self.topo_resid_stream[:,:], self.routing_area[:,:], \
        #                      vtyp, nbbins, tcst)

363
        #
364 365 366 367 368
        # Geometric properties of HTU
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "riverlen", "Mean river length within HTU", "m", self.topo_rlen[:,:], vtyp)
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "riverdz", "Mean river elevation change within HTU", "m", self.topo_rdz[:,:], vtyp)
369
        # Inflow number
370 371
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, orig_type="int")
372 373 374
        #
        # Inflow grid
        #gingrid = part.l2glandindex(self.inflow_grid[:,:,:inflow_size])
375 376
        self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow','z','y','x'), \
                          "route_ingrid", "Grid from which the water flows", "-", \
377
                          self.route_ingrid[:,:,:], vtyp, orig_type="int")
378 379
        #
        # Inflow basin
380 381
        self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow','z','y','x'), \
                          "route_inbasin", "Basin from which the water flows", "-", \
382
                          self.route_inbasin[:,:,:], vtyp, orig_type="int")
383 384 385 386 387

        #
        # Save centre of gravity of HTU
        #
        # Check if it works
388 389
        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)
390

391 392
        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)
393 394 395
        #
        # Save the fetch of each basin
        #
396 397
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp)
398 399 400
        #
        # Build and save a map of all the structures located on the graph
        #
401
        locations.addtonetcdf(self.nbpt, self.nbasmax, outnf, procgrid, part, ('y','x'), ('locations',), vtyp)
402 403
        #
        # Close the file
404
        #
405 406 407 408
        if part.rank == 0:
            outnf.close()
        #
        return