From a776322db64636ecc33267940a73e581dfdfe0a8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Falk=20He=C3=9Fe?= Date: Thu, 1 Jun 2017 13:45:36 +0000 Subject: [PATCH] trunk: post-proc: added travel time distribution and BUGFIX mhm class --- post-proc/mhm.py | 106 ++++---- post-proc/sas/__init__.py | 66 +++++ post-proc/sas/aux_fun.py | 129 ++++++++++ post-proc/sas/get_U_num.py | 12 + post-proc/sas/get_p.py | 139 +++++++++++ post-proc/sas/get_theta.py | 72 ++++++ post-proc/sas/get_validity_range.py | 67 ++++++ post-proc/sas/plot_fun.py | 123 ++++++++++ post-proc/sas/plot_study_domain.py | 339 ++++++++++++++++++++++++++ post-proc/sas/sas_base.py | 359 ++++++++++++++++++++++++++++ post-proc/test_mhm_sas_class.py | 46 ++++ 11 files changed, 1416 insertions(+), 42 deletions(-) create mode 100644 post-proc/sas/__init__.py create mode 100644 post-proc/sas/aux_fun.py create mode 100644 post-proc/sas/get_U_num.py create mode 100644 post-proc/sas/get_p.py create mode 100644 post-proc/sas/get_theta.py create mode 100644 post-proc/sas/get_validity_range.py create mode 100644 post-proc/sas/plot_fun.py create mode 100644 post-proc/sas/plot_study_domain.py create mode 100644 post-proc/sas/sas_base.py create mode 100644 post-proc/test_mhm_sas_class.py diff --git a/post-proc/mhm.py b/post-proc/mhm.py index f8e85c81..18887b01 100644 --- a/post-proc/mhm.py +++ b/post-proc/mhm.py @@ -62,13 +62,13 @@ class MHM(object): ##-- init function ------------------------------------------------------------ - def __init__(self, model_path, file_name = 'mhm.nml'): + def __init__(self, model_path, file_name = 'mhm.nml', + rel_path_name = False): # to do: - option for restart file # - both implicite and explicite # - adapt for use with several basins # - onyl simple upscaling - # - only absolute paths # - documentation # - doctests @@ -103,14 +103,7 @@ def __init__(self, model_path, file_name = 'mhm.nml'): f_id.close() self.nBasins = int(self.nml['nBasins']) - - self.cellsize['L0'] = 100 - self.cellsize['L1'] = int(self.nml['resolution_Hydrology(1)']) - self.cellsize['L2'] = 4000 - - self.cellsize_ratio['L0'] = 1 - self.cellsize_ratio['L1'] = float(self.nml['resolution_Hydrology(1)'])/self.cellsize['L0'] - self.cellsize_ratio['L2'] = self.cellsize['L2']/self.cellsize['L0'] + if self.nml['timeStep_sm_input'] == '-1': self.t_stepsize = 'daily' @@ -118,10 +111,29 @@ def __init__(self, model_path, file_name = 'mhm.nml'): self.t_stepsize = 'monthly' elif self.nml['timeStep_sm_input'] == '-3': self.t_stepsize = 'yearly' - - # reading data from the 'dem.asc' file for domain information - + + if rel_path_name: + for k,v in self.nml.items(): + if ('dir_' in k) or ('file_' in k): + self.nml[k] = model_path + self.nml[k] + + # reading data from the 'dem.asc' file for domain information + f_id = open(self.nml['dir_Morpho(1)'] + 'dem.asc') + + for i, line in enumerate(f_id): + if i == 4: + self.cellsize['L0'] = int(line[9:]) + self.cellsize['L1'] = int(self.nml['resolution_Hydrology(1)']) + self.cellsize['L2'] = 4000 + f_id.close() + + self.cellsize_ratio['L0'] = 1 + self.cellsize_ratio['L1'] = float(self.nml['resolution_Hydrology(1)'])/self.cellsize['L0'] + self.cellsize_ratio['L2'] = self.cellsize['L2']/self.cellsize['L0'] + + f_id = open(self.nml['dir_Morpho(1)'] + 'dem.asc') + for line_i in range(6): line = f_id.next().strip() line = line.split() @@ -152,18 +164,18 @@ def __init__(self, model_path, file_name = 'mhm.nml'): self.mask['L0'][line_i - 6][col_i] = 1 for row_i in range(0, self.nrows['L1']): for col_i in range(0, self.ncols['L1']): - row_a = row_i*self.cellsize_ratio['L1'] - row_e = (row_i + 1)*self.cellsize_ratio['L1'] - col_a = col_i*self.cellsize_ratio['L1'] - col_e = (col_i + 1)*self.cellsize_ratio['L1'] + row_a = int(row_i*self.cellsize_ratio['L1']) + row_e = int((row_i + 1)*self.cellsize_ratio['L1']) + col_a = int(col_i*self.cellsize_ratio['L1']) + col_e = int((col_i + 1)*self.cellsize_ratio['L1']) if np.mean(self.mask['L0'][row_a:row_e,col_a:col_e]) == 1: self.mask['L1'][row_i][col_i] = 1 for row_i in range(0, self.nrows['L2']): for col_i in range(0, self.ncols['L2']): - row_a = row_i*self.cellsize_ratio['L2'] - row_e = (row_i + 1)*self.cellsize_ratio['L2'] - col_a = col_i*self.cellsize_ratio['L2'] - col_e = (col_i + 1)*self.cellsize_ratio['L2'] + row_a = int(row_i*self.cellsize_ratio['L2']) + row_e = int((row_i + 1)*self.cellsize_ratio['L2']) + col_a = int(col_i*self.cellsize_ratio['L2']) + col_e = int((col_i + 1)*self.cellsize_ratio['L2']) if np.mean(self.mask['L0'][row_a:row_e,col_a:col_e]) == 1: self.mask['L2'][row_i][col_i] = 1 f_id.close() @@ -174,13 +186,17 @@ def import_data(self, data_type, *kwargs ): if data_type == 'states_and_fluxes': self.import_states_and_fluxes( *kwargs ) elif data_type == 'lat_lon_L0': - data_path = self.nml['dir_LatLon(1)'] - self.lon_L0 = self.import_lat_lon(data_path, 'lon_l0') - self.lat_L0 = self.import_lat_lon(data_path, 'lat_l0') + data_path = str(kwargs[0]) + self.lon_L0 = self.import_lat_lon(data_path, 'lon') + self.lat_L0 = self.import_lat_lon(data_path, 'lat') elif data_type == 'lat_lon': - data_path = self.nml['dir_LatLon(1)'] + data_path = self.nml['file_LatLon(1)'] self.lon_L1 = self.import_lat_lon(data_path, 'lon') self.lat_L1 = self.import_lat_lon(data_path, 'lat') + elif data_type == 'lat_lon_L2': + data_path = str(kwargs[0]) + self.lon_L2 = self.import_lat_lon(data_path, 'lon') + self.lat_L2 = self.import_lat_lon(data_path, 'lat') elif data_type == 'landcover': data_path = self.nml['dir_LCover(1)'] + str(kwargs[0]) setattr(self, data_type, self.import_L0_data(data_path)) @@ -213,6 +229,7 @@ def import_precipitation(self, f_path): def import_states_and_fluxes(self, *kwargs): f_path = self.nml['dir_Out(1)'] + 'mHM_Fluxes_States.nc' +# print(f_path) if kwargs[0] == 'all': var_list = [str(i) for i in readnetcdf(f_path, variables=True)] @@ -224,8 +241,8 @@ def import_states_and_fluxes(self, *kwargs): 'recharge', 'aET_L01', 'aET_L02', 'aET_L03', 'preEffect'] else: - var_list = kwargs - + var_list = kwargs[0] +# print(var_list) for var_i in range(0, len(var_list)): var = var_list[var_i] setattr(self, var, readnetcdf(f_path, var=var)) @@ -251,14 +268,14 @@ def import_L0_data(self, f_path): ##-- importing restart file --------------------------------------------------- -# def import_restart_file(self, *kwargs): -# -# f_path = self.nml['dir_Out(1)'] + 'mHM_restart_001.nc' -# var_list = [str(i) for i in readnetcdf(f_path, variables=True)] + def import_restart_file(self, *kwargs): + + f_path = self.nml['dir_Out(1)'] + 'mHM_restart_001.nc' + var_list = [str(i) for i in readnetcdf(f_path, variables=True)] # print(var_list) -# for var_i in range(0, len(var_list)): -# var = var_list[var_i] -# setattr(self, var, readnetcdf(f_path, var=var)) + for var_i in range(0, len(var_list)): + var = var_list[var_i] + setattr(self, var, readnetcdf(f_path, var=var)) ##-- upscaling and downscaling functions -------------------------------------- @@ -277,11 +294,11 @@ def upscale_L1_data(self, L0_array, flag): for col_i in range(0, self.ncols['L1']): if self.mask['L1'][row_i, col_i]: continue - row_s = row_i*self.cellsize_ratio['L1'] - row_e = (row_i + 1)*self.cellsize_ratio['L1'] + row_s = int(row_i*self.cellsize_ratio['L1']) + row_e = int((row_i + 1)*self.cellsize_ratio['L1']) row_c = int((row_s + row_e)/2) - col_s = col_i*self.cellsize_ratio['L1'] - col_e = (col_i + 1)*self.cellsize_ratio['L1'] + col_s = int(col_i*self.cellsize_ratio['L1']) + col_e = int((col_i + 1)*self.cellsize_ratio['L1']) col_c = int((col_s + col_e)/2) num += 1 if flag == 'any': @@ -301,14 +318,18 @@ def upscale_L1_data(self, L0_array, flag): def downscale_data(self, data_type): pre_L1 = np.zeros((self.nrows['L1'], self.ncols['L1'])) - pre = np.mean(self.pre, axis = 0) - ratio = int(self.cellsize_ratio['L2']/self.cellsize_ratio['L1']) + pre = np.mean(self.pre, axis=0) + L2 = self.cellsize_ratio['L2'] +# L2 = 10 + ratio = int(L2/self.cellsize_ratio['L1']) for row_i in range(0, self.nrows['L1']): for col_i in range(0, self.ncols['L1']): pre_L1[row_i, col_i] = pre[row_i/ratio, col_i/ratio] - self.pre_L1 = np.ma.array( np.mean(self.preEffect, axis = 0), mask = self.mask['L1']) + self.pre_L1 = np.ma.array(pre_L1, mask = self.mask['L1']) +# self.pre_L1 = np.ma.array( np.mean(self.preEffect, axis = 0), mask = self.mask['L1']) + ##-- flow direction functions ------------------------------------------------- @@ -399,6 +420,7 @@ def combine_variables(self, my_list): f_path = self.nml['dir_Out(1)'] + 'mHM_Fluxes_States.nc' var_list = [str(i) for i in readnetcdf(f_path, variables=True)] +# print(var_list) var_dict = {} for elem in my_list: var_dict[elem] = [s for s in var_list if elem + '_L0' in s] @@ -406,4 +428,4 @@ def combine_variables(self, my_list): for sub_elem in var_dict[elem]: tmp += getattr(self, sub_elem ) setattr(self, elem, tmp ) - + \ No newline at end of file diff --git a/post-proc/sas/__init__.py b/post-proc/sas/__init__.py new file mode 100644 index 00000000..4a69f4c1 --- /dev/null +++ b/post-proc/sas/__init__.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python +""" + Python Utilities for computing StorAge-Selection (SAS) functions + + Get help on each function by typing + >>> import sas + >>> help(sas.function) + + + License + ------- + + This file is part of the UFZ Python package. + + Not all files in the package are free software. The license is given in the + 'License' section of the docstring of each routine. + + The package is released under the GNU Lesser General Public License. The + following applies: The SAS Python package is free software: you can + redistribute it and/or modify it under the terms of the GNU Lesser + General Public License as published by the Free Software Foundation, + either version 3 of the License, or (at your option) any later version. + + The SAS Python package is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with the UFZ makefile project (cf. gpl.txt and lgpl.txt). + If not, see . + + Copyright 2015 Falk He"sse + + + History + ------- + Written, FH, Apr 2015 + +""" +#from __future__ import print_function + +# Routines + +from sas_base import SAS +from aux_fun import * +from get_p import * +from get_theta import * +from get_U_num import * +from get_validity_range import * +#from plot_fun import * + +# Information +__author__ = 'Falk Hesse' +__version__ = '0.1.0' +#__revision__ = +__date__ = 'Date: 01.04.2015' + +# Main +#if __name__ == '__main__': +# print('\nMAD Python Package.') +# print("Version {:s} from {:s}.".format(__version__,__date__)) +# print('\nThis is the README file. See als the license file LICENSE.\n\n') +# f = open('README','r') +# for line in f: print(line,end='') +# f.close() diff --git a/post-proc/sas/aux_fun.py b/post-proc/sas/aux_fun.py new file mode 100644 index 00000000..03737151 --- /dev/null +++ b/post-proc/sas/aux_fun.py @@ -0,0 +1,129 @@ + +import numpy as np +from get_theta import * +from get_p import * + +class MapPoint(object): + + def __init__(self, x_array, y_array, lon, lat, x_pos, y_pos): + + self.x_i = (np.abs(x_array - x_pos)).argmin() + self.y_i = (np.abs(y_array - y_pos)).argmin() + self.x = x_array[self.x_i] + self.y = y_array[self.y_i] + self.lon = lon[self.y_i, self.x_i] + self.lat = lat[self.y_i, self.x_i] + + def set_soil_moisture(self, numpy_array): + + self.soil_moisture = numpy_array[:, self.y_i, self.x_i] + + def set_ET(self, numpy_array): + + self.ET = numpy_array[:, self.y_i, self.x_i] + + def set_I(self, numpy_array): + + self.I = numpy_array[:, self.y_i, self.x_i] + + def set_groundwater(self, numpy_array): + + self.groundwater = numpy_array[:, self.y_i, self.x_i] + + def get_ttd(self, U_array, Q_out_array, ET_array, t): + + U = U_array.data[:, self.y_i, self.x_i] + Q_out = Q_out_array.data[:, self.y_i, self.x_i] + ET = ET_array.data[:, self.y_i, self.x_i] + + + t_no = len(t) + tau = np.linspace(1, t_no, t_no) + tau_no = len(tau) + + self.theta = np.zeros(tau_no) + self.eta = np.zeros(tau_no) + self.p_tt_mean = np.zeros(tau_no) + self.p_tt_var = np.zeros(tau_no) + self.mean_p_tt = np.zeros(tau_no) + + for t_in in range(0, t_no): + + self.theta[t_in] = get_theta(U, Q_out, ET, t, t_in) + self.eta[t_in] = get_theta(U, ET, Q_out, t, t_in) + p_tt = get_pdf_tt(U, Q_out, ET, t, t_in) + self.p_tt_mean[t_in] = np.sum(p_tt*np.linspace(1, t_no - t_in, t_no - t_in)) + self.p_tt_var[t_in] = np.sum(p_tt*np.linspace(1, t_no - t_in, t_no - t_in)**2) - self.p_tt_mean[t_in]**2 + self.mean_p_tt[0:t_no-t_in] = self.mean_p_tt[0:t_no-t_in] + p_tt + + self.mean_p_tt = self.mean_p_tt/t_no + +# t_in = 10 +# p_tt = get_pdf_tt(U, Q_out, ET, t, t_in) +# print( np.sum(p_tt*np.linspace(1, t_no - t_in, t_no - t_in)) ) +# plot(self.theta) +# show() + +class DataStruct( object ): + + def __init__(self, numpy_array): + + tmp = numpy_array.shape + y_no = tmp[0] + x_no = tmp[1] + + mask = np.zeros( ( y_no, x_no ), dtype=bool) + + for x_i in range(0, x_no): + for y_i in range(0, y_no): + if numpy_array[y_i][x_i] == -9999: + mask[y_i][x_i] = True + else: + mask[y_i][x_i] = False + + self.field = np.ma.array(numpy_array, mask = mask) +# print(mask) +# print(tmp) + +def get_nearest_value(numpy_array, value): + + nearest_index = (np.abs(numpy_array - value)).argmin() + nearest_value = numpy_array[nearest_index] + + return nearest_value + +def get_nearest_index(numpy_array, value): + + nearest_index = (np.abs(numpy_array - value)).argmin() + + return nearest_index + + +def write_data(data): + + eol = '\n' + eov = '\t' + + f_path = 'elevation.txt' + shape = data.shape + + row_no = shape[1] + col_no = shape[0] + + + f_id = open(f_path, 'w') + + for row_i in range (0, row_no): + + output_line = '' + for col_i in range(0, col_no): +# if data[col_i, row_i] == False: +# output_line = output_line + str(1) + eov +# elif data[col_i, row_i] == True: +# output_line = output_line + str(0) + eov + output_line = output_line + str(data[col_i, row_i]) + eov + + f_id.write( output_line + eol) + + f_id.close() + \ No newline at end of file diff --git a/post-proc/sas/get_U_num.py b/post-proc/sas/get_U_num.py new file mode 100644 index 00000000..ab28eb62 --- /dev/null +++ b/post-proc/sas/get_U_num.py @@ -0,0 +1,12 @@ +#!/usr/bin/env python + +#import numpy as np + +from scipy import integrate + +def get_U_num(t, Q_in, Q_out, U_0): + +# t_no = len(t) +# t = np.linspace(1,t_no,t_no) + + return U_0 + integrate.cumtrapz( Q_in - Q_out, t, initial=0 ) \ No newline at end of file diff --git a/post-proc/sas/get_p.py b/post-proc/sas/get_p.py new file mode 100644 index 00000000..d5d693f1 --- /dev/null +++ b/post-proc/sas/get_p.py @@ -0,0 +1,139 @@ +#!/usr/bin/env python + +import numpy as np +from scipy import integrate + +from get_theta import get_theta +#from get_validity_range import get_validity_range + +def get_p_forward(S, Q_out_1, Q_out_2, t, t_in): + + """ + Computes the forward travel-time pdf within a single reservoir with + one input flux and two output fluxes according to Botter et al. + 'Transport in the hydrological response: Travel time distributions, + soil moisture dynamics, and the old water paradox' WRR (2010), + Equation (35) and Botter et al. 'Catchtment residence time and travel + time distributions: The master equation' GRL (2011), Equation (9) + + + Definition + ---------- + def get_p_forward(S, Q_out_1, Q_out_2, t, t_in): + + + Input + ----- + S volumetric water content of the reservoir + Q_out_1 first output flux of reservoir (contributing to discharge) + Q_out_2 second output flux of reservoir + t chronological time of reservoir dynamcis + t_in point in time when water is entering reservoir + + + Output + ------ + Float array of time-dependent travel-time distribution, representing + the PDF of the random variable T, i.e. the time from entering the + reservoir at t_in to leaving the reservoir as discharge + + + License + ------- + This file is part of the UFZ Python package. + + The UFZ Python package is free software: you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public License as + published by the Free Software Foundation, either version 3 of the + License, or (at your option) any later version. + + The UFZ Python package is distributed in the hope that it will be + useful, but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the UFZ makefile project (cf. gpl.txt and lgpl.txt). + If not, see . + + Copyright 2015 Falk He"sse + + + History + ------- + Written, FH, Mar 2015 """ + + theta = get_theta(S, Q_out_1, Q_out_2, t, t_in) +# eta = get_theta(U, Q_out_2, Q_out_1, t, t_in) + + # needs to be amended later + if theta == 0: + theta = 0.00001 + +# t_end = get_validity_range(t[t_in:], theta + eta) + + f = (Q_out_1[t_in:] + Q_out_2[t_in:])/S[t_in:] + F = integrate.cumtrapz( f, t[t_in:], initial=0 ) + return Q_out_1[t_in:]/(S[t_in:]*theta)*np.exp( -F ) + +def get_p_backward(S, Q_in, Q_out_1, Q_out_2, t, t_ex): + + """ + Computes the backward travel-time pdf within a single reservoir with + one input flux and two output fluxes according to Botter et al. + 'Catchtment residence time and travel time distributions: The master + equation' GRL (2011), Equation (8) + + + Definition + ---------- + def get_p_forward(S, Q_out_1, Q_out_2, t, t_in): + + + Input + ----- + S volumetric water content of the reservoir + Q_in input flux of reservoir (effective precipitation) + Q_out_1 first output flux of reservoir (contributing to discharge) + Q_out_2 second output flux of reservoir + t chronological time of reservoir dynamcis + t_in point in time when water is entering reservoir + + + Output + ------ + Float array of time-dependent travel-time distribution, representing + the PDF of the age distribution of the outflow flux at a given point in + time. + + + License + ------- + This file is part of the UFZ Python package. + + The UFZ Python package is free software: you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public License as + published by the Free Software Foundation, either version 3 of the + License, or (at your option) any later version. + + The UFZ Python package is distributed in the hope that it will be + useful, but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the UFZ makefile project (cf. gpl.txt and lgpl.txt). + If not, see . + + Copyright 2015 Falk He"sse + + + History + ------- + Written, FH, Aug 2015 """ + + + f = (Q_out_1[t_ex::-1] + Q_out_2[t_ex::-1])/S[t_ex::-1] + F = integrate.cumtrapz( f, t[t_ex::-1], initial=0 ) + + return Q_in[t_ex::-1]/S[t_ex::-1]*np.exp(F) diff --git a/post-proc/sas/get_theta.py b/post-proc/sas/get_theta.py new file mode 100644 index 00000000..97145c40 --- /dev/null +++ b/post-proc/sas/get_theta.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python + +import numpy as np + +from scipy import integrate +##from get_indef_integral import get_indef_integral + +def get_theta(U, Q_out_1, Q_out_2, t, t_in): + + """ + Computes the partitioning function theta according to Botter et al. + 'Transport in the hydrological response: Travel time distributions, + soil moisture dynamics, and the old water paradox' WRR (2010), + Equation (B2) + + + Definition + ---------- + def get_theta(U, Q_out_1, Q_out_2, t, t_in): + + + Input + ----- + U volumetric water content of the reservoir + Q_out_1 first output flux of reservoir (contributing to discharge) + Q_out_2 second output flux of reservoir + t chronological time of reservoir dynamcis + t_in time when water is entering reservoir + + + Output + ------ + Float array of time-dependent partitioning function theta, representing + the percentage of Q_in that is entering the reservoir at t_in that is + leaving the reservoir as discharge eventually + + + License + ------- + This file is part of the UFZ Python package. + + The UFZ Python package is free software: you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public License as + published by the Free Software Foundation, either version 3 of the + License, or (at your option) any later version. + + The UFZ Python package is distributed in the hope that it will be + useful, but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the UFZ makefile project (cf. gpl.txt and lgpl.txt). + If not, see . + + Copyright 2015 Falk He"sse + + + History + ------- + Written, FH, Mar 2015 """ + +# t_num = len(t) + tau = t[t_in:] + + f_1 = (Q_out_1[t_in:] + Q_out_2[t_in:])/U[t_in:] + f_2 = Q_out_1[t_in:]/U[t_in:] + + aux_fun = np.exp( - integrate.cumtrapz( f_1, tau, initial=0 )) + theta = np.trapz( f_2*aux_fun, tau ) + + return theta diff --git a/post-proc/sas/get_validity_range.py b/post-proc/sas/get_validity_range.py new file mode 100644 index 00000000..792e92e4 --- /dev/null +++ b/post-proc/sas/get_validity_range.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python + +import numpy as np +#import matplotlib.pyplot as plt + +def get_validity_range(t, f): + + """ + Computes the range of the validity of the approach of Botter et al. + 'Transport in the hydrological response: Travel time distributions, + soil moisture dynamics, and the old water paradox' WRR (2010), + Equation (35) + + + Definition + ---------- + def get_validity_range(t, f): + + + Input + ----- + + t chronological time of reservoir dynamcis + f function whose value should be equal to unity + + + Output + ------ + range within f is actually close enough to unity, i.e. within the + approach of Botter et al. is valid + + + License + ------- + This file is part of the UFZ Python package. + + The UFZ Python package is free software: you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public License as + published by the Free Software Foundation, either version 3 of the + License, or (at your option) any later version. + + The UFZ Python package is distributed in the hope that it will be + useful, but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the UFZ makefile project (cf. gpl.txt and lgpl.txt). + If not, see . + + Copyright 2015 Falk He"sse + + + History + ------- + Written, FH, Mar 2015 """ + + threshold = 0.8 + f = f - threshold + f = np.sign(f) + + range_end = np.where(f==(-1.0)) + range_end = range_end[0] + +# print(range_end[0]) + + return range_end[0] diff --git a/post-proc/sas/plot_fun.py b/post-proc/sas/plot_fun.py new file mode 100644 index 00000000..9f23e3e5 --- /dev/null +++ b/post-proc/sas/plot_fun.py @@ -0,0 +1,123 @@ +#from pylab import * +import ufz +import matplotlib as mpl +import matplotlib.pyplot as plt +#from ufz import readnetcdf +from sas import * + +from cartopy.feature import ShapelyFeature +import cartopy.crs as ccrs +import cartopy.io.shapereader as shpreader + +def plot_routing_network(lon, lat, routing_array): + + ax = axes() + +# ax = axes( projection=ccrs.Mercator()) +# ax.pcolormesh(lon, lat, routing_array[ :, :], transform=ccrs.PlateCarree(), cmap='RdYlGn', edgecolors='face') +# ax.gridlines(draw_labels=True) + +# ax.arrow(0, 0, 1.5, 0.5, head_width=0.05, head_length=0.1, fc='k', ec='k') + ax.annotate('axes center', xy=(.5, .5), xycoords='axes fraction', + horizontalalignment='center', verticalalignment='center') + +# fig = figure(1) +# ax = fig.add_axes(ufz.position(1,1,1, wspace=0.01, hspace=0.01, left=0.14, right=0.97, top=0.95, bottom=0.05), projection=ccrs.Mercator()) + + + show() + +def plot_basin(lon, lat, basin_data, path = None, color_map = None): + + mpl.rc('text', usetex=True) + mpl.rc('text.latex', unicode=True) + +# shpf_river_unstrut = './shape/unstrut_river.shp' +# shpf_basin_unstrut = './shape/unstrut.shp' +# shpf_thubasin = './shape/watershed_WGS84.shp' +# shpf_thubasin = './shape/Grenzen_WGS84' + +# if path: +# x = readnetcdf(path, var='easting') +# y = readnetcdf(path, var='northing') + +# nagelstedt = MapPoint(x, y, lon, lat, 4409672, 5664533) +# cammerbach = MapPoint(x, y, lon, lat, 4392737, 5668127) +# suthbach = MapPoint(x, y, lon, lat, 4396321, 5667869) + +# reckenbuhl = MapPoint(x, y, lon, lat, 4387837, 5664518) +# bechstedt = MapPoint(x, y, lon, lat, 4388942, 5664961) +# heuberg = MapPoint(x, y, lon, lat, 4390215, 5665475) +# fuchsloch = MapPoint(x, y, lon, lat, 4391430, 5665590) +# ruspelsweg = MapPoint(x, y, lon, lat, 4392925, 5666035) + +# test = np.append(lon.data, np.zeros((1, 90)), 0) +# print(test.shape) + + fig = plt.figure(1) + + font = {'size' : 36} + mpl.rc('font', **font) +# cmap = mpl.cm.get_cmap('YlOrRd') + mp = fig.add_axes(ufz.position(1,1,1), projection=ccrs.Mercator()) + pcm = mp.pcolormesh(lon, lat, basin_data, transform=ccrs.PlateCarree(), + cmap=color_map, edgecolors='face') + ixticks = np.round(np.arange((lon.min()), (lon.max()), 0.1)*10)/10 + iyticks = np.round(np.arange((lat.min()), (lat.max()), 0.1)*10)/10 + + mp.set_yticks(iyticks, crs=ccrs.PlateCarree()) + ynames = [ r'$\mathrm{'+str((i))+'\,^{\circ} N}$' for i in iyticks] + yticknames = plt.setp(mp, yticklabels=ynames) + mp.set_xticks(ixticks, crs=ccrs.PlateCarree()) + xnames = [] + for i in list(ixticks): + if (i < 0): + xnames.append(r'$\mathrm{'+str((i))+'\,^{\circ} W}$') + elif (i == 0): + xnames.append(r'$\mathrm{'+str((i))+'\,^{\circ}}$') + elif (i > 0): + xnames.append(r'$\mathrm{'+str((i))+'\,^{\circ} E}$') + xticknames = plt.setp(mp, xticklabels=xnames) +# plt.axis([lon.min(), lon.max(), lat.min(), lat.max()]) +# plt.xticks(ixticks) +# mp.tick_params(labelbottom=True, labelleft=False) + mp.grid(True) + cb = plt.colorbar(pcm, orientation='vertical') +# mpl.rcParams['xtick.labelsize'] = font_size + +# fig = plt.figure(1) +# ax = fig.add_axes(ufz.position(1,1,1, wspace=0.01, hspace=0.01, left=0.14, right=0.97, top=0.95, bottom=0.05), projection=ccrs.Mercator()) + +# reader = shpreader.Reader(shpf_river_unstrut) +# riv_unstrut = ShapelyFeature(reader.geometries(), ccrs.PlateCarree()) +# ax.add_feature(riv_unstrut, facecolor='none', edgecolor='#045a8d', lw=2.5) + +# reader = shpreader.Reader(shpf_basin_unstrut) +# unstrut_basin = ShapelyFeature(reader.geometries(), ccrs.PlateCarree()) +# ax.add_feature(unstrut_basin, facecolor='none', edgecolor='0.5', lw=1.0) + +# reader = shpreader.Reader(shpf_thubasin) +# thube_basin = ShapelyFeature(reader.geometries(), ccrs.PlateCarree()) +# ax.add_feature(thube_basin, facecolor='none', edgecolor='k', lw=1.0) + +# if path: +## ax.scatter(nagelstedt.lon, nagelstedt.lat, transform=ccrs.PlateCarree()) +# ax.scatter(reckenbuhl.lon, reckenbuhl.lat, c='b', +# transform=ccrs.PlateCarree()) +# ax.scatter(bechstedt.lon, bechstedt.lat, c='g', +# transform=ccrs.PlateCarree()) +# ax.scatter(heuberg.lon, heuberg.lat, c='r', +# transform=ccrs.PlateCarree()) +# ax.scatter(fuchsloch.lon, fuchsloch.lat, c='c', +# transform=ccrs.PlateCarree()) +# ax.scatter(ruspelsweg.lon, ruspelsweg.lat, c='m', +# transform=ccrs.PlateCarree()) + +# ax.set_autoscaley_on(False) +# ax.set_ylim([0,10]) +# ax.set_extent( [ 9.5, 12.5, 50.5, 51.7 ] ) + + +# ax.title('mean travel time', fontsize=font_size) + + plt.show() \ No newline at end of file diff --git a/post-proc/sas/plot_study_domain.py b/post-proc/sas/plot_study_domain.py new file mode 100644 index 00000000..eac00f89 --- /dev/null +++ b/post-proc/sas/plot_study_domain.py @@ -0,0 +1,339 @@ +#!/usr/bin/env python +# +# purpose: creating netdcf for small subbasins in catch having the same +# georeferencing (lower left corner) +# +# input: catch mask, lut subbasin masks, subbasin mask (ASCII files) +# +# restrictions: only integer values in all input data excpected/accepted +# +# created by Matthias Zink Feb. 2014 +# +import numpy as np # array manipulation +import ufz +import matplotlib as mpl +import cartopy.crs as ccrs +import cartopy.feature as cfeature +import cartopy.io.shapereader as shpreader +from cartopy.feature import ShapelyFeature +from cartopy.io import srtm +# +# ------------------------------------------------------------------------- +# Command line arguments +# +pdffile = '' +pngfile = '' +import optparse +parser = optparse.OptionParser(usage='%prog [options]', + description="Plotting file following template of MC.") +parser.add_option('-p', '--pdffile', action='store', dest='pdffile', type='string', + default=pdffile, metavar='File', + help='Name of pdf output file (default: open X-window).') +parser.add_option('-g', '--pngfile', action='store', + default=pngfile, dest='pngfile', metavar='pngfile', + help='Name basis for png output files (default: open screen window).') +(opts, args) = parser.parse_args() + +pdffile = opts.pdffile +pngfile = opts.pngfile +del parser, opts, args + +if (pdffile != '') & (pngfile != ''): + raise ValueError('PDF and PNG are mutually exclusive. Only either -p or -g possible.') + +# ------------------------------------------------------------------------- +# Customize plots +# +if (pdffile == ''): + if (pngfile == ''): + outtype = 'x' + else: + outtype = 'png' +else: + outtype = 'pdf' + + +# Plot - paper_plots, but also all if not otherwise defined +nrow = 1 # # of rows per figure +ncol = 1 # # of columns per figure +hspace = 0.01 # x-space between plots +wspace = 0.01 # y-space between plots +figleft = 0.14 # left border ogf the figure +figright = 0.97 # right border ogf the figure +figbottom = 0.05 # lower border ogf the figure +figtop = 0.95 # upper border ogf the figure +textsize = 18 # Standard text size +dt = 4 # # of hours between tick marks on plots +dxabc = 0.90 # % shift from left y-axis of a,b,c,... labels +dyabc = 0.90 # % shift from lower x-axis of a,b,c,... labels +dyabcdown = 0.05 # y-shift if abc in lower right corner +lwidth = 0.5 # linewidth +elwidth = 1.0 # errorbar line width +alwidth = 1.0 # axis line width +msize = 10. # marker size +mwidth = 4 # marker edge width +bxwidth = 0.85 # boxlplot width +# color: 'b'|'g'|'r'|'c'|'m'|'y'|'k'|'w' +# 'blue'|'green'|'red'|'cyan'|'magenta'|'yellow'|'black'|'white' +# hex string '#eeefff' | RGB tuple (1,0.5,1) | html names 'burlywod', 'chartreuse', ... +# grayscale intensity, e.g. '0.7', 'k'='0.0' +mcol1 = '#67A9CF' # color of second markers +mcol2 = '#A1D99B' # color of third markers +mcol3 = '#EF8A62' # primary line colour +mcol4 = 'r' +lcol2 = '0.5' # color of second lines +lcol3 = '0.0' # color of third lines + +llxbbox = 0.5 # x-anchor legend bounding box +llybbox = 0.87 # y-anchor legend bounding box +llrspace = 0.02 # spacing between rows in legend +llcspace = 1.0 # spacing between columns in legend +llhtextpad = 0.2 # the pad between the legend handle and text +llhlength = 0.9 # the length of the legend handles +frameon = True # if True, draw a frame around the legend. If None, use rc +llxbbox2 = 0.60 # Tight bounding of symbol and text (w/o lines) +llhtextpad2= 0. # " +llhlength2 = 1.0 # " + +# PNG +dpi = 300 +transparent = False +bbox_inches = 'tight' +pad_inches = 0.02 +# +if (outtype == 'pdf'): + mpl.use('PDF') # set directly after import matplotlib + import matplotlib.pyplot as plt + from matplotlib.backends.backend_pdf import PdfPages + # Customize: http://matplotlib.sourceforge.net/users/customizing.html + mpl.rc('ps', papersize='a4', usedistiller='xpdf') # ps2pdf +elif (outtype == 'png'): + mpl.use('Agg') # set directly after import matplotlib + import matplotlib.pyplot as plt + mpl.rc('text.latex', unicode=True) + mpl.rc('savefig', dpi=dpi, format='png') +else: + import matplotlib.pyplot as plt +# +# scale figsize for pdf +latexwidth = 470 # invoke latex document textwidth with \the\textwidth +factor = 1.00 # factor used in latex document like width=factor\textwidth + +fig_width_in = latexwidth * factor * 0.01388888 # figure width in inches * inches per point +fig_height_in = fig_width_in #* 1.318 # figure height in inches * golden ratio + +mpl.rc('figure', figsize=( fig_width_in, fig_height_in)) # half side A4potrait, golden ratio + +#mpl.rc('figure', figsize=( 8, 5)) # for AGU 2012 Poster +#for paper half of the size because of 2 columns +# one column of 2 column paper +#mpl.rc('figure', figsize=( 8.27/2., 8.27/2./1.618)) # half side A4potrait, golden ratio +#mpl.rc('figure', figsize=( 8.27, 11.69)) # a4 portrait +#mpl.rc('figure', figsize=(11.69, 8.27)) # a4 landscape +#mpl.rc('figure', figsize=( 8.27, 11.69/2)) # a4 portrait +mpl.rc('font', **{'family':'serif','serif':['times']}) +mpl.rc('font', size=textsize) +mpl.rc('legend', fontsize=textsize) +mpl.rc('lines', linewidth=lwidth, color='black') +mpl.rc('axes', linewidth=alwidth, labelcolor='black') +mpl.rc('path', simplify=False) # do not remove +mpl.rc('text', usetex=True) +mpl.rc('text.latex', unicode=True) + +############################################################################################## +if (outtype == 'pdf'): + print 'Plot PDF ', pdffile + pdf_pages = PdfPages(pdffile) +elif (outtype == 'png'): + print('Plot PNG ', pngfile) +else: + print 'Plot X' + +figsize = mpl.rcParams['figure.figsize'] +ifig = 0 + +shpf_eu_adm = '../wb_paper_de/study_domain/data/eu_admin_bound.shp' +shpf_unstr = '../wb_paper_de/study_domain/data/BuLaender.shp' +shpf_river = '../wb_paper_de/study_domain/data/neckar_mulde_saale.shp' +shpf_river_unstrut = './shape/unstrut_river.shp' +shpf_basin_unstrut = './shape/unstrut.shp' +shpf_thubasin = './shape/watershed_WGS84.shp' +shpf_thubasin = './shape/Grenzen_WGS84' +txt_ecstat = './EC_wb_paper.txt' + +demfile = '../wb_paper_de/study_domain/data/dem_de/dem_1km.nc' + +inpath = '../wb_paper_de/study_domain/data/watersheds/' + +catchlist = ['saale','weser'] +label_cat = ['Saale','Weser'] +lon_cat = [ 10.90, 9.20] +lat_cat = [ 51.85, 52.40] +# MAPPLOT ######################################### +ifig += 1 +print 'Plot - Fig ', ifig +fig = plt.figure(ifig) +iplot = 1 + +# data PROCESSING ############################################################# +data = ufz.readnetcdf(demfile, var='mask')[0,:,:] +lons = ufz.readnetcdf(demfile, var='lon') +lats = ufz.readnetcdf(demfile, var='lat') +dem = np.ma.array(data, mask=~(data>-9999.)) +shaded = np.ma.array(srtm.add_shading(data*10,315.,45.), mask=~(data>-9999.)) +test = dem / 5. + shaded +smooth = np.ma.array(ufz.savitzky_golay2d(shaded,7,2), mask=~(data>-9999.)) + +# ec stations +data = ufz.sread(txt_ecstat, skip=1, strarr=True) +data_header = ufz.sread(txt_ecstat, skip=1, header=True, strarr=True) + +lat_ec = data[:,np.where(data_header=='Latitude')[0][0]].astype('float') +lon_ec = data[:,np.where(data_header=='Longitude')[0][0]].astype('float') +label_ec = data[:,np.where(data_header=='Stat_ID')[0][0]] +id_ec = data[:,np.where(data_header=='Id')[0][0]] +# plotting ############################################################# + +mp = fig.add_axes(ufz.position(nrow,ncol,1, wspace=0.01, hspace=0.01, + left=figleft, right=figright, top=figtop, bottom=figbottom), + projection=ccrs.Mercator()) + +# coordinate transformation +transform = ccrs.PlateCarree()._as_mpl_transform(mp) +# coloring +#cmap = 'binary' +#cmap2 = 'gist_gray' +colors = ufz.get_brewer('gsdtol',rgb=True)[1:] +cmap = mpl.colors.ListedColormap(colors[:4:-1]) +#cmap = mpl.colors.ListedColormap(ufz.get_brewer('OceanLakeLandSnow',rgb=True)[1:]) +cmap2 = mpl.colors.ListedColormap(colors) + +breaks = list(np.linspace(-50,1000,22)) + [1200,1400,1600,1800,2000,2200] +norm = mpl.colors.BoundaryNorm(breaks, cmap.N) +norm2 = mpl.colors.BoundaryNorm(np.linspace(-20,250,100),cmap.N) + +# transparency +cmap2._init() +cmap2._lut[:-3:,-1] = 0.99 +cmap._init() +cmap._lut[:-3:,-1] = 0.015 + + +#pcm2 = mp.pcolormesh(lons, lats, test, transform=ccrs.PlateCarree(), cmap=cmap, norm=norm2, +# edgecolors='face') +# pcm1 = mp.pcolormesh(lons, lats, smooth, transform=ccrs.PlateCarree(), cmap=cmap2, +# edgecolors='face')#, antialiased=True) +# pcm3 = mp.pcolormesh(lons, lats, dem, transform=ccrs.PlateCarree(), cmap=cmap, norm=norm, +# edgecolors='face') # alpha=0.08, +#pcm3 = mp.contourf(lons, lats, dem, transform=ccrs.PlateCarree(), +# levels=breaks, cmap=cmap, norm=norm, alpha=0.5)#, antialiased=True) + +#plt.show() +#stop +# catchments +for icat in range(len(catchlist)): + # federal countries, bundeslaender + if (icat==5 or icat==2): + color='r' + elif (icat==0 or icat==4): + color='g' + else: + color='b' + reader = shpreader.Reader(inpath + 'basin_' + catchlist[icat] + '.shp') + fed_count = ShapelyFeature(reader.geometries(), ccrs.PlateCarree()) + mp.add_feature(fed_count, facecolor=color, alpha=0.10, lw=0.5, edgecolor='0') #, alpha=0.75) + # catchment labels + mp.annotate(label_cat[icat],(lon_cat[icat], lat_cat[icat]), xycoords=transform, + va='center', ha='left', color=color) + + +# eddy stations +mp.plot(lon_ec, lat_ec, ls='None',marker='s',mfc='None', mec='k',mew=1.2, ms=3,transform=ccrs.PlateCarree()) +for i in range(len(lon_ec)): + if ( label_ec[i] == 'E1'): + xp=lon_ec[i]+0.03; yp=lat_ec[i]#-0.15 + va='top'; ha='left' + elif ( label_ec[i] == 'E2'): + xp=lon_ec[i]; yp=lat_ec[i]#-0.15 + va='top'; ha='center' + elif ( label_ec[i] == 'E3'): + xp=lon_ec[i]-0.1; yp=lat_ec[i] + va='bottom'; ha='right' + elif ( label_ec[i] == 'E4'): + xp=lon_ec[i]; yp=lat_ec[i]#-0.15 + va='top'; ha='right' + mp.annotate(id_ec[i].split('-')[1],(xp, yp), xycoords=transform, va=va, ha=ha, fontsize=0.95*textsize) + #plt.text(lon_ec[i], lat_ec[i], label_ec[i], transform=ccrs.PlateCarree()) + +# mp.set_global() +ixticks = np.arange( 2,17.5, 1) +iyticks = np.arange( 46, 60, 1) + +mp.set_yticks(iyticks, crs=ccrs.PlateCarree()) +mp.set_xticks(ixticks, crs=ccrs.PlateCarree()) + +xnames = [] +for i in list(ixticks): + if (i < 0): + xnames.append(r'$\mathrm{'+str(int(i))+'\,^{\circ} W}$') + elif (i == 0): + xnames.append(r'$\mathrm{'+str(int(i))+'\,^{\circ}}$') + elif (i > 0): + xnames.append(r'$\mathrm{'+str(int(i))+'\,^{\circ} E}$') +ynames = [ r'$\mathrm{'+str(i)+'\,^{\circ} N}$' for i in iyticks] + +yticknames = plt.setp(mp, yticklabels=ynames) +plt.setp(yticknames, fontsize=.9*textsize) # fontsize='medium') +xticknames = plt.setp(mp, xticklabels=xnames) +plt.setp(xticknames, fontsize=.9*textsize) + +# properties and features +mp.set_extent([7.95, 12.95, 50.1, 53.95]) #, crs=ccrs.Mercator()) + +# geographical features +ocean = cfeature.NaturalEarthFeature(category='physical', name='ocean', + scale='50m', facecolor=cfeature.COLORS['water']) +mp.add_feature(ocean, edgecolor='none') +# state_borders = cfeature.NaturalEarthFeature(category='cultural', name='admin_0_countries', +# scale='10m', facecolor='none') +# mp.add_feature(state_borders, edgecolor='k', lw=0.5) +reader = shpreader.Reader(shpf_eu_adm) +state_borders = ShapelyFeature(reader.geometries(), ccrs.PlateCarree()) +mp.add_feature(state_borders, facecolor='none', edgecolor='k', lw=0.5) +# rivers +rivers = cfeature.NaturalEarthFeature(category='physical', name='rivers_lake_centerlines', + scale='10m', facecolor='none') +mp.add_feature(rivers, edgecolor='#045a8d', lw=0.5) + +reader = shpreader.Reader(shpf_river) +add_rivers = ShapelyFeature(reader.geometries(), ccrs.PlateCarree()) +mp.add_feature(add_rivers, facecolor='none', edgecolor='#045a8d', lw=0.5) + +reader = shpreader.Reader(shpf_basin_unstrut) +unstrut_basin = ShapelyFeature(reader.geometries(), ccrs.PlateCarree()) +mp.add_feature(unstrut_basin, facecolor='0.5', edgecolor='0.5', lw=1.0) + +reader = shpreader.Reader(shpf_thubasin) +thube_basin = ShapelyFeature(reader.geometries(), ccrs.PlateCarree()) +mp.add_feature(thube_basin, facecolor='none', edgecolor='k', lw=1.0) + +reader = shpreader.Reader(shpf_river_unstrut) +river_unstrut = ShapelyFeature(reader.geometries(), ccrs.PlateCarree()) +mp.add_feature(river_unstrut, facecolor='none', edgecolor='#045a8d', lw=0.5) + + +if (outtype == 'pdf'): + pdf_pages.savefig(fig) + plt.close() +elif (outtype == 'png'): + fig.savefig(pngfile+str(ifig).zfill(3)+'.png', transparent=transparent, bbox_inches=bbox_inches, pad_inches=pad_inches) + plt.close(fig) +else: + plt.show() + +if (outtype == 'pdf'): + pdf_pages.close() +elif (outtype == 'png'): + pass + + diff --git a/post-proc/sas/sas_base.py b/post-proc/sas/sas_base.py new file mode 100644 index 00000000..5d928385 --- /dev/null +++ b/post-proc/sas/sas_base.py @@ -0,0 +1,359 @@ + +import numpy as np + +from mhm import MHM +from sas.get_p import get_p_forward +from sas.get_p import get_p_backward +from sas.get_theta import get_theta + +class SAS(MHM): + + def __init__(self, model_path, file_name = 'mhm.nml', + rel_path_name = False): + MHM.__init__(self, model_path, file_name, rel_path_name) + + def map_variables(self): + + var_dict = {} + var_dict['x_3_array'] = ['SWC'] + var_dict['x_5_array'] = ['unsatSTW'] + var_dict['x_6_array'] = ['satSTW'] + var_dict['ET_array'] = ['aET'] + var_dict['QI'] = ['QIf', 'QIs'] + var_dict['BF_array'] = ['QB'] + var_dict['R_array'] = ['recharge'] +# var_dict['root'] = ['SWC'] +# var_dict['soil'] = ['SWC', 'unsatSTW'] +# var_dict['subsurface'] = ['SWC', 'unsatSTW', 'satSTW'] +# var_dict['Q_in_array'] = ['preEffect'] + + for elem in var_dict: + tmp = 0 + for sub_elem in var_dict[elem]: + tmp += getattr(self, sub_elem ) + setattr(self, elem, tmp ) + +##-- determining pdf's -------------------------------------------------------- + + def get_p(self, p_type = 'forward', flag = 'soil', flux_type = 'Q', + time_series = 'all', fit_p = True): + + self.map_variables() + self.set_spatial_mask('L1') + + if time_series == 'all': + self.time_series_flag = np.ones((self.t_no)) + elif time_series == 'storm_events': + self.get_storm_events() + elif time_series == 'dry_years': + self.get_wet_dry_years(time_series) + elif time_series == 'wet_years': + self.get_wet_dry_years(time_series) + + if flag == 'root': + S = self.SWC + Q_1 = self.J + Q_2 = self.aET + elif flag == 'soil': + S = self.SWC + self.unsatSTW + Q_1 = self.QIf + self.QIs + self.recharge + Q_2 = self.aET + elif flag == 'subsurface': + S = self.SWC + self.unsatSTW + self.satSTW + Q_1 = self.QIf + self.QIs + self.OB + Q_2 = self.aET + + if p_type == 'forward': + if flux_type == 'Q': + self.get_p_Q_forward(S, Q_1, Q_2) + elif flux_type == 'ET': + self.get_p_Q_forward(S, Q_2, Q_1) + elif p_type == 'backward': + if flux_type == 'Q': + self.get_p_Q_backward(flag) + elif flux_type == 'ET': + self.get_p_Q_backward(flag) + elif p_type == 'theta': + self.get_p_theta(S, Q_1, Q_2) + + if fit_p == True: + self.fit_p() + + def get_p_Q_forward(self, S, Q_1, Q_2): + + self.mean_p_Q = np.zeros((self.t_no, self.y_no, self.x_no)) + self.my_p_Q = np.zeros((self.t_no)) + + for x_i in range(0, self.x_no): + for y_i in range(0, self.y_no): + if self.mask['L1_special'][y_i, x_i]: + continue + p_Q = np.zeros((self.t_no)) + counter = 0 + for t_in in range(0, self.t_no): + T = self.t_no - t_in + if self.time_series_flag[t_in] > 0: + p_Q[:T] += get_p_forward(S[:, y_i, x_i], + Q_1[:, y_i, x_i], + Q_2[:, y_i, x_i], + self.t, t_in) + counter += 1 + if (t_in == 10) and (x_i == 40) and (y_i == 40): +# print(t_in) + self.my_p_Q = get_p_forward(S[:, y_i, x_i], + Q_1[:, y_i, x_i], + Q_2[:, y_i, x_i], + self.t, t_in) + if (x_i == 40) and (y_i == 40): + self.p_Q = p_Q/counter + self.mean_p_Q[:,y_i,x_i] = p_Q/counter + + def get_p_Q_backward(self, flag): + + self.mean_p_Q = np.zeros((self.t_no, self.y_no, self.x_no)) + + for x_i in range(0, self.x_no): + for y_i in range(0, self.y_no): + if self.mask['L1'][y_i, x_i]: + continue + if flag == 'root': + S = self.x_3_array.data[:, y_i, x_i] + J = self.Q_in_array[:, y_i, x_i] + Q_1 = self.Q_out_33_array[:, y_i, x_i] + Q_2 = self.ET_array[:, y_i, x_i] + elif flag == 'soil': + S = self.SWC[:, y_i, x_i] + self.unsatSTW[:, y_i, x_i] + J = self.preEffect[:, y_i, x_i] + Q_1 = self.QI[:, y_i, x_i] + self.recharge[:, y_i, x_i] + Q_2 = self.aET[:, y_i, x_i] + p_Q = np.zeros((self.t_no)) +# counter = 0 + for t_ex in range(1, self.t_no): + p_Q[:t_ex+1] += get_p_backward(S, J, Q_1, Q_2, self.t, t_ex) + self.mean_p_Q[:,y_i,x_i] = p_Q/(self.t_no - 1) +# if (x_i == 20) and (y_i == 20): +# self.p_Q = p_Q/(self.t_no - 1) + + def get_p_rt(self): + + self.map_variables() + + self.p_rt_tot = np.zeros((self.t_no, self.t_no)) + self.p_a_tot = np.zeros((self.t_no, self.t_no)) + self.omega = np.zeros((self.t_no, self.t_no)) + + soil_moisture = self.SWC + self.unsatSTW + discharge = self.QI + self.recharge +# soil_moisture_tot = np.mean(np.mean(soil_moisture, axis = 1), axis = 1) +# self.soil_moisture_rel = np.zeros((self.t_no, self.y_no, self.x_no)) + self.soil_moisture_rel = np.ma.copy((self.SWC)) + self.discharge_rel = np.ma.copy((self.SWC)) + for t_i in range(0, self.t_no): + soil_moisture_tot = np.sum(np.sum(soil_moisture[t_i,:,:])) + self.soil_moisture_rel.data[t_i, :, :] = soil_moisture[t_i, :, :]/soil_moisture_tot + discharge_tot = np.sum(np.sum(discharge[t_i,:,:])) + self.discharge_rel.data[t_i, :, :] = discharge[t_i, :, :]/discharge_tot +# print(soil_moisture.shape) +# print(soil_moisture_tot.shape) +# discharge = self.QIf + self.QIs + self.recharge +# discharge_tot = np.mean(np.mean(discharge, axis = 1), axis = 1) + + for t_ex in range(1, self.t_no): + p_Q = np.zeros((self.t_no)) + p_rt = np.zeros((self.t_no)) + p_a = np.zeros((self.t_no)) + for x_i in range(0, self.x_no): + for y_i in range(0, self.y_no): + if self.mask['L1'][y_i, x_i]: + continue + S = self.SWC[:, y_i, x_i] + self.unsatSTW[:, y_i, x_i] + J = self.preEffect[:, y_i, x_i] + Q_1 = self.QI[:, y_i, x_i] + self.recharge[:, y_i, x_i] + Q_2 = self.aET[:, y_i, x_i] + +# T = self.t_no - t_ex +# p_Q[:T] = get_p_forward(S, Q_1, Q_2, self.t, t_ex) + p_Q[:t_ex+1] = get_p_backward(S, J, Q_1, Q_2, self.t, t_ex) +# ratio = soil_moisture[t_ex, y_i, x_i]/soil_moisture_tot[t_ex] +# print(ratio.shape) +# print(p_Q.shape) + p_rt += p_Q*self.soil_moisture_rel[t_ex,y_i,x_i] + p_a += p_Q*self.discharge_rel[t_ex,y_i,x_i] + + self.p_rt_tot[:,t_ex] = p_rt + self.p_a_tot[:,t_ex] = p_a +# self.omega[:,t_ex] = p_a/p_rt + + def get_p_theta(self, S, Q_1, Q_2): + + self.theta = np.zeros((self.t_no, self.y_no, self.x_no)) + self.eta = np.zeros((self.t_no, self.y_no, self.x_no)) + + for x_i in range(0, self.x_no): + for y_i in range(0, self.y_no): + if self.mask['L1_special'][y_i, x_i]: + continue + for t_in in range(0, self.t_no): + self.theta[t_in, y_i, x_i] = get_theta(S[:, y_i, x_i], Q_1[:, y_i, x_i], Q_2[:, y_i, x_i], self.t, t_in) + self.eta[t_in, y_i, x_i] = get_theta(S[:, y_i, x_i], Q_2[:, y_i, x_i], Q_1[:, y_i, x_i], self.t, t_in) + + + +##-- fitting pdf to analytical functions -------------------------------------- + + def fit_p(self, func = 'exp'): + + from scipy.optimize import curve_fit + from scipy.special import gamma + + def exp_dist(x, a, b): + return a*np.exp(-x/b) + def gamma_dist(x, a, b): + return (1/b)**a/gamma(a)*x**(a-1)*np.exp(-x/b) + + func = func[:3].lower() + tmp = np.zeros((self.nrows['L1'], self.ncols['L1'] )) + self.mean_tt_opt = np.ma.array(tmp, mask = self.mask['L1_special']) + + for x_i in range(0, self.x_no): + for y_i in range(0, self.y_no): + if self.mask['L1_special'][y_i][x_i]: + continue + if func == 'exp': + popt, pcov = curve_fit(exp_dist, self.t[:100], + self.mean_p_Q[:100, y_i, x_i]) + if func == 'gam': + popt, pcov = curve_fit(gamma_dist, self.t[:100], + self.mean_p_Q[:100, y_i, x_i]) + self.mean_tt_opt[y_i, x_i] = popt[1] + +##-- determining mask for evaluating the TTDs --------------------------------- + + def set_spatial_mask(self, mask_array, value = 0): + + data_level = 'L1' + if 'L1_special' not in self.mask: + self.mask['L1_special'] = np.copy(self.mask[data_level]) + +# index = np.where( self.landcover_L1 == mask_id) + if mask_array == 'non_urban': + data_type = 'landcover' + self.upscale_data(data_type, data_level) + self.mask['L1_special'][np.where(self.landcover_L1 == 2)] = 1 + if mask_array == 'urban': + data_type = 'landcover' + self.upscale_data(data_type, data_level, flag = 'mean') +# self.mask['L1_special'][np.where(self.landcover_L1 == 1)] = 1 +# self.mask['L1_special'][np.where(self.landcover_L1 == 3)] = 1 + self.mask['L1_special'][np.where(self.mask['L1_special'] == 0)] = 1 + self.mask['L1_special'][np.where(self.landcover_L1 == 2)] = 0 + if mask_array == 'forest': + data_type = 'landcover' + self.upscale_data(data_type, data_level) + self.mask['L1_special'][np.where(self.landcover_L1 == 2)] = 1 + self.mask['L1_special'][np.where(self.landcover_L1 == 3)] = 1 + if mask_array == 'grass': + data_type = 'landcover' + self.upscale_data(data_type, data_level) + self.mask['L1_special'][np.where(self.landcover_L1 == 1)] = 1 + self.mask['L1_special'][np.where(self.landcover_L1 == 2)] = 1 + if mask_array == 'lai_class_1': + data_type = 'lai_class' + self.upscale_data(data_type, data_level) + self.mask['L1_special'][:,:] = 1 + self.mask['L1_special'][np.where(self.lai_class_L1 == 1)] = 0 + if mask_array == 'lai_class_10': + data_type = 'lai_class' + self.upscale_data(data_type, data_level) + self.mask['L1_special'][:,:] = 1 + self.mask['L1_special'][np.where(self.lai_class_L1 == 10)] = 0 + if mask_array == 'soil_class': +# self.mask['L1_special'] = np.copy(self.mask[data_level]) + data_type = 'soil_class' + self.upscale_data(data_type, data_level) + self.mask['L1_special'][:,:] = 1 + self.mask['L1_special'][np.where(self.soil_class_L1 == value)] = 0 + + +# self.mask['L1_special'] = np.copy(self.mask['L1']) + +##-- determining weather extremes --------------------------------------------- + + def get_wet_dry_years(self, time_series): + + time_step = 12 + self.averaged_x_3 = np.zeros((self.t_no)) +# self.dry_year = np.zeros((self.t_no)) +# self.wet_year = np.zeros((self.t_no)) + self.time_series_flag = np.zeros((self.t_no)) + +# tmp = np.mean(self.SWC, axis=1) + self.mean_soil_moisture = np.mean(np.mean(self.SWC, axis=1), axis=1) + my_soil_moisture = (self.mean_soil_moisture - np.mean(self.mean_soil_moisture))/np.std(self.mean_soil_moisture) + + for t_i in range(0, self.t_no/time_step): + a = t_i*time_step + e = t_i*time_step + time_step + tmp = np.mean(my_soil_moisture[a:e]) + tmp = np.ones(( time_step ))*tmp + self.averaged_x_3[a:e] = tmp + if time_series == 'wet_years': + if tmp[0] > 0.605: + self.time_series_flag[a:e] = np.ones(( time_step )) + else: + self.time_series_flag[a:e] = np.zeros(( time_step )) + elif time_series == 'dry_years': + if tmp[0] < -0.426: + self.time_series_flag[a:e] = np.ones(( time_step )) + else: + self.time_series_flag[a:e] = np.zeros(( time_step )) + + def get_storm_events(self): + + self.time_series_flag = np.zeros((self.t_no)) + self.mean_preEffect = np.mean(np.mean(self.preEffect, axis=1), axis=1) + + for t_i in range(0, self.t_no): + if self.mean_preEffect[t_i] > 100: + self.time_series_flag[t_i] = 1 + +##-- compressing time series -------------------------------------------------- + + def map_to_month(self, flux_data_list, state_data_list): + chunk_no = 33 + self.compress_flux(flux_data_list, chunk_no) + self.compress_state(state_data_list, chunk_no) + self.t_no = self.t_no/chunk_no + self.t = range(1, self.t_no + 1) + + + def compress_flux(self, data_list, chunk_no): + for data_i in data_list: + setattr(self, data_i, self.sum_chunk(getattr(self, data_i), + chunk_no, 0)) + + def compress_state(self, data_list, chunk_no): + for data_i in data_list: + setattr(self, data_i, self.mean_chunk(getattr(self, data_i), + chunk_no, 0)) + + def sum_chunk(self, x, chunk_size, axis): + + shape = x.shape + if axis < 0: + axis += x.ndim + shape = shape[:axis] + (-1, chunk_size) + shape[axis+1:] + x = x.reshape(shape) + return x.sum(axis=axis+1) + + def mean_chunk(self, x, chunk_size, axis): + + shape = x.shape + if axis < 0: + axis += x.ndim + shape = shape[:axis] + (-1, chunk_size) + shape[axis+1:] + x = x.reshape(shape) + return x.mean(axis=axis+1) + + + \ No newline at end of file diff --git a/post-proc/test_mhm_sas_class.py b/post-proc/test_mhm_sas_class.py new file mode 100644 index 00000000..9c7e5f29 --- /dev/null +++ b/post-proc/test_mhm_sas_class.py @@ -0,0 +1,46 @@ +#!/usr/bin/env python + +#import os +#import sys + +#sys.path.append( os.environ['HOME'] + '/Source/lib/chs-svn/PYTHON_chs_lib/' ) +#sys.path.append( './sas/' ) + +#import ufz +import numpy as np +from sas import SAS +#from sas.plot_fun import plot_basin +#import matplotlib as mpl +#import matplotlib.pyplot as plt +#import cartopy.crs as ccrs +#import cartopy.feature as cfeature + +#-- paths --------------------------------------------------------------------- + +project_path = '../' + +#-- initialize mhm class ------------------------------------------------------ + +mhm = SAS(model_path = project_path, file_name = 'mhm.nml', + rel_path_name = True) + +#-- import mhm data (input and output) ---------------------------------------- + +mhm.import_data('lat_lon') +mhm.import_data('slope') +mhm.import_data('soil_class') +mhm.import_data('lai_class') +mhm.import_data('dem') +mhm.import_data('pre') + +mhm.import_data('states_and_fluxes', 'all') +mhm.combine_variables(['SWC']) + +#-- postprocessing mhm -------------------------------------------------------- + +mhm.get_p(p_type = 'forward', time_series = 'all') + +#-- plotting ------------------------------------------------------------------ + +#plot_basin(mhm.lon, mhm.lat, np.mean(mhm.recharge, axis = 0)) +#plot_basin(mhm.lon, mhm.lat, mhm.mean_tt_opt)