# -*- coding: utf-8 -*- """ util.py -- Utilities for working with PRMS data or other functionality that aren't appropriate to put elsewhere at this time. """ import warnings import os, shutil, json import numpy as np import pandas as pd def calc_emp_CDF(data): # changed function name for PEP 8 style warnings.warn("calc_emp_CDF is deprecated, please use "+\ "util.calc_emp_cdf instead", DeprecationWarning) return calc_emp_cdf(data) def calc_emp_cdf(data): """ Create empirical CDF of arbitrary data Arguments: data (array_like) : array to calculate CDF on Returns: X (numpy.ndarray) : array of x values of CDF (sorted data) F (numpy.ndarray) : array of CDF values for each X value or cumulative exceedence probability, in [0,1]. """ n_bins = len(data) X = np.sort(data) F = np.array(range(n_bins))/float(n_bins) return X,F def Kolmogorov_Smirnov(uncond, cond, n_bins=10000): # changed function name for PEP 8 style warnings.warn("Kolmogorov_Smirnov is deprecated, please use "+\ "util.komogorov_smirnov instead", DeprecationWarning) return kolmogorov_smirnov(uncond, cond, n_bins=10000) def kolmogorov_smirnov(uncond, cond, n_bins=10000): """ Calculate the Kolmogorov-Smirnov statistic between two datasets by first computing their empirical CDFs Arguments: uncond (array_like) : data for creating the unconditional CDF. cond (array_like) : data for creating the conditional CDF n_bins (int) : number of bins for both CDFs, note if n_bins > length of either dataset then CDF values are interpolated by numpy Returns: KS (float) : Kolmogorov-Smirnov statistic, i.e. absolute max distance between uncond and cond CDFs """ # create unconditional CDF (F_Uc) H,X = np.histogram(uncond, bins=n_bins, normed=True) dx = X[1] - X[0] F_Uc = np.cumsum(H)*dx # create conditional CDF (F_C) H,X = np.histogram(cond, bins=n_bins, normed=True) dx = X[1] - X[0] F_C = np.cumsum(H)*dx # Calc max absolulte divergence KS = np.max(np.abs(F_Uc - F_C)) return KS def remove_all_optimization_sims_of_other_stage(work_directory, stage): """ Track number of simulation directories not tracked by a specific stage and recursively delete them and their contents. This was created to avoid having nutracked simulations in an optimizer working directory for example when an optimization method was interupted before data was saved to a meta data file. Arguments: work_directory (str) : Directory to look for Optimization metadata json files and simulation directories to keep or remove. stage (str) : Optimization stage that will not have its simulation data deleted. All other stages if any are found in metadata files will have their associated simulation directories deleted. Returns: None """ from .optimizer import OptimizationResult # avoid circular import try: result = OptimizationResult(work_directory,stage=stage) tracked_dirs = [] for f in result.metadata_json_paths[stage]: with open(f) as fh: json_data = json.load(fh) tracked_dirs.extend(json_data.get('sim_dirs')) count = 0 for d in os.listdir(result.working_dir): path = os.path.join(result.working_dir, d) if path in tracked_dirs: continue elif os.path.isdir(path) and '_archived' not in path: count+=1 for dirpath, dirnames, filenames in os.walk(path,\ topdown=False): shutil.rmtree(dirpath, ignore_errors=True) # if no json file in working dir for given stage, delete any other sim dirs except: count = 0 for d in os.listdir(work_directory): path = os.path.join(work_directory, d) if os.path.isdir(path) and '_archived' not in path: count+=1 for dirpath, dirnames, filenames in os.walk(path,\ topdown=False): shutil.rmtree(dirpath, ignore_errors=True) print('deleted {} simulations that were either not tracked by a JSON file'\ .format(count) + ' or were not part of {} optimization stage'\ .format(stage)) def delete_files(work_directory, file_name=''): """ Recursively delete all files of a certain name from multiple PRMS simulations that are within a given directory. Can be useful to removw large files that are no longer needed. For example initial condition output files are often large and not always used, similarly animation, data, control, ... files may no longer be needed. Arguments: work_directory (str) : path to directory with simulations. file_name (str) : Name of the PRMS input or output file(s) to be removed, default = '' empty string- nothing will be deleted. Example: e.g. if you have several simulation directories: >>> "test/results/intcp:-26.50_slope:0.49", "test/results/intcp:-11.68_slope:0.54", "test/results/intcp:-4.70_slope:0.51", "test/results/intcp:-35.39_slope:0.39", "test/results/intcp:-20.91_slope:0.41" each of these contains an '/inputs' folder with a duplicate data file that you would like to delete. In this case, delete all data files like so: >>> work_dir = 'test/results/' >>> delete_files(work_dir, file_name='data') Returns: None """ for dirpath, dirnames, filenames in os.walk(work_directory, topdown=False): paths = (os.path.join(dirpath, filename) for filename in filenames\ if filename == file_name) for path in paths: os.remove(path) def load_statvar(statvar_file): """ Read the statvar file and load into a datetime indexed Pandas dataframe object Arguments: statvar_file (str): statvar file path Returns: (pandas.DataFrame) Pandas DataFrame of PRMS variables date indexed from statvar file """ # make list of statistical output variables for df header column_list = ['index', 'year', 'month', 'day', 'hh', 'mm', 'sec'] # append to header list the variables present in the file with open(statvar_file, 'r') as inf: for idx, l in enumerate(inf): # first line is always number of stat variables if idx == 0: n_statvars = int(l) elif idx <= n_statvars and idx != 0: column_list.append(l.rstrip().replace(' ', '_')) else: break # arguments for read_csv function missing_value = -999 skiprows = n_statvars+1 df = pd.read_csv( statvar_file, delim_whitespace=True, skiprows=skiprows, header=None, na_values=[missing_value] ) # apply correct header names using metadata retrieved from file df.columns = column_list date = pd.Series( pd.to_datetime(df.year*10000+df.month*100+df.day, format='%Y%m%d'), index=df.index ) # make the df index the datetime for the time series data df.index = pd.to_datetime(date) # drop unneeded columns df.drop(['index', 'year', 'month', 'day', 'hh', 'mm', 'sec'], axis=1, inplace=True) # name dataframe axes (index,columns) df.columns.name = 'statistical_variables' df.index.name = 'date' return df def load_data_file(data_file): # changed function name for PEP 8 style warnings.warn("load_data_file is deprecated, please use "+\ "util.load_data instead", DeprecationWarning) return load_data(data_file) def load_data(data_file): """ Read the data file and load into a datetime indexed Pandas dataframe object. Arguments: data_file (str): data file path Returns: df (pandas.DataFrame): Pandas dataframe of input time series data from data file with datetime index """ # valid input time series that can be put into a data file valid_input_variables = ('gate_ht', 'humidity', 'lake_elev', 'pan_evap', 'precip', 'rain_day', 'runoff', 'snowdepth', 'solrad', 'tmax', 'tmin', 'wind_speed') # starting list of names for header in dataframe column_list = ['year', 'month', 'day', 'hh', 'mm', 'sec'] # append to header list the variables present in the file with open(data_file, 'r') as inf: for idx, l in enumerate(inf): # first line always string identifier of the file- may use later if idx == 0: data_head = l.rstrip() elif l.startswith('/'): # comment lines continue # header lines with name and number of input variables if l.startswith(valid_input_variables): # split line into list, first element name and # second number of columns h = l.split() # more than one input time series of a particular variable if int(h[1]) > 1: for el in range(int(h[1])): tmp = '{var_name}_{var_ind}'.format(var_name=h[0], var_ind=el+1) column_list.append(tmp) elif int(h[1]) == 1: column_list.append(h[0]) # end of header info and begin time series input data if l.startswith('#'): skip_line = idx+1 break # read data file into pandas dataframe object with correct header names missing_value = -999 # missing data representation df = pd.read_csv(data_file, header=-1, skiprows=skip_line, delim_whitespace=True, na_values=[missing_value]) # apply correct header names using metadata retrieved from file df.columns = column_list # create date column date = pd.Series( pd.to_datetime(df.year*10000+df.month*100+df.day, format='%Y%m%d'), index=df.index ) df.index = pd.to_datetime(date) # make the df index the datetime # drop unneeded columns df.drop(['year', 'month', 'day', 'hh', 'mm', 'sec'], axis=1, inplace=True) df.columns.name = 'input variables' df.index.name = 'date' # name dataframe axes (index,columns) return df def nash_sutcliffe(observed, modeled): """ Calculates the Nash-Sutcliffe Goodness-of-fit Arguments: observed (numpy.ndarray): historic observational data modeled (numpy.ndarray): model output with matching time index """ numerator = sum((observed - modeled)**2) denominator = sum((observed - np.mean(observed))**2) return 1 - (numerator/denominator) def percent_bias(observed, modeled): """ Calculates percent bias Arguments: observed (numpy.ndarray): historic observational data modeled (numpy.ndarray): model output with matching time index """ return 100 * ( sum( modeled - observed ) / sum( observed ) ) def rmse(observed, modeled): """ Calculates root mean squared error Arguments: observed (numpy.ndarray): historic observational data modeled (numpy.ndarray): model output with matching time index """ return np.sqrt( sum((observed - modeled)**2) / len(observed) )