GraphHydro.py 18.3 KB
Newer Older
1 2 3 4 5 6 7 8 9
#
import os
import sys
from inspect import currentframe, getframeinfo
localdir=os.path.dirname(getframeinfo(currentframe()).filename)
sys.path.append(localdir+'/F90subroutines')
#
import numpy as np
from netCDF4 import Dataset
10
import datetime
11
#
12 13 14 15 16 17 18
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)
#
19
log_master, log_world = getargs.getLogger(__name__)
20 21 22
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
#
23 24 25 26
import routing_interface
import NcHelp as NH
import RPPtools as RPP
#
27 28
NCFillValue=1.0e20
#
29 30
# Add time constants and possibly other prameters to the graph file
#
31
def addparameters(outnf, version, part, tcst, vtyp, NCFillValue) :
32
    if part.rank == 0 :
33 34 35 36
        # Add global attributes to the netCDF file
        outnf.setncattr("RoutingPPVersion", version)
        outnf.setncattr("GenerationDate", datetime.datetime.utcnow().isoformat(sep=" "))
        # Add routing parameters
37 38 39
        outnf.createDimension('dimpara', 1)
        stream = outnf.createVariable("StreamTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        stream.title = "Time constant for the stream reservoir"
40
        stream.units = "s/km"
41 42 43
        stream[:] = tcst.stream_tcst
        fast = outnf.createVariable("FastTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        fast.title = "Time constant for the fast reservoir"
44
        fast.units = "s/km"
45 46 47
        fast[:] = tcst.fast_tcst
        slow = outnf.createVariable("SlowTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        slow.title = "Time constant for the slow reservoir"
48
        slow.units = "s/km"
49 50 51
        slow[:] = tcst.slow_tcst
        flood = outnf.createVariable("FloodTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        flood.title = "Time constant for the flood reservoir"
52
        flood.units = "s/km"
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
        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, \
72
            self.routing_floodp,self.routing_beta, self.routing_cg, self.topo_resid, self.route_nbbasin,\
73
            self.route_togrid, self.route_tobasin, self.route_nbintobas, self.global_basinid, \
74
            self.route_outlet, self.route_type, self.origin_nbintobas, self.routing_fetch, self.floodcri = \
75 76 77 78 79 80 81
                                    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, \
82 83
                                                                      basin_cg = hydrosuper.basin_cg, basin_topoind = hydrosuper.basin_topoind, \
                                                                      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 167 168 169 170 171 172 173 174
        return
    #
    #
    #
    def add_tstepdistrib(self, outnf, procgrid, NCFillValue, part, data, vtyp, nbbins, tcst, arrayorder = 'F') :

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

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

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

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

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

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

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

282
        # route number into basin
283 284
        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), \
                          "routenbintobas", "Number of basin into current one", "-", self.route_nbintobas[:], vtyp, orig_type="int")
285 286
        #
        # original number into basin
287 288
        self.add_variable(outnf, procgrid, NCFillValue, part, ( 'y','x'), \
                          "originnbintobas", "Number of sub-grid basin into current one before truncation", "-", \
289
                          self.origin_nbintobas[:], vtyp, orig_type="int")
290 291
        #
        # latitude of outlet
292 293
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp)
294
        # longitude of outlet
295 296
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "outletlon", "Longitude of Outlet", "degrees east", self.route_outlet[:,:,1], vtyp)
297
        # type of outlet
298 299
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp)
300 301
        #
        # topographic index
302 303 304
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "topoindex", "Topographic index of the retention time", "km", self.topo_resid[:,:], vtyp)
        self.add_tstepdistrib(outnf, procgrid, NCFillValue, part, self.topo_resid[:,:], vtyp, nbbins, tcst)
305 306 307
        #

        # Inflow number
308 309
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, orig_type="int")
310 311 312
        #
        # Inflow grid
        #gingrid = part.l2glandindex(self.inflow_grid[:,:,:inflow_size])
313 314
        self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow','z','y','x'), \
                          "route_ingrid", "Grid from which the water flows", "-", \
315
                          self.route_ingrid[:,:,:], vtyp, orig_type="int")
316 317
        #
        # Inflow basin
318 319
        self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow','z','y','x'), \
                          "route_inbasin", "Basin from which the water flows", "-", \
320
                          self.route_inbasin[:,:,:], vtyp, orig_type="int")
321 322 323 324 325

        #
        # Save centre of gravity of HTU
        #
        # Check if it works
326 327
        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)
328

329 330
        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)
331 332 333
        #
        # Save the fetch of each basin
        #
334 335
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
                          "fetch", "Fetch contributing to each HTU", "m^2", self.routing_fetch[:,:], vtyp)
336 337 338
        #
        # Build and save a map of all the structures located on the graph
        #
339
        locations.addtonetcdf(self.nbpt, self.nbasmax, outnf, procgrid, part, ('y','x'), ('locations',), vtyp)
340 341
        #
        # Close the file
342
        #
343 344 345 346
        if part.rank == 0:
            outnf.close()
        #
        return