Source code for cwatm.hydrological_modules.routing_reservoirs.routing_sub

#-------------------------------------------------------------------------------
# Name:        routing_sub
# Purpose:     subroutines for routing kinematic waves
#
# Author:      burekpe
#
# Created:     17/01/2017
# Copyright:   (c) burek 2017
#-------------------------------------------------------------------------------


import numpy as np
import math
from cwatm.management_modules.data_handling import *

"""
ROUTING subroutines
partly using C++ for speeding up
"""


[docs]def Compress (map,mask): """ compressing map from 2D to 1D without missing values :param map: input map :param mask: mask map :return: compressed map """ maskmap = np.ma.masked_array(map,mask) compmap = np.ma.compressed(maskmap) return compmap
[docs]def decompress1(map): """ Decompressing map from 1D to 2D with missing values :param map: compressed map :return: decompressed 2D map """ dmap = maskinfo['maskall'].copy() dmap[~maskinfo['maskflat']] = map[:] dmap = dmap.reshape(maskinfo['shape']) return dmap
[docs]def postorder(dirUp,catchment,node,catch,dirDown): """ Routine to run a postorder tree traversal :param dirUp: :param catchment: :param node: :param catch: :param dirDown: :return: dirDown and catchment """ if dirUp[node] != []: postorder(dirUp,catchment,dirUp[node][0],catch,dirDown) catchment[dirUp[node][0]]=catch dirDown.append(dirUp[node][0]) if len(dirUp[node])>1: postorder(dirUp,catchment,dirUp[node][1],catch,dirDown) catchment[dirUp[node][1]]=catch dirDown.append(dirUp[node][1]) if len(dirUp[node])>2: postorder(dirUp,catchment,dirUp[node][2],catch,dirDown) catchment[dirUp[node][2]]=catch dirDown.append(dirUp[node][2]) if len(dirUp[node])>3: postorder(dirUp,catchment,dirUp[node][3],catch,dirDown) catchment[dirUp[node][3]]=catch dirDown.append(dirUp[node][3]) if len(dirUp[node])>4: postorder(dirUp,catchment,dirUp[node][4],catch,dirDown) catchment[dirUp[node][4]]=catch dirDown.append(dirUp[node][4]) if len(dirUp[node])>5: postorder(dirUp,catchment,dirUp[node][5],catch,dirDown) catchment[dirUp[node][5]]=catch dirDown.append(dirUp[node][5]) if len(dirUp[node])>6: postorder(dirUp,catchment,dirUp[node][6],catch,dirDown) catchment[dirUp[node][6]]=catch dirDown.append(dirUp[node][6]) if len(dirUp[node])>7: postorder(dirUp,catchment,dirUp[node][7],catch,dirDown) catchment[dirUp[node][7]]=catch dirDown.append(dirUp[node][7])
[docs]def dirUpstream(dirshort): """ runs the network tree upstream from outlet to source :param dirshort: :return: direction upstream """ # -- up direction dirUp = list([] for i in range(maskinfo['mapC'][0])) for i in range(dirshort.shape[0]): value=dirshort[i] if value > -1: dirUp[value].append(i) dirupLen=[0] dirupID=[] j = 0 for i in range(dirshort.shape[0]): j += len(dirUp[i]) dirupLen.append(j) for k in range(len(dirUp[i])): dirupID.append(dirUp[i][k]) return dirUp,np.array(dirupLen).astype(np.int64),np.array(dirupID).astype(np.int64)
[docs]def dirDownstream(dirUp,lddcomp,dirDown): """ runs the river network tree downstream - from source to outlet :param dirUp: :param lddcomp: :param dirDown: :return: direction downstream """ catchment = np.array(np.zeros(maskinfo['mapC'][0]), dtype=np.int64) j=0 for pit in range(maskinfo['mapC'][0]): if lddcomp[pit]==5: j+=1 postorder(dirUp,catchment,pit,j,dirDown) dirDown.append(pit) catchment[pit]=j return np.array(dirDown).astype(np.int64),np.array(catchment).astype(np.int64)
[docs]def upstreamArea(dirDown,dirshort,area): """ calculates upstream area :param dirDown: array which point from each cell to the next downstream cell :param dirshort: :param area: area in m2 for a single gridcell :return: upstream area """ ups = area.copy() lib2.ups(dirDown,dirshort, ups,len(dirDown)) return ups
[docs]def upstream1(downstruct, weight): """ Calculates 1 cell upstream :param downstruct: :param weight: :return: upstream 1cell """ return np.bincount(downstruct, weights=weight)[:-1]
[docs]def downstream1(dirUp,weight): """ calculated 1 cell downstream :param dirUp: :param weight: :return: dowmnstream 1 cell """ downstream = weight.copy() k = 0 for i in dirUp: for j in i: downstream[j]=weight[k] k += 1 return downstream
[docs]def catchment1(dirUp, points): """ calculates all cells which belongs to a catchment from point onward :param dirUp: :param points: :return: subcatchment """ subcatch = np.array(np.zeros(maskinfo['mapC'][0]), dtype=np.int64) # if subcatchment = true -> calculation of subcatchment: every point is calculated # if false : calculation of catchment: only point calculated which are not inside a bigger catchment from another point for cell in range(maskinfo['mapC'][0]): j = points[cell] if (j > 0) and (subcatch[cell] < 1): dirDown =[] postorder(dirUp, subcatch, cell, j,dirDown) subcatch[cell] = j return subcatch
[docs]def subcatchment1(dirUp, points,ups): """ calculates subcatchments of points :param dirUp: :param points: :param ups: :return: subcatchment """ subcatch = np.array(np.zeros(maskinfo['mapC'][0]), dtype=np.int64) # if subcatchment = true -> calculation of subcatchment: every point is calculated # if false : calculation of catchment: only point calculated which are not inside a bigger catchment from another point subs = {} #sort waterbodies of reverse upstream area for cell in range(maskinfo['mapC'][0]): if points[cell] > 0: subs[points[cell]] = [cell, ups[cell]] subsort = sorted(list(subs.items()), key=lambda x: x[1][1], reverse=True) #for cell in xrange(maskinfo['mapC'][0]): for sub in subsort: j = sub[0] cell = sub[1][0] dirDown=[] postorder(dirUp, subcatch, cell, j,dirDown) subcatch[cell] = j return subcatch
# ----------------------------------------------------------------
[docs]def defLdd2(ldd): """ defines river network :param ldd: river network :return: ldd variables """ # decompressing ldd from 1D -> 2D dmap = maskinfo['maskall'].copy() dmap[~maskinfo['maskflat']] = ldd[:] ldd2D = dmap.reshape(maskinfo['shape']).astype(np.int64) ldd2D[ldd2D.mask] = 0 # every cell gets an order starting from 0 ... lddshortOrder =np.arange(maskinfo['mapC'][0]) # decompress this map to 2D lddOrder = decompress1(lddshortOrder) lddOrder[maskinfo['mask']]=-1 lddOrder = np.array(lddOrder.data, dtype=np.int64) lddCompress, dirshort = lddrepair(ldd2D, lddOrder) dirUp, dirupLen, dirupID = dirUpstream(dirshort) # for upstream calculation inAr = np.arange(maskinfo['mapC'][0], dtype=np.int64) # each upstream pixel gets the id of the downstream pixel downstruct = downstream1(dirUp, inAr).astype(np.int64) # all pits gets a high number downstruct[lddCompress == 5] = maskinfo['mapC'][0] # self.var.dirDown: direction downstream - from each cell the pointer to a downstream cell (can only be 1) # self.var.catchment: each catchment with a pit gets a own ID dirDown = [] dirDown, catchment = dirDownstream(dirUp, lddCompress,dirDown) lendirDown = len(dirDown) return lddCompress, dirshort, dirUp, dirupLen, dirupID, downstruct, catchment, dirDown, lendirDown
[docs]def lddshort(lddnp,lddOrder): """ return short for calculating a catchment from a river network :param lddnp: rivernetwork as 1D array :param lddOrder: :return: short ldd """ yi=lddnp.shape[0] xi=lddnp.shape[1] # direction downstream naming the order of the cell dir = np.array(np.empty(maskinfo['shape']), dtype=np.int64) dir.fill(-1) #lib2.repairLdd1(lddnp, yi,xi) lddcomp = compressArray(lddnp).astype(np.int64) lib2.dirID(lddOrder, lddnp, dir,yi,xi) dirshort = compressArray(dir).astype(np.int64) return dirshort
[docs]def lddrepair(lddnp,lddOrder): """ repairs a river network * eliminate unsound parts * add pits at points with no connections :param lddnp: rivernetwork as 1D array :param lddOrder: :return: repaired ldd """ yi=lddnp.shape[0] xi=lddnp.shape[1] #print "========= repair 1" # direction downstream naming the order of the cell dir = np.array(np.empty(maskinfo['shape']), dtype=np.int64) dir.fill(-1) lib2.repairLdd1(lddnp, yi,xi) lddcomp = compressArray(lddnp).astype(np.int64) lib2.dirID(lddOrder, lddnp, dir,yi,xi) dirshort = compressArray(dir).astype(np.int64) check = np.array(np.zeros(maskinfo['mapC'][0]),dtype=np.int64) lib2.repairLdd2(lddcomp, dirshort, check,maskinfo['mapC'][0] ) ii=2 """ for i in range(maskinfo['mapC'][0]): path=[] j=i while 1: if j in path: lddcomp[path[-1]] = 5 dirshort[path[-1]] = -1 break if (lddcomp[j]==5) or (check[j]==1): break path.append(j) j = dirshort[j] check[path]=1 """ return lddcomp, dirshort