GraphHydro.py 12.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10
#
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
#
11 12 13 14 15 16 17
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)
#
18 19 20 21 22 23 24 25 26 27 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 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
import routing_interface
import NcHelp as NH
import RPPtools as RPP
#
# Add time constants and possibly other prameters to the graph file
#
def addparameters(outnf, part, tcst, vtyp, NCFillValue) :
    if part.rank == 0 :
        outnf.createDimension('dimpara', 1)
        stream = outnf.createVariable("StreamTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        stream.title = "Time constant for the stream reservoir"
        stream.units = "day/m"
        stream[:] = tcst.stream_tcst
        fast = outnf.createVariable("FastTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        fast.title = "Time constant for the fast reservoir"
        fast.units = "day/m"
        fast[:] = tcst.fast_tcst
        slow = outnf.createVariable("SlowTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        slow.title = "Time constant for the slow reservoir"
        slow.units = "day/m"
        slow[:] = tcst.slow_tcst
        flood = outnf.createVariable("FloodTimeCst", vtyp, ('dimpara'), fill_value=NCFillValue)
        flood.title = "Time constant for the flood reservoir"
        flood.units = "day/m"
        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, \
            self.routing_floodp, self.routing_cg, self.topo_resid, 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 = \
                                    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, \
                                                                      basin_cg = hydrosuper.basin_cg, \
                                                                      basin_topoind = hydrosuper.basin_topoind, fetch_basin = hydrosuper.fetch_basin, \
                                                                      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, \
                                                                      inflow_basin = hydrosuper.inflow_basin)

        #
        #
        #
        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])

        return
    #
    #
    #
98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
    def position(self, modelgrid, locations) :
        npt, nbas = self.routing_fetch.shape
        if len(locations.id) > 0 :
            for i,id in enumerate(locations.id) :
                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)
                    locations.addhtupos(i, ij, ib, np.nanmax(cost))
    #
    #
    #
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233
    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')

        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 :
            ncvar = np.zeros((1,1,1))
        ncvar[:] = part.gather(var)
    #
    #
    #
    def dumpnetcdf(self, filename, globalgrid, procgrid, part, tcst) :
        #
        NCFillValue=1.0e20
        vtyp=np.float64
        cornerind=[0,2,4,6]
        nbcorners = len(cornerind)
        #
        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)
        else :
            outnf = None
        #
        NH.addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
        NH.addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt)
        addparameters(outnf, part, tcst, vtyp, NCFillValue)
        #
        # 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)
        self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "nbpt_glo", "Land point indices on global grid", "-", nbpt_glo[:,0], vtyp)
        #
        ################
        #
        # 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.
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, "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")
        # basin id
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z', 'y','x'), "basinid", "ID of basin", "-", self.global_basinid[:,:], vtyp, "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_orog_mean", "Mean orography", "m", self.routing_orog_mean[:,:], vtyp, "float")

        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_max", "Max orography", "m", self.routing_orog_max[:,:], vtyp, "float")

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

        # 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")
        #
        # 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")
        #
        # 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")
        # 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")
        # type of outlet
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp, "float")
        #
        # 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")
        #

        # Inflow number
        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, "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")
        #
        # 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")

        #
        # 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_lat", "Latitude of centre of gravity of HTU", "degrees north", self.routing_cg[:,:,0], vtyp, "float")
        #
        # 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")
        #
        # Close the file
        if part.rank == 0:
            outnf.close()
        #
        return