GraphHydro.py 19.4 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
# Add time constants and possibly other prameters to the graph file
#
30
def addparameters(outnf, version, part, tcst, vtyp, NCFillValue) :
31
    if part.rank == 0 :
32 33 34 35
        # Add global attributes to the netCDF file
        outnf.setncattr("RoutingPPVersion", version)
        outnf.setncattr("GenerationDate", datetime.datetime.utcnow().isoformat(sep=" "))
        # Add routing parameters
36 37 38
        outnf.createDimension('dimpara', 1)
        stream = outnf.createVariable("StreamTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        stream.title = "Time constant for the stream reservoir"
39
        stream.units = "s/km"
40 41 42
        stream[:] = tcst.stream_tcst
        fast = outnf.createVariable("FastTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        fast.title = "Time constant for the fast reservoir"
43
        fast.units = "s/km"
44 45 46
        fast[:] = tcst.fast_tcst
        slow = outnf.createVariable("SlowTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        slow.title = "Time constant for the slow reservoir"
47
        slow.units = "s/km"
48 49 50
        slow[:] = tcst.slow_tcst
        flood = outnf.createVariable("FloodTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        flood.title = "Time constant for the flood reservoir"
51
        flood.units = "s/km"
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
        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, \
71 72 73
            self.routing_floodp,self.routing_beta, self.routing_cg, self.topo_resid, self.topo_rlen, \
            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 = \
74 75 76 77 78 79 80
                                    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, \
81
                                                                      basin_cg = hydrosuper.basin_cg, basin_topoind = hydrosuper.basin_topoind, \
82
                                                                      basin_rlen = hydrosuper.basin_rlen, basin_rdz = hydrosuper.basin_rdz, \
83
                                                                      basin_beta_fp = hydrosuper.basin_beta_fp , fetch_basin = hydrosuper.fetch_basin, \
84 85 86 87 88
                                                                      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, \
89 90
                                                                      inflow_basin = hydrosuper.inflow_basin, floodcri = hydrosuper.floodcri, \
                                                                      rfillval = NCFillValue, ifillval = RPP.IntFillValue)
91 92 93 94 95 96 97 98 99 100 101
        #
        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])
102
        #
103 104 105 106
        return
    #
    #
    #
107 108
    def position(self, modelgrid, locations) :
        npt, nbas = self.routing_fetch.shape
109 110
        if len(locations.lid) > 0 :
            for i,lid in enumerate(locations.lid) :
111 112 113 114 115 116 117 118 119 120 121
                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)
122 123
                    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)
124
                else :
125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142
                    # 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)+" % ")
143
    #
144 145
    #
    #
146 147 148 149
    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)
            
150 151 152 153 154 155 156 157 158 159 160
        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 :
161 162
            ncvar = np.zeros((1,)*len(coord))
        ncvar[:] = part.gather(var, default=NCFillValue)
163 164 165 166
        return
    #
    #
    #
167
    def add_tstepdistrib(self, outnf, procgrid, NCFillValue, part, data, area, vtyp, nbbins, tcst, arrayorder = 'F') :
168 169 170 171

        
        var = procgrid.landscatter(data.astype(vtyp), order=arrayorder)
        gvar = part.gather(var, default=NCFillValue)
172 173
        larea = procgrid.landscatter(area.astype(vtyp), order=arrayorder)
        garea = part.gather(larea, default=NCFillValue)
174 175 176
        
        if part.rank == 0:
            gvar[gvar >= NCFillValue] = np.nan
177
            #
POLCHER Jan's avatar
POLCHER Jan committed
178
            hist, bin_edges = np.histogram(gvar[~np.isnan(gvar)], bins=nbbins, range=(0,min(10000, np.nanmax(gvar))))
179 180 181
            whist, bin_edges = np.histogram(gvar[~np.isnan(gvar)], weights=garea[~np.isnan(gvar)], \
                                             bins=nbbins, range=(0,min(10000, np.nanmax(gvar))))
            #
182 183 184 185 186 187 188 189 190
            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[:]
191 192 193 194
            wtopohist = outnf.createVariable("wtopocount", vtyp, ("nbbins",), fill_value=NCFillValue)
            wtopohist.title = "Area weighted counts of HTU by topoindex"
            wtopohist.units = "counts"
            wtopohist[:] = whist[:]
195
            #
196
            thist, tbin_edges = np.histogram(gvar[~np.isnan(gvar)]*tcst.stream_tcst, bins=nbbins, range=(0,RPP.OneDay/8))
197 198 199 200 201 202
            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)
203 204
            ncvar.title = "Stable timestep distribution"
            ncvar.units = "count"
205
            ncvar[:] = thist[:]
206 207 208 209 210 211
            #
            maxtstep = np.quantile(gvar[~np.isnan(gvar)], 0.1)
            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
212
        return
213 214 215
    #
    #
    #
216
    def dumpnetcdf(self, filename, version, globalgrid, procgrid, part, tcst, locations) :
217 218 219 220
        #
        vtyp=np.float64
        cornerind=[0,2,4,6]
        nbcorners = len(cornerind)
221
        nbbins = 2000
222
        #
223
        nlmax, nblocated = locations.maxlocations(self.nbasmax, part)
224
        #
225 226 227 228 229 230 231 232 233
        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)
234
            outnf.createDimension('stnperhtu', nlmax)
235
            outnf.createDimension('nbbins', nbbins)
236
            outnf.createDimension('nbbnds',2)
237
            outnf.createDimension('locations', nblocated)
238 239 240 241 242
        else :
            outnf = None
        #
        NH.addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
        NH.addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt)
243
        addparameters(outnf, version, part, tcst, vtyp, NCFillValue)
244 245 246 247 248 249
        #
        # 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)
250 251
        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), \
                          "nbpt_glo", "Land point indices on global grid", "-", nbpt_glo[:,0], vtyp)
252 253 254 255 256 257 258 259
        #
        ################
        #
        # 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.
260 261
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, orig_type="int")
262
        # route to basin
263 264
        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")
265
        # basin id
266 267
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "basinid", "ID of basin", "-", self.global_basinid[:,:], vtyp, orig_type="int")
268 269
        #
        # basin area
270 271
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "basin_area", "area of basin", "m^2", self.routing_area[:,:], vtyp)
272

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

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

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

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

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

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

291
        # route number into basin
292 293
        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), \
                          "routenbintobas", "Number of basin into current one", "-", self.route_nbintobas[:], vtyp, orig_type="int")
294 295
        #
        # original number into basin
296 297
        self.add_variable(outnf, procgrid, NCFillValue, part, ( 'y','x'), \
                          "originnbintobas", "Number of sub-grid basin into current one before truncation", "-", \
298
                          self.origin_nbintobas[:], vtyp, orig_type="int")
299 300
        #
        # latitude of outlet
301 302
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp)
303
        # longitude of outlet
304 305
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "outletlon", "Longitude of Outlet", "degrees east", self.route_outlet[:,:,1], vtyp)
306
        # type of outlet
307 308
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp)
309 310
        #
        # topographic index
311 312
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "topoindex", "Topographic index of the retention time", "km", self.topo_resid[:,:], vtyp)
313 314
        self.add_tstepdistrib(outnf, procgrid, NCFillValue, part, self.topo_resid[:,:], self.routing_area[:,:], \
                              vtyp, nbbins, tcst)
315
        #
316 317 318 319 320
        # 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)
321
        # Inflow number
322 323
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, orig_type="int")
324 325 326
        #
        # Inflow grid
        #gingrid = part.l2glandindex(self.inflow_grid[:,:,:inflow_size])
327 328
        self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow','z','y','x'), \
                          "route_ingrid", "Grid from which the water flows", "-", \
329
                          self.route_ingrid[:,:,:], vtyp, orig_type="int")
330 331
        #
        # Inflow basin
332 333
        self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow','z','y','x'), \
                          "route_inbasin", "Basin from which the water flows", "-", \
334
                          self.route_inbasin[:,:,:], vtyp, orig_type="int")
335 336 337 338 339

        #
        # Save centre of gravity of HTU
        #
        # Check if it works
340 341
        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)
342

343 344
        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)
345 346 347
        #
        # Save the fetch of each basin
        #
348 349
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp)
350 351 352
        #
        # Build and save a map of all the structures located on the graph
        #
353
        locations.addtonetcdf(self.nbpt, self.nbasmax, outnf, procgrid, part, ('y','x'), ('locations',), vtyp)
354 355
        #
        # Close the file
356
        #
357 358 359 360
        if part.rank == 0:
            outnf.close()
        #
        return