From 2023d8beaab4a6a7a572e48177b0bd3a26147482 Mon Sep 17 00:00:00 2001 From: navass11 Date: Thu, 16 Jul 2020 12:32:14 +0200 Subject: [PATCH] Upgrade to python 3.6 Various lines of code have been updated to work correctly in python 3.6 --- docs/user-guide/parameters.md | 2 +- .../parameters-checkpoint.py | 654 ++++++++++++++++ .../.ipynb_checkpoints/history-checkpoint.py | 725 ++++++++++++++++++ .../param_file-checkpoint.py | 469 +++++++++++ .../.ipynb_checkpoints/share-checkpoint.py | 345 +++++++++ rvic/core/history.py | 2 +- rvic/core/param_file.py | 31 +- rvic/core/share.py | 12 +- rvic/parameters.py | 2 +- rvic/version.py | 4 +- samples/configs/rvic.convolution.rasm.cfg | 2 +- ...columbia_sample_pour_points-checkpoint.csv | 3 + .../rasm_sample_pour_points-checkpoint.csv | 3 + .../columbia_sample_pour_points.csv | 4 +- .../pour_points/rasm_sample_pour_points.csv | 4 +- scripts/.ipynb_checkpoints/rvic-checkpoint | 62 ++ 16 files changed, 2292 insertions(+), 32 deletions(-) create mode 100644 rvic/.ipynb_checkpoints/parameters-checkpoint.py create mode 100644 rvic/core/.ipynb_checkpoints/history-checkpoint.py create mode 100644 rvic/core/.ipynb_checkpoints/param_file-checkpoint.py create mode 100644 rvic/core/.ipynb_checkpoints/share-checkpoint.py create mode 100644 samples/pour_points/.ipynb_checkpoints/columbia_sample_pour_points-checkpoint.csv create mode 100644 samples/pour_points/.ipynb_checkpoints/rasm_sample_pour_points-checkpoint.csv create mode 100644 scripts/.ipynb_checkpoints/rvic-checkpoint diff --git a/docs/user-guide/parameters.md b/docs/user-guide/parameters.md index bb0ad64..ef8526d 100644 --- a/docs/user-guide/parameters.md +++ b/docs/user-guide/parameters.md @@ -146,7 +146,7 @@ A csv file that describes the routing of flow to the edge of the origin grid cel - Type: bool - Valid values: True, False - *Note: In most cases, this should be set to True unless routing is being done to nested basins, in which case, this should be set to False.* + Note: *True when routing to coastal grid cells, else False.* ###POUR_POINTS 1. **FILE_NAME** diff --git a/rvic/.ipynb_checkpoints/parameters-checkpoint.py b/rvic/.ipynb_checkpoints/parameters-checkpoint.py new file mode 100644 index 0000000..e98c7aa --- /dev/null +++ b/rvic/.ipynb_checkpoints/parameters-checkpoint.py @@ -0,0 +1,654 @@ +# -*- coding: utf-8 -*- +''' +RVIC parameter file development driver +''' +import os +import numpy as np +import pandas as pd +from logging import getLogger +from .core.log import init_logger, close_logger, LOG_NAME +from .core.multi_proc import error +from .core.utilities import make_directories, copy_inputs, strip_invalid_char +from .core.utilities import read_netcdf, tar_inputs, latlon2yx +from .core.utilities import check_ncvars, clean_file, read_domain +from .core.utilities import search_for_channel +from .core.aggregate import make_agg_pairs, aggregate +from .core.make_uh import rout +from .core.share import NcGlobals +from .core.write import write_agg_netcdf +from .core.variables import Point +from .core.param_file import finish_params +from .core.config import read_config +from .core.pycompat import OrderedDict, iteritems, pyrange, basestring + +try: + from .core.remap import remap + remap_available = True +except ImportError: + remap_available = False + +# global multiprocessing results container +results = {} + + +# -------------------------------------------------------------------- # +# Top level driver +def parameters(config, numofproc=1): + ''' + Top level function for RVIC parameter generation function. + + Parameters + ---------- + config : str or dict + Path to RVIC parameters configuration file or dictionary of + configuration options. + numofproc : int + Number of processors to use when developing RVIC parameters. + ''' + + # ---------------------------------------------------------------- # + # Initilize + uh_box, fdr_data, fdr_vatts, dom_data, \ + outlets, config_dict, directories = gen_uh_init(config) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Get main logger + log = getLogger(LOG_NAME) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Run + if numofproc > 1: + from multiprocessing import Pool + pool = Pool(processes=numofproc) + status = [] + + for i, (key, outlet) in enumerate(iteritems(outlets)): + log.info('On Outlet #%s of %s', i + 1, len(outlets)) + stat = pool.apply_async(gen_uh_run, + (uh_box, fdr_data, fdr_vatts, + dom_data, outlet, config_dict, + directories), + callback=store_result, + error_callback=error) + # Store the result + status.append(stat) + + # Close the pool + pool.close() + + # Check that everything worked + [stat.get() for stat in status] + + pool.join() + + outlets = OrderedDict(sorted(results.items(), reverse=True)) + else: + for i, (key, outlet) in enumerate(iteritems(outlets)): + outlet = gen_uh_run(uh_box, fdr_data, fdr_vatts, dom_data, outlet, + config_dict, directories) + + if not outlets: + raise ValueError('outlets in parameters are empty') + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Finally, make the parameter file + gen_uh_final(outlets, dom_data, config_dict, directories) + # ---------------------------------------------------------------- # + return +# -------------------------------------------------------------------- # + + +def gen_uh_init(config): + '''Initialize RVIC parameters script. + + This function: + - Reads the configuration file + - Sets up the RVIC case directories + - Copies all input files to the case directory + - Initializes the logging + - Reads the pour-points, uh-box, FDR, and domain files + - Aggregates pour points into outlet grid cells + + Parameters + ---------- + config : str or dict + Path to RVIC parameters configuration file or dictionary of + configuration options. + + Returns + ---------- + uh_box : numpy.ndarray + UH-box array. + fdr_data : dict + Dictionary of arrays of flow direction, velocity, diffusion, etc. + This dictionary includes all the variables from the FDR netCDF file. + fdr_vatts : dict + Dictionary of attributes from the FDR netCDF file. + dom_data : dict + Dictionary of arrays of mask, fraction, lats, lons, etc. + This dictionary includes all the variables from the domain netCDF file. + outlets : dict + Dictionary of outlet `Point` objects. + config_dict : dict + Dictionary of values from the configuration file. + directories : dict + Dictionary of directories created by this function. + ''' + + # ---------------------------------------------------------------- # + # Read Configuration files + if isinstance(config, dict): + config_dict = config + else: + config_dict = read_config(config) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Import optional modules + if config_dict['OPTIONS']['REMAP'] and not remap_available: + raise ValueError('Problem importing remap module ' + 'check to make sure cdo.py is available)') + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Setup Directory Structures + directories = make_directories(config_dict['OPTIONS']['CASE_DIR'], + ['plots', 'logs', 'params', 'inputs']) + directories.update(make_directories(config_dict['OPTIONS']['TEMP_DIR'], + ['aggregated', 'remapped'])) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # copy inputs to $case_dir/inputs and update configuration + if not isinstance(config, dict): + config_dict = copy_inputs(config, directories['inputs']) + options = config_dict['OPTIONS'] + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Start Logging + log = init_logger(directories['logs'], options['LOG_LEVEL'], + options['VERBOSE']) + + for direc in directories: + log.info('%s directory is %s', direc, directories[direc]) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Read Pour Points files + try: + pour_points = pd.read_csv(config_dict['POUR_POINTS']['FILE_NAME'], + comment='#') + log.info('Opened Pour Points File: %s', + config_dict['POUR_POINTS']['FILE_NAME']) + if not (all(x in list(pour_points.keys()) for x in ['lons', 'lats']) or + all(x in list(pour_points.keys()) for x in ['x', 'y'])): + raise ValueError('Pour Points File must include ' + 'variables (lons, lats) or (x, y)') + if 'names' in pour_points: + pour_points.fillna(inplace=True, value='unknown') + for i, name in enumerate(pour_points.names): + pour_points.ix[i, 'names'] = strip_invalid_char(name) + + pour_points.drop_duplicates(inplace=True) + pour_points.dropna() + except Exception as e: + log.error('Error opening pour points file: %s', + config_dict['POUR_POINTS']['FILE_NAME']) + log.exception(e) + raise + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Read uh box file + uh_file = config_dict['UH_BOX']['FILE_NAME'] + uh_header = int(config_dict['UH_BOX']['HEADER_LINES']) + uh_box = {} + try: + uh_box['time'], uh_box['func'] = np.genfromtxt(uh_file, + skip_header=uh_header, + delimiter=',', + unpack=True) + log.info('Opened UHbox File: %s', + config_dict['UH_BOX']['FILE_NAME']) + except: + log.exception('Error opening uh_box file: %s', + config_dict['POUR_POINTS']['FILE_NAME']) + raise + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Read FDR file + fdr_file = config_dict['ROUTING']['FILE_NAME'] + fdr_var = config_dict['ROUTING']['FLOW_DIRECTION_VAR'] + fdr_lat = config_dict['ROUTING']['LATITUDE_VAR'] + fdr_lon = config_dict['ROUTING']['LONGITUDE_VAR'] + fdr_vel = config_dict['ROUTING']['VELOCITY'] + fdr_dif = config_dict['ROUTING']['DIFFUSION'] + try: + fdr_data, fdr_vatts, _ = read_netcdf(fdr_file) + fdr_shape = fdr_data[fdr_var].shape + + # ---------------------------------------------------------------- # + # Check latitude order, flip if necessary. + if fdr_data[fdr_lat][-1] > fdr_data[fdr_lat][0]: + log.debug('Flow Direction inputs came in upside down, flipping ' + 'everything now.') + + remove_vars = [] + + for var, data in iteritems(fdr_data): + log.debug('flipping %s', var) + if data.ndim >= 1 and var != fdr_lon: + fdr_data[var] = np.flipud(data) + elif data.ndim == 0: + remove_vars.append(var) + + if remove_vars: + for var in remove_vars: + del fdr_data[var] + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Add velocity and/or diffusion grids if not present yet + if not isinstance(fdr_vel, basestring): + fdr_data['velocity'] = \ + np.zeros(fdr_shape, dtype=np.float64) + fdr_vel + config_dict['ROUTING']['VELOCITY'] = 'velocity' + log.info('Added velocity grid to fdr_data') + if not isinstance(fdr_dif, basestring): + fdr_data['diffusion'] = \ + np.zeros(fdr_shape, dtype=np.float64) + fdr_dif + config_dict['ROUTING']['DIFFUSION'] = 'diffusion' + log.info('Added diffusion grid to fdr_data') + if ('SOURCE_AREA_VAR' not in config_dict['ROUTING'] or + config_dict['ROUTING']['SOURCE_AREA_VAR'] not in fdr_data): + log.warning('Upstream `SOURCE_AREA` was not provided, output ' + 'source area will be zero.') + config_dict['ROUTING']['SOURCE_AREA_VAR'] = 'src_area' + fdr_data[config_dict['ROUTING']['SOURCE_AREA_VAR']] = \ + np.zeros(fdr_shape, dtype=np.float64) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + fdr_data['resolution'] = np.abs(fdr_data[fdr_lon][1] - + fdr_data[fdr_lon][0]) + check_ncvars(config_dict['ROUTING'], list(fdr_data.keys())) + # ---------------------------------------------------------------- # + + log.info('Opened FDR File: %s', fdr_file) + except: + log.exception('Error opening FDR file') + raise + # ---------------------------------------------------------------- # + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Read domain file + domain = config_dict['DOMAIN'] + dom_data = read_domain(domain)[0] + log.info('Opened Domain File: %s', domain['FILE_NAME']) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # If remap is False, domain coordinates needs to be in the fdr coordinates + # We can move the unit hydrographs to the domain grid later + if options['AGGREGATE'] and not options['REMAP']: + log.error('RVIC parameter generation requires REMAP option to be True' + ' if AGGREGATE is True') + raise ValueError('Invalid option') + + # If remap is False, then the resolution needs to match the routing data + if not options['REMAP']: + domain_res = np.abs(dom_data[domain['LONGITUDE_VAR']][0, 1] - + dom_data[domain['LONGITUDE_VAR']][0, 0]) + if not np.isclose(fdr_data['resolution'], domain_res): + log.error('routing grid resolution: %s', fdr_data['resolution']) + log.error('domain grid resolution: %s', domain_res) + raise ValueError('If remap is false, domain and routing grid ' + 'resolutions must match.') + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Group pour points (if aggregate) + if options['AGGREGATE']: + outlets = make_agg_pairs(pour_points, dom_data, fdr_data, config_dict) + + log.info('Finished making agg pairs of ' + 'pour points and outlet grid cells') + + else: + outlets = OrderedDict() + if all(x in list(pour_points.keys()) for x in ['x', 'y', + 'lons', 'lats']): + lats = pour_points['lats'].values + lons = pour_points['lons'].values + routys = pour_points['y'].values + routxs = pour_points['x'].values + elif all(x in list(pour_points.keys()) for x in ['x', 'y']): + # use x and y (assume from routing inputs grid) + # find lons and lats from xs and ys + routys = pour_points['y'].values + routxs = pour_points['x'].values + lats = fdr_data[fdr_lat][routys] + lons = fdr_data[fdr_lon][routxs] + else: + # use and lats to find xs and ys + lats = pour_points['lats'].values + lons = pour_points['lons'].values + + # find x and y on routing grid + routys, routxs = latlon2yx(plats=lats, + plons=lons, + glats=fdr_data[fdr_lat], + glons=fdr_data[fdr_lon]) + + if options['SEARCH_FOR_CHANNEL']: + routys, routxs = search_for_channel( + fdr_data[config_dict['ROUTING']['SOURCE_AREA_VAR']], + routys, routxs, tol=10, search=5) + + # update lats and lons + lats = fdr_data[fdr_lat][routys] + lons = fdr_data[fdr_lon][routxs] + + # Find location on domain grid + domys, domxs = latlon2yx(plats=lats, + plons=lons, + glats=dom_data[domain['LATITUDE_VAR']], + glons=dom_data[domain['LONGITUDE_VAR']]) + + for i in pyrange(len(lats)): + if 'names' in list(pour_points.keys()): + name = pour_points['names'].values[i] + name = name.replace("'", '').replace(' ', '_') + else: + # fill name filed with p-outlet_num + name = 'p-{0}'.format(i) + + outlets[i] = Point(lat=lats[i], + lon=lons[i], + domx=domxs[i], + domy=domys[i], + routx=routxs[i], + routy=routys[i], + name=name, + cell_id=dom_data['cell_ids'][domys[i], domxs[i]]) + + outlets[i].pour_points = [outlets[i]] + # ---------------------------------------------------------------- # + + log.debug(outlets) + log.info('Finished with gen_uh_init') + log.info('-------------------------------------------------------------\n') + + return (uh_box, fdr_data, fdr_vatts, dom_data, outlets, + config_dict, directories) +# -------------------------------------------------------------------- # + + +def gen_uh_run(uh_box, fdr_data, fdr_vatts, dom_data, outlet, config_dict, + directories): + ''' + Develop unit hydrographs for one outlet `Point`. + + Parameters + ---------- + uh_box : numpy.ndarray + UH-box array. + fdr_data : dict + Dictionary of arrays of flow direction, velocity, diffusion, etc. + This dictionary includes all the variables from the FDR netCDF file. + fdr_vatts : dict + Dictionary of attributes from the FDR netCDF file. + dom_data : dict + Dictionary of arrays of mask, fraction, lats, lons, etc. + This dictionary includes all the variables from the domain netCDF file. + outlet : Point + Outlet Point + config_dict : dict + Dictionary of values from the configuration file. + directories : dict + Dictionary of directories created by gen_uh_init. + + Returns + ---------- + outlet: Point + Point object with unit hydrographs. + ''' + log = getLogger(LOG_NAME) + + log.info('Running outlet cell id %s', outlet.cell_id) + + agg_data = {} + domain = config_dict['DOMAIN'] + dom_lat = domain['LATITUDE_VAR'] + dom_lon = domain['LONGITUDE_VAR'] + dom_mask = domain['LAND_MASK_VAR'] + options = config_dict['OPTIONS'] + + # ------------------------------------------------------------ # + # netCDF variable options + ncvaropts = {} + if 'NETCDF_ZLIB' in options: + ncvaropts['zlib'] = options['NETCDF_ZLIB'] + if 'NETCDF_COMPLEVEL' in options: + ncvaropts['complevel'] = options['NETCDF_COMPLEVEL'] + if 'NETCDF_SIGFIGS' in options: + ncvaropts['least_significant_digit'] = options['NETCDF_SIGFIGS'] + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Loop over pour points + n_pour_points = len(outlet.pour_points) + for j, pour_point in enumerate(outlet.pour_points): + + log.info('On pour_point #%s of %s', j + 1, n_pour_points) + + # -------------------------------------------------------- # + # Make the Unit Hydrograph Grid + rout_data = rout(pour_point, uh_box, fdr_data, fdr_vatts, + config_dict['ROUTING']) + + log.debug('Done routing to pour_point') + log.debug('rout_data: %s, %s', rout_data['unit_hydrograph'].min(), + rout_data['unit_hydrograph'].max()) + log.debug('rout_data sum: %s, %s', rout_data['unit_hydrograph'].sum(), + rout_data['fraction'].sum()) + + # -------------------------------------------------------- # + + # -------------------------------------------------------- # + # aggregate + if options['AGGREGATE']: + if j != len(outlet.pour_points) - 1: + agg_data = aggregate(rout_data, agg_data, + res=fdr_data['resolution']) + else: + agg_data = aggregate(rout_data, agg_data, + res=fdr_data['resolution'], + pad=options['AGG_PAD'], + maskandnorm=True) + + log.debug('agg_data: %s, %s', + agg_data['unit_hydrograph'].min(), + agg_data['unit_hydrograph'].max()) + else: + agg_data = rout_data + # -------------------------------------------------------- # + + # ------------------------------------------------------------ # + # write temporary file #1 + if options['REMAP']: + glob_atts = NcGlobals( + title='RVIC Unit Hydrograph Grid File', + RvicPourPointsFile=os.path.split( + config_dict['POUR_POINTS']['FILE_NAME'])[1], + RvicUHFile=os.path.split(config_dict['UH_BOX']['FILE_NAME'])[1], + RvicFdrFile=os.path.split(config_dict['ROUTING']['FILE_NAME'])[1], + RvicDomainFile=os.path.split(domain['FILE_NAME'])[1]) + + temp_file_1 = os.path.join( + directories['aggregated'], + 'aggUH_{0}.nc'.format(outlet.name.replace(' ', '_'))) + + write_agg_netcdf(temp_file_1, agg_data, glob_atts, + options['NETCDF_FORMAT'], **ncvaropts) + + # -------------------------------------------------------- # + # Remap temporary file #1 to temporary file #2 + temp_file_2 = os.path.join( + directories['remapped'], + 'remapUH_{0}.nc'.format(outlet.name.replace(' ', '_'))) + + remap(domain['FILE_NAME'], temp_file_1, temp_file_2) + + # -------------------------------------------------------- # + # Read temporary file #2 + final_data = read_netcdf( + temp_file_2, variables=['unit_hydrograph', 'fraction', dom_lat])[0] + + # -------------------------------------------------------- # + # Check latitude order, flip if necessary. + if final_data[dom_lat].ndim == 1: + if final_data[dom_lat][-1] > final_data[dom_lat][0]: + var_list = list(final_data.keys()) + + log.debug('Remapped inputs came in upside down, flipping %s', + ', '.join(var_list)) + # flip lattiutude and fraction along y axis (axis 0) + final_data[dom_lat] = final_data[dom_lat][::-1] + final_data['fraction'] = final_data['fraction'][::-1, :] + # flip unit hydrograph along y axis (axis 1) + final_data['unit_hydrograph'] = \ + final_data['unit_hydrograph'][:, ::-1, :] + assert dom_data['cord_lats'][0] == final_data[dom_lat][0] + # -------------------------------------------------------- # + + # -------------------------------------------------------- # + # Clean temporary file #2 (if applicable) + if config_dict['OPTIONS']['CLEAN']: + clean_file(temp_file_1) + clean_file(temp_file_2) + # -------------------------------------------------------- # + + else: + # -------------------------------------------------------- # + # Put the agg data back onto the original grid + uh_shape = (agg_data['unit_hydrograph'].shape[0], ) + \ + dom_data[dom_mask].shape + final_data = {} + final_data['unit_hydrograph'] = np.zeros(uh_shape, dtype=np.float64) + final_data['fraction'] = np.zeros(dom_data[dom_mask].shape, + dtype=np.float64) + + bys, bxs = np.nonzero(agg_data['fraction']) + + ys = bys + outlet.domy - agg_data['basiny'] + xs = bxs + outlet.domx - agg_data['basinx'] + + if (ys < 0).any() or (xs < 0).any(): + raise ValueError('Negative indicies found when mapping ' + '`non-remapped` rout_data to domain grid.') + + final_data['unit_hydrograph'][:, ys, xs] = \ + agg_data['unit_hydrograph'][:, bys, bxs] + final_data['fraction'][ys, xs] = agg_data['fraction'][bys, bxs] + # -------------------------------------------------------- # + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Add to 'adjust fractions structure' + y, x = np.nonzero((final_data['fraction'] > 0.) * + (dom_data[dom_mask] > np.finfo(np.float).resolution)) + + # From final data + outlet.time = np.arange(final_data['unit_hydrograph'].shape[0]) + outlet.fractions = final_data['fraction'][y, x] + outlet.unit_hydrograph = final_data['unit_hydrograph'][:, y, x] + + # From domain data + outlet.lon_source = dom_data[dom_lon][y, x] + outlet.lat_source = dom_data[dom_lat][y, x] + outlet.cell_id_source = dom_data['cell_ids'][y, x] + outlet.x_source = x + outlet.y_source = y + + # Make sure the inds are all greater than zero, ref: Github #79 + assert all(outlet.cell_id_source >= 0) + assert all(outlet.x_source >= 0) + assert all(outlet.y_source >= 0) + # ---------------------------------------------------------------- # + return outlet +# -------------------------------------------------------------------- # + + +def gen_uh_final(outlets, dom_data, config_dict, directories): + ''' + Make the RVIC Parameter File + + Parameters + ---------- + outlets : dict + Dictionary of outlet `Point` objects. + dom_data : dict + Dictionary of arrays of mask, fraction, lats, lons, etc. + This dictionary includes all the variables from the domain netCDF file. + config_dict : dict + Dictionary of values from the configuration file. + directories : dict + Dictionary of directories created by gen_uh_init. + ''' + log = getLogger(LOG_NAME) + + log.info('In gen_uh_final') + + if not len(outlets) > 0: + raise ValueError('outlets in gen_uh_final are empty') + + # ---------------------------------------------------------------- # + # Write the parameter file + param_file, today = finish_params(outlets, dom_data, config_dict, + directories) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # tar the inputs directory / log file + inputs_tar = tar_inputs(directories['inputs'], suffix=today) + log_tar = tar_inputs(log.filename) + + log.info('Done with RvicGenParam.') + log.info('Location of Inputs: %s', inputs_tar) + log.info('Location of Log: %s', log_tar) + log.info('Location of Parmeter File %s', param_file) + + close_logger() + # ---------------------------------------------------------------- # + return +# -------------------------------------------------------------------- # + + +# -------------------------------------------------------------------- # +# store_result helper function +def store_result(result): + ''' + Store values returned by a multiprocessing.pool member. + + This is called whenever foo_pool(i) returns a result. + result_list is modified only by the main process, not the pool workers. + + Parameters + ---------- + result : object + Result to append to the global `results` list + + Globals + ---------- + results : dict + Global results container for multiprocessing results to be appended to. + ''' + results[result.cell_id] = result +# -------------------------------------------------------------------- # diff --git a/rvic/core/.ipynb_checkpoints/history-checkpoint.py b/rvic/core/.ipynb_checkpoints/history-checkpoint.py new file mode 100644 index 0000000..df11bec --- /dev/null +++ b/rvic/core/.ipynb_checkpoints/history-checkpoint.py @@ -0,0 +1,725 @@ +# -*- coding: utf-8 -*- +''' +history.py + +Summary: + This is the core history file module for the rvic model. + The core of the module is the Tape class. The basic procedure is as + follows: + - initialization: sets tape options, determines filenames, etc. + - update: method that incorporates new fluxes into the history tape. + - __next_update_out_data: method to determine when to update the + out_data container + - __next_write_out_data: method to determine when to write the out_data + container + - finish: method to close all remaining history tapes. +''' + +import os +import numpy as np +from netCDF4 import Dataset, date2num, num2date, stringtochar +from datetime import datetime +from .time_utility import ord_to_datetime +from logging import getLogger +from .log import LOG_NAME +from .share import SECSPERDAY, HOURSPERDAY, TIMEUNITS +from .share import NC_INT, NC_FLOAT, NC_CHAR +from .share import NC_DOUBLE, WATERDENSITY, MONTHSPERYEAR +from .pycompat import iteritems +from . import share + + +# -------------------------------------------------------------------- # +# create logger +log = getLogger(LOG_NAME) +# -------------------------------------------------------------------- # + + +# -------------------------------------------------------------------- # +# RVIC History File Object +class Tape(object): + ''' History Tape Object''' + + # ---------------------------------------------------------------- # + # Init + def __init__(self, time_ord, caseid, rvar, tape_num=0, + fincl=('streamflow'), mfilt=1, ndens=2, nhtfrq=0, + avgflag='A', units='kg m-2 s-1', + file_format='NETCDF4_CLASSIC', outtype='grid', + grid_lons=False, grid_lats=False, grid_area=None, out_dir='.', + calendar='standard', glob_ats=None, zlib=True, complevel=4, + least_significant_digit=None): + self._tape_num = tape_num + self._time_ord = time_ord # Days since basetime + self._caseid = caseid # Case ID and prefix for outfiles + self._fincl = list(fincl) # Fields to include in history file + self._mfilt = mfilt # Maximum number of time samples + self._ndens = ndens + if self._ndens == 1: # Output file precision + self._ncprec = NC_FLOAT + else: + self._ncprec = NC_DOUBLE + self._nhtfrq = nhtfrq # Write frequency + self._avgflag = avgflag # Average Flag (A,I,X,M) + self._outtype = outtype # Outfile type (grid, array) + self._count = 0 + self.files_count = 0 + self._file_format = file_format + self._calendar = calendar + self._out_dir = out_dir + self._glob_ats = glob_ats + + self.__get_rvar(rvar) # Get the initial rvar fields + self._grid_shape = grid_area.shape + self._out_data = {} + + # ------------------------------------------------------------ # + # calculate the step size for each out_data timestep (units=days) + if self._nhtfrq > 0: + # If some number of timesteps + self._out_data_stepsize = self._nhtfrq * self._dt / SECSPERDAY + elif self._nhtfrq < 0: + # If some number hours + self._out_data_stepsize = -1 * self._nhtfrq / HOURSPERDAY + else: + # If monthly + self._out_data_stepsize = None # varies by month + log.debug('_out_data_stepsize: %s', self._out_data_stepsize) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Get Grid Lons/Lats if outtype is grid + if outtype.lower() == 'grid': + self._out_data_shape = self._grid_shape + if type(grid_lons) == np.ndarray and type(grid_lats) == np.ndarray: + self._grid_lons = grid_lons + self._grid_lats = grid_lats + else: + raise ValueError('Must include grid lons / lats if ' + 'outtype == grid') + elif outtype.lower() == 'array': + self._out_data_shape = (self._num_outlets, ) + else: + raise ValueError('Unknown value for outtype: {0}'.format(outtype)) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Get units multiplier (size of noutlets) + self._units = units + if units in ['kg/m2/s', 'kg m-2 s-1', 'kg m^-2 s^-1', + 'kg*m-2*s-1', 'kg s-1 m-2']: + self._units_mult = np.ones_like(self._outlet_y_ind, + dtype=np.float64) + elif units in ['m3/s', 'm^3/s', 'm3 s-1']: + # kg/m2/s --> m3/s + self._units_mult = grid_area[self._outlet_y_ind, + self._outlet_x_ind] + self._units_mult /= WATERDENSITY + elif units in ['mm/day', 'mm d-1', 'mm d^-1', 'mm/day']: + # kg/m2/s --> mm/day over basin area + self._units_mult = grid_area[self._outlet_y_ind, + self._outlet_x_ind] + self._units_mult *= SECSPERDAY + self._units_mult /= WATERDENSITY + self._units_mult /= self._outlet_upstream_area + elif units in ['gal/day', 'gpd', 'gal d-1']: + self._units_mult = grid_area[self._outlet_y_ind, + self._outlet_x_ind] + self._units_mult /= WATERDENSITY + self._units_mult *= 2.28E7 + elif units in ['cfs', 'ft^3 s-1', 'f3/s']: + self._units_mult = grid_area[self._outlet_y_ind, + self._outlet_x_ind] + self._units_mult /= WATERDENSITY + self._units_mult *= 35.3 + elif units in ['acre-ft/d']: + self._units_mult = grid_area[self._outlet_y_ind, + self._outlet_x_ind] + self._units_mult /= WATERDENSITY + self._units_mult *= 70.0 + else: + raise ValueError('{0} is not a valid units string'.format(units)) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # netCDF variable options + self.ncvaropts = {'zlib': zlib, + 'complevel': complevel, + 'least_significant_digit': least_significant_digit} + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # get current timestamp + self._timestamp = ord_to_datetime(self._time_ord, TIMEUNITS, + self._calendar) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Initialize the temporary history fields + self._temp_data = {} + for field in self._fincl: + self._temp_data[field] = np.zeros(self._num_outlets, + dtype=np.float64) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Determine the format of the output filename + if self._avgflag == 'I': + self._fname_format = os.path.join( + out_dir, '%s.rvic.h%s%s.%%Y-%%m-%%d-%%H-%%M-%%S.nc' % + (self._caseid, self._tape_num, self._avgflag.lower())) + else: + if self._nhtfrq == 0: + self._fname_format = os.path.join( + out_dir, + '%s.rvic.h%s%s.%%Y-%%m.nc' % + (self._caseid, self._tape_num, self._avgflag.lower())) + elif (self._nhtfrq == -24) or (nhtfrq * self._dt == SECSPERDAY): + self._fname_format = os.path.join( + out_dir, + '%s.rvic.h%s%s.%%Y-%%m-%%d.nc' % + (self._caseid, self._tape_num, self._avgflag.lower())) + else: + self._fname_format = os.path.join( + out_dir, + '%s.rvic.h%s%s.%%Y-%%m-%%d-%%H.nc' % + (self._caseid, self._tape_num, self._avgflag.lower())) + self._rest_fname_format = os.path.join( + out_dir, + '%s.rvic.rh%s.%%Y-%%m-%%d-%%H-%%M-%%S.nc' % + (self._caseid, self._tape_num)) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Determine when the next write should be and initialize out_data + self.__next_write_out_data() + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Determine when the update of out_data should be + self.__next_update_out_data() + # ------------------------------------------------------------ # + + log.debug(self.__repr__()) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Write a summary + def __str__(self): + return 'History Tape - {0}'.format(self.filename) + + def __repr__(self): + parts = ['------- Summary of History Tape Settings -------', + '\t# caseid: {0}'.format(self._caseid), + '\t# fincl: {0}'.format(','.join(self._fincl)), + '\t# nhtfrq: {0}'.format(self._nhtfrq), + '\t# mfilt: {0}'.format(self._mfilt), + '\t# ncprec: {0}'.format(self._ncprec), + '\t# avgflag: {0}'.format(self._avgflag), + '\t# fname_format: {0}'.format(self._fname_format), + '\t# file_format: {0}'.format(self._file_format), + '\t# outtype: {0}'.format(self._outtype), + '\t# out_dir: {0}'.format(self._out_dir), + '\t# calendar: {0}'.format(self._calendar), + '\t# units: {0}'.format(self._units), + ' ------- End of History Tape Settings -------'] + return '\n'.join(parts) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Update the history tapes with new fluxes + def update(self, data2tape, time_ord): + ''' Update the tape with new data''' + + # ------------------------------------------------------------ # + # Check that the time_ord is in sync + if self._time_ord != time_ord: + raise ValueError('rout_var.time_ord does not match the time_ord ' + 'passed in by the convolution call') + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Get the next timestamp + self._time_ord += self._dt / SECSPERDAY + self._timestamp = ord_to_datetime(self._time_ord, TIMEUNITS, + calendar=self._calendar) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Advance the Count + self._count += 1 + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Update the fields + for field in self._fincl: + tracer = 'LIQ' + log.debug('updating %s', field) + fdata = data2tape[field][tracer] + if self._avgflag == 'A': + self._temp_data[field] += fdata + elif self._avgflag == 'I': + if self._count == self._update_count: + self._temp_data[field] = fdata[:] + elif self._avgflag == 'X': + self._temp_data[field] = np.maximum(self._temp_data[field], + fdata) + elif self._avgflag == 'M': + self._temp_data[field] = np.minimum(self._temp_data[field], + fdata) + else: + raise ValueError('Average flag ({0}) does not match any of' + ' (A,I,X,M)'.format(self._avgflag)) + + # ------------------------------------------------------------ # + # If count == _update_count, add to _out_data + # Average first, if necessary + if (self._avgflag == 'A' and self._count == self._update_count): + self.__average() + + if self._count == self._update_count: + # move the data to the out_data structure + self.__update_out_data() + # Determine next update + self.__next_update_out_data() + + # zero out temp_data + for field in self._fincl: + self._temp_data[field][:] = 0.0 + # ------------------------------------------------------------ # + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + def write_initial(self): + pass + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + def __next_write_out_data(self): + '''determine the maximum size of out_data''' + + log.debug('determining size of out_data') + + self._out_data_i = 0 # position counter for out_data array + + # ------------------------------------------------------------ # + # b0 is first timestep of next period + # b1 is end of last timestep of next period + + # time when out_data will start (now) + b0 = self._time_ord + self._begtime = b0 + + # determine time when out_data will be full + if self._mfilt == 'year': + if self._nhtfrq == 0: + mfilt = MONTHSPERYEAR + else: + t1 = datetime(self._timestamp.year + 1, 1, 1) + b1 = date2num(t1, TIMEUNITS, calendar=self._calendar) + + # calculate the mfilt value + mfilt = int(round((b1 - b0) / self._out_data_stepsize)) + elif self._mfilt == 'month': + if self._nhtfrq == 0: + mfilt = 1 + else: + if self._timestamp.month == 12: + t1 = datetime(self._timestamp.year + 1, 2, 1) + else: + t1 = datetime(self._timestamp.year, + self._timestamp.month + 1, 1) + b1 = date2num(t1, TIMEUNITS, calendar=self._calendar) + + # calculate the mfilt value + mfilt = int(round((b1 - b0) / self._out_data_stepsize)) + elif self._mfilt == 'day': + if self._nhtfrq != 0: + b1 = b0 + 1.0 + else: + raise ValueError('Incompatable values for NHTFRQ and MFILT') + + # calculate the mfilt value + mfilt = int(round((b1 - b0) / self._out_data_stepsize)) + else: + mfilt = int(self._mfilt) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + if mfilt < 1: + mfilt = 1 + + self._out_data_write = mfilt - 1 + self._out_times = np.empty(mfilt, dtype=np.float64) + if self._avgflag != 'I': + self._out_time_bnds = np.empty((mfilt, 2), dtype=np.float64) + + shape = (mfilt, ) + self._out_data_shape + + log.debug('out_data shape: %s', shape) + log.debug('_out_data_write: %s', self._out_data_write) + + for field in self._fincl: + self._out_data[field] = np.zeros(shape, dtype=np.float64) + + self._out_data_has_values = False + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # fill in out_data + def __update_out_data(self): + + self._out_data_has_values = True + + # ------------------------------------------------------------ # + # Update the _out_data fields + for field in self._fincl: + if self._outtype == 'grid': + # ---------------------------------------------------- # + # Grid the fields + self._out_data[field][self._out_data_i, + self._outlet_y_ind, + self._outlet_x_ind] = \ + self._temp_data[field][:] * self._units_mult + # ---------------------------------------------------- # + else: + self._out_data[field][self._out_data_i, :] = \ + self._temp_data[field] * self._units_mult + + # Check that all values are valid, if not, exit gracefully + if np.isnan(self._out_data[field][self._out_data_i].sum()): + raise ValueError('nan found in output field: {0}, most likely ' + 'there is a nan/missing/fill value in the' + 'input forcings'.format(field)) + # ------------------------------------------------------------ # + + self._out_times[self._out_data_i] = self._write_ord + + if self._avgflag != 'I': + self._out_time_bnds[self._out_data_i, :] = self._time_bnds + + # ------------------------------------------------------------ # + # if out_data is full, write + if self._out_data_i == self._out_data_write: + self.finish() + self._out_data_i = 0 + + # Determine when the next write should be and initialize out_data + self.__next_write_out_data() + else: + self._out_data_i += 1 + log.debug('out_data counter is %s of %s', self._out_data_i, + self._out_data_write) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + def finish(self): + '''write out_data''' + log.debug('finishing tape %s', self._tape_num) + if self._out_data_has_values: + if self._outtype == 'grid': + self.__write_grid() + else: + self.__write_array() + + self.files_count += 1 + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Get import rvar fields + def __get_rvar(self, rvar): + ''' Get the rvar Fields that are useful for writing output ''' + self._dt = rvar.unit_hydrograph_dt + self._num_outlets = rvar.n_outlets + self._outlet_decomp_ind = rvar.outlet_decomp_ind + self._outlet_x_ind = rvar.outlet_x_ind + self._outlet_y_ind = rvar.outlet_y_ind + self._outlet_lon = rvar.outlet_lon + self._outlet_lat = rvar.outlet_lat + self._outlet_name = rvar.outlet_name + self._outlet_upstream_area = rvar.outlet_upstream_area + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Determine next write time + def __next_update_out_data(self): + ''' Determine the count for when the next write should occur ''' + # ------------------------------------------------------------ # + # If monthly, write at (YYYY,MM,1,0,0) + # b0 is first timestep of next period + # b1 is end of last timestep of next period + + b0 = self._time_ord + self._begtime = b0 + + if self._nhtfrq == 0: + if self._timestamp.month == 12: + b1 = date2num(datetime(self._timestamp.year + 1, 2, 1), + TIMEUNITS, calendar=self._calendar) + else: + b1 = date2num(datetime(self._timestamp.year, + self._timestamp.month + 1, 1), + TIMEUNITS, calendar=self._calendar) + + else: + b1 = b0 + self._out_data_stepsize + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Get the number of timesteps and datestamp for the next write + # next_ord is the ord_time when the write will happen + self._update_count = int(round((b1 - b0) / (self._dt / SECSPERDAY))) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Get next file names and timeord + if self._avgflag == 'I': + self._write_ord = b1 + self.filename = num2date( + b1, TIMEUNITS, + calendar=self._calendar).strftime(self._fname_format) + else: + self._time_bnds = np.array([[b0, b1]]) + self._write_ord = np.average(self._time_bnds) + self.filename = num2date( + b0, TIMEUNITS, + calendar=self._calendar).strftime(self._fname_format) + self.rest_filename = num2date( + b1, TIMEUNITS, + calendar=self._calendar).strftime(self._fname_format) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Set the count to zero + self._count = 0 + # ------------------------------------------------------------ # + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Average fields + def __average(self): + ''' Take the average based on the number of accumulated timesteps ''' + for field in self._fincl: + self._temp_data[field] /= self._count + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Write grid style history file + def __write_grid(self): + ''' Write history file ''' + + # ------------------------------------------------------------ # + # Open file + f = Dataset(self.filename, 'w', self._file_format) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Time Variable + f.createDimension('time', None) + + time = f.createVariable('time', self._ncprec, ('time',)) + time[:] = self._out_times[:self._out_data_i] + for key, val in iteritems(share.time): + if val: + setattr(time, key, val.encode()) + time.calendar = self._calendar.encode() + + if self._avgflag != 'I': + f.createDimension('nv', 2) + + time.bounds = 'time_bnds' + + time_bnds = f.createVariable('time_bnds', self._ncprec, + ('time', 'nv',), **self.ncvaropts) + time_bnds[:, :] = self._out_time_bnds[:self._out_data_i] + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Setup Coordinate Variables + if self._grid_lons.ndim > 1: + coords = ('yc', 'xc',) + + # Grid is not regular + xc = f.createDimension('xc', self._grid_lons.shape[1]) + yc = f.createDimension('yc', self._grid_lons.shape[0]) + + xc = f.createVariable('xc', self._ncprec, coords, **self.ncvaropts) + yc = f.createVariable('yc', self._ncprec, coords, **self.ncvaropts) + xc[:, :] = self._grid_lons + yc[:, :] = self._grid_lats + + for key, val in iteritems(share.xc): + if val: + setattr(xc, key, val.encode()) + + for key, val in iteritems(share.yc): + if val: + setattr(yc, key, val.encode()) + + else: + coords = ('lat', 'lon',) + + lon = f.createDimension('lon', len(self._grid_lons)) + lat = f.createDimension('lat', len(self._grid_lats)) + + lon = f.createVariable('lon', self._ncprec, ('lon',), + **self.ncvaropts) + lat = f.createVariable('lat', self._ncprec, ('lat',), + **self.ncvaropts) + lon[:] = self._grid_lons + lat[:] = self._grid_lats + + for key, val in iteritems(share.lon): + if val: + setattr(lon, key, val.encode()) + + for key, val in iteritems(share.lat): + if val: + setattr(lat, key, val.encode()) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Write Fields + tcoords = ('time',) + coords + + for field in self._fincl: + var = f.createVariable(field, self._ncprec, tcoords, + **self.ncvaropts) + var[:, :] = self._out_data[field][:self._out_data_i] + + for key, val in iteritems(getattr(share, field)): + if val: + setattr(var, key, val.encode()) + var.units = self._units.encode() + if self._grid_lons.ndim > 1: + var.coordinates = ' '.join(coords).encode() + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # write global attributes + self._glob_ats.update() + for key, val in iteritems(self._glob_ats.atts): + if val: + setattr(f, key, val.encode()) + # ------------------------------------------------------------ # + f.close() + log.info('Finished writing %s', self.filename) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Write array style history file + def __write_array(self): + ''' Write history file ''' + + # ------------------------------------------------------------ # + # Open file + f = Dataset(self.filename, 'w', self._file_format) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Time Variable + f.createDimension('time', None) + + time = f.createVariable('time', self._ncprec, ('time',), + **self.ncvaropts) + time[:] = self._out_times[:self._out_data_i] + for key, val in iteritems(share.time): + if val: + setattr(time, key, val.encode()) + time.calendar = self._calendar.encode() + + if self._avgflag != 'I': + f.createDimension('nv', 2) + + time.bounds = 'time_bnds' + + time_bnds = f.createVariable('time_bnds', self._ncprec, + ('time', 'nv',), **self.ncvaropts) + time_bnds[:, :] = self._out_time_bnds[:self._out_data_i] + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Setup Coordinate Variables + coords = ('outlets',) + + f.createDimension('outlets', self._num_outlets) + + nocoords = coords + ('nc_chars',) + char_names = stringtochar(self._outlet_name) + f.createDimension(nocoords[1], char_names.shape[1]) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Variables + outlet_lon = f.createVariable('lon', self._ncprec, coords, + **self.ncvaropts) + outlet_lat = f.createVariable('lat', self._ncprec, coords, + **self.ncvaropts) + outlet_x_ind = f.createVariable('outlet_x_ind', NC_INT, coords, + **self.ncvaropts) + outlet_y_ind = f.createVariable('outlet_y_ind', NC_INT, coords, + **self.ncvaropts) + outlet_decomp_ind = f.createVariable('outlet_decomp_ind', NC_INT, + coords, **self.ncvaropts) + onm = f.createVariable('outlet_name', NC_CHAR, nocoords, + **self.ncvaropts) + + outlet_lon[:] = self._outlet_lon + outlet_lat[:] = self._outlet_lat + outlet_x_ind[:] = self._outlet_x_ind + outlet_y_ind[:] = self._outlet_y_ind + outlet_decomp_ind[:] = self._outlet_decomp_ind + onm[:, :] = char_names + + for key, val in iteritems(share.outlet_lon): + if val: + setattr(outlet_lon, key, val.encode()) + + for key, val in iteritems(share.outlet_lat): + if val: + setattr(outlet_lat, key, val.encode()) + + for key, val in iteritems(share.outlet_y_ind): + if val: + setattr(outlet_y_ind, key, val.encode()) + + for key, val in iteritems(share.outlet_x_ind): + if val: + setattr(outlet_x_ind, key, val.encode()) + + for key, val in iteritems(share.outlet_decomp_ind): + if val: + setattr(outlet_decomp_ind, key, val.encode()) + + for key, val in iteritems(share.outlet_name): + if val: + setattr(onm, key, val.encode()) + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # Write Fields + tcoords = ('time',) + coords + + for field in self._fincl: + var = f.createVariable(field, self._ncprec, tcoords, + **self.ncvaropts) + var[:, :] = self._out_data[field][:self._out_data_i] + + for key, val in iteritems(getattr(share, field)): + if val: + setattr(var, key, val.encode()) + var.units = self._units.encode() + # ------------------------------------------------------------ # + + # ------------------------------------------------------------ # + # write global attributes + self._glob_ats.update() + for key, val in iteritems(self._glob_ats.atts): + if val: + setattr(f, key, val.encode()) + f.featureType = 'timeSeries'.encode() + # ------------------------------------------------------------ # + f.close() + log.info('Finished writing %s', self.filename) + # ---------------------------------------------------------------- # +# -------------------------------------------------------------------- # diff --git a/rvic/core/.ipynb_checkpoints/param_file-checkpoint.py b/rvic/core/.ipynb_checkpoints/param_file-checkpoint.py new file mode 100644 index 0000000..56f28c5 --- /dev/null +++ b/rvic/core/.ipynb_checkpoints/param_file-checkpoint.py @@ -0,0 +1,469 @@ +# -*- coding: utf-8 -*- +''' +param_file.py +''' +import numpy as np +import logging +from .log import LOG_NAME +from .write import write_param_file +from .share import NcGlobals, SECSPERDAY, MAX_NC_CHARS +from .pycompat import iteritems, pyrange, pyzip +from . import plots +import os +from datetime import date + + +# -------------------------------------------------------------------- # +# create logger +log = logging.getLogger(LOG_NAME) +# -------------------------------------------------------------------- # + + +# -------------------------------------------------------------------- # +# Wrap up functiions to finish the parameter file +def finish_params(outlets, dom_data, config_dict, directories): + ''' + Adjust the unit hydrographs and pack for parameter file + ''' + options = config_dict['OPTIONS'] + routing = config_dict['ROUTING'] + domain = config_dict['DOMAIN'] + dom_area = domain['AREA_VAR'] + dom_frac = domain['FRACTION_VAR'] + + if not len(outlets) > 0: + raise ValueError('outlets in finish_params are empty') + + # ------------------------------------------------------------ # + # netCDF variable options + ncvaropts = {} + if 'NETCDF_ZLIB' in options: + ncvaropts['zlib'] = options['NETCDF_ZLIB'] + if 'NETCDF_COMPLEVEL' in options: + ncvaropts['complevel'] = options['NETCDF_COMPLEVEL'] + if 'NETCDF_SIGFIGS' in options: + ncvaropts['least_significant_digit'] = options['NETCDF_SIGFIGS'] + # ------------------------------------------------------------ # + + # ---------------------------------------------------------------- # + # subset (shorten time base) + if options['SUBSET_DAYS'] and \ + options['SUBSET_DAYS'] < routing['BASIN_FLOWDAYS']: + subset_length = (options['SUBSET_DAYS'] * + SECSPERDAY / routing['OUTPUT_INTERVAL']) + outlets, full_time_length, \ + before, after = subset(outlets, subset_length=subset_length) + + slc = slice(min(len(before), 1000)) + + log.debug('plotting unit hydrograph timeseries now for before' + ' / after subseting') + + title = 'UHS before subset' + pfname = plots.uhs(before[slc], title, options['CASEID'], + directories['plots']) + log.info('%s Plot: %s', title, pfname) + + title = 'UHS after subset' + pfname = plots.uhs(after[slc], title, options['CASEID'], + directories['plots']) + log.info('%s Plot: %s', title, pfname) + else: + log.info('Not subsetting because either SUBSET_DAYS is null or ' + 'SUBSET_DAYS= 0, source_decomp_ind + assert outlet_decomp_ind.min() >= 0, outlet_decomp_ind + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Adjust Unit Hydrographs for differences in source/outlet areas and + # fractions + area = dom_data[domain['AREA_VAR']] + + if outlet_y_ind.ndim == 0 or outlet_x_ind.ndim == 0: + for source, outlet in enumerate(source2outlet_ind): + unit_hydrograph[:, source] *= area[source_y_ind[source], + source_x_ind[source]] + unit_hydrograph[:, source] /= area[outlet_y_ind[()], + outlet_x_ind[()]] + unit_hydrograph[:, source] *= frac_sources[source] + else: + for source, outlet in enumerate(source2outlet_ind): + unit_hydrograph[:, source] *= area[source_y_ind[source], + source_x_ind[source]] + unit_hydrograph[:, source] /= area[outlet_y_ind[outlet], + outlet_x_ind[outlet]] + unit_hydrograph[:, source] *= frac_sources[source] + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Make diagnostic plots + sum_after = np.zeros(dom_data[domain['FRACTION_VAR']].shape) + for i, (y, x) in enumerate(pyzip(source_y_ind, source_x_ind)): + sum_after[y, x] += unit_hydrograph[:, i].sum() + + plot_dict['Sum UH Final'] = sum_after + + dom_y = dom_data[domain['LATITUDE_VAR']] + dom_x = dom_data[domain['LONGITUDE_VAR']] + + for title, data in iteritems(plot_dict): + pfname = plots.fractions(data, dom_x, dom_y, title, options['CASEID'], + directories['plots']) + log.info('%s Plot: %s', title, pfname) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # fill in some misc arrays + if outlet_y_ind.ndim == 0: + numoutlets = 1 + else: + numoutlets = len(outlet_lon) + outlet_mask = np.zeros(numoutlets) + newshape = unit_hydrograph.shape + (1, ) + unit_hydrograph = unit_hydrograph.reshape(newshape) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Write parameter file + today = date.today().strftime('%Y%m%d') + param_file = os.path.join(directories['params'], + '{0}.rvic.prm.{1}.{2}.' + 'nc'.format(options['CASEID'], + options['GRIDID'], + today)) + + if 'NEW_DOMAIN' in list(config_dict.keys()): + dom_file_name = config_dict['NEW_DOMAIN']['FILE_NAME'] + else: + dom_file_name = config_dict['DOMAIN']['FILE_NAME'] + + param_file_name = \ + os.path.split(config_dict['POUR_POINTS']['FILE_NAME'])[1] + glob_atts = NcGlobals( + title='RVIC parameter file', + RvicPourPointsFile=param_file_name, + RvicUHFile=os.path.split(config_dict['UH_BOX']['FILE_NAME'])[1], + RvicFdrFile=os.path.split(routing['FILE_NAME'])[1], + RvicDomainFile=os.path.split(dom_file_name)[1]) + + log.debug('UH Range: (%f %f)', unit_hydrograph.min(), unit_hydrograph.max()) + + write_param_file(param_file, + nc_format=options['NETCDF_FORMAT'], + glob_atts=glob_atts, + full_time_length=full_time_length, + subset_length=subset_length, + unit_hydrograph_dt=routing['OUTPUT_INTERVAL'], + outlet_lon=outlet_lon, + outlet_lat=outlet_lat, + outlet_x_ind=outlet_x_ind, + outlet_y_ind=outlet_y_ind, + outlet_decomp_ind=outlet_decomp_ind, + outlet_number=outlet_number, + outlet_mask=outlet_mask, + outlet_name=outlet_name, + outlet_upstream_gridcells=outlet_upstream_gridcells, + outlet_upstream_area=outlet_upstream_area, + source_lon=source_lon, + source_lat=source_lat, + source_x_ind=source_x_ind, + source_y_ind=source_y_ind, + source_decomp_ind=source_decomp_ind, + source_time_offset=source_time_offset, + source2outlet_ind=source2outlet_ind, + unit_hydrograph=unit_hydrograph, + **ncvaropts) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # write a summary of what was done to the log file. + log.info('Parameter file includes %i outlets', len(outlets)) + log.info('Parameter file includes %i Source Points', len(source_lon)) + # ---------------------------------------------------------------- # + + return param_file, today +# -------------------------------------------------------------------- # + + +# -------------------------------------------------------------------- # +def adjust_fractions(outlets, dom_fractions, adjust=True): + ''' + Constrain the fractions in the outles. + The basic idea is that the sum of fraction from the outlets should not + exceed the domain fractions. + ''' + + log.info('Adjusting fractions now') + + fractions = np.zeros(dom_fractions.shape, dtype=np.float64) + ratio_fraction = np.ones(fractions.shape, dtype=np.float64) + adjusted_fractions = np.zeros(dom_fractions.shape, dtype=np.float64) + sum_uh_fractions = np.zeros(dom_fractions.shape, dtype=np.float64) + + # ---------------------------------------------------------------- # + # Aggregate the fractions + for key, outlet in iteritems(outlets): + y = outlet.y_source + x = outlet.x_source + + fractions[y, x] += outlet.fractions + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # First set fractions to zero where there is no land in the domain + yd, xd = np.nonzero(dom_fractions == 0.0) + fractions[yd, xd] = 0.0 + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Only adjust fractions where the aggregated fractions are gt the domain + # fractions + yi, xi = np.nonzero(fractions > dom_fractions) + log.info('Adjust fractions for %s grid cells', len(yi)) + ratio_fraction[yi, xi] = dom_fractions[yi, xi] / fractions[yi, xi] + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Adjust fracs based on ratio_fraction + for key, outlet in iteritems(outlets): + y = outlet.y_source + x = outlet.x_source + if adjust: + outlet.fractions *= ratio_fraction[y, x] + # For Diagnostics only + adjusted_fractions[y, x] += outlet.fractions + sum_uh_fractions[y, x] += outlet.unit_hydrograph.sum(axis=0) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # Make Fractions Dict for plotting + plot_dict = {'Domain Fractions': dom_fractions, + 'Aggregated Fractions': fractions, + 'Ratio Fractions': ratio_fraction, + 'Adjusted Fractions': adjusted_fractions, + 'Sum UH Before': sum_uh_fractions} + # ---------------------------------------------------------------- # + return outlets, plot_dict +# -------------------------------------------------------------------- # + + +# -------------------------------------------------------------------- # +# Shorten the unit hydrograph +def subset(outlets, subset_length=None): + ''' Shorten the Unit Hydrograph''' + + log.info('subsetting unit-hydrographs now...') + log.debug('Subset Length: %s', subset_length) + log.debug(outlets) + for i, (key, outlet) in enumerate(iteritems(outlets)): + if i == 0: + full_time_length = outlet.unit_hydrograph.shape[0] + log.debug('Subset Length: %s', subset_length) + log.debug('full_time_length: %s', full_time_length) + if not subset_length: + subset_length = full_time_length + log.debug('No subset_length provided, using full_time_length') + before = outlet.unit_hydrograph + else: + before = np.append(before, outlet.unit_hydrograph, + axis=1) + + outlet.offset = np.empty(outlet.unit_hydrograph.shape[1], + dtype=np.int32) + out_uh = np.zeros((subset_length, outlet.unit_hydrograph.shape[1]), + dtype=np.float64) + + d_left = -1 * subset_length / 2 + d_right = subset_length / 2 + + for j in pyrange(outlet.unit_hydrograph.shape[1]): + # find index position of maximum + maxind = np.argmax(outlet.unit_hydrograph[:, j]) + + # find bounds + left = maxind + d_left + right = maxind + d_right + + # make sure left and right fit in unit hydrograph array, + # if not adjust + if left < 0: + left = 0 + right = subset_length + if right > full_time_length: + right = full_time_length + left = full_time_length - subset_length + + log.warning('Subset centered on UH max extends beyond length ' + 'of unit hydrograph.') + log.warning('--> Outlet %s', outlet) + log.warning('----> Max Index is %s', maxind) + log.warning('----> Last value in subset ' + 'is %s', outlet.unit_hydrograph[-1, j]) + if maxind == full_time_length: + log.warning('maxind == full_time_length, not able to ' + 'resolve unithydrograph') + if left < 0 or right > full_time_length: + raise ValueError('Subsetting failed left:{0} or right {1} does' + ' not fit inside bounds'.format(left, right)) + + outlet.offset[j] = left + + # clip and normalize + tot = outlet.unit_hydrograph[left:right, j].sum() + out_uh[:, j] = outlet.unit_hydrograph[left:right, j] / tot + + outlet.unit_hydrograph = out_uh + + if i == 0: + after = outlet.unit_hydrograph + else: + after = np.append(after, outlet.unit_hydrograph, axis=1) + + log.info('Done subsetting') + + return outlets, full_time_length, before, after +# -------------------------------------------------------------------- # + + +# -------------------------------------------------------------------- # +def group(outlets, subset_length): + ''' + group the outlets into one set of arrays + ''' + + n_outlets = len(outlets) + n_sources = 0 + for key, outlet in iteritems(outlets): + n_sources += len(outlet.y_source) + + gd = {} + + log.debug('n_outlets: %s', n_outlets) + log.debug('n_sources: %s', n_sources) + log.debug('subset_length: %s', subset_length) + + # ---------------------------------------------------------------- # + # Source specific values + gd['unit_hydrograph'] = np.empty((subset_length, n_sources), + dtype=np.float64) + gd['frac_sources'] = np.empty(n_sources, dtype=np.float64) + gd['source_lon'] = np.empty(n_sources, dtype=np.float64) + gd['source_lat'] = np.empty(n_sources, dtype=np.float64) + gd['source_x_ind'] = np.empty(n_sources, dtype=np.int32) + gd['source_y_ind'] = np.empty(n_sources, dtype=np.int32) + gd['source_decomp_ind'] = np.empty(n_sources, dtype=np.int32) + gd['source_time_offset'] = np.empty(n_sources, dtype=np.int32) + gd['source2outlet_ind'] = np.empty(n_sources, dtype=np.int32) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # outlet specific inputs + gd['outlet_lon'] = np.empty(n_outlets, dtype=np.float64) + gd['outlet_lat'] = np.empty(n_outlets, dtype=np.float64) + gd['outlet_x_ind'] = np.empty(n_outlets, dtype=np.int32) + gd['outlet_y_ind'] = np.empty(n_outlets, dtype=np.int32) + gd['outlet_decomp_ind'] = np.empty(n_outlets, dtype=np.int32) + gd['outlet_number'] = np.empty(n_outlets, dtype=np.int32) + gd['outlet_name'] = np.empty(n_outlets, dtype='S{0}'.format(MAX_NC_CHARS)) + gd['outlet_upstream_gridcells'] = np.empty(n_outlets, dtype=np.int32) + gd['outlet_upstream_area'] = np.empty(n_outlets, dtype=np.float64) + # ---------------------------------------------------------------- # + + # ---------------------------------------------------------------- # + # place outlet and source vars into gd dictionary + a = 0 + for i, (key, outlet) in enumerate(iteritems(outlets)): + b = a + len(outlet.y_source) + log.debug('%s unit_hydrograph.shape %s', outlet.name, + outlet.unit_hydrograph.shape) + # -------------------------------------------------------- # + # Point specific values + gd['unit_hydrograph'][:, a:b] = outlet.unit_hydrograph + gd['frac_sources'][a:b] = outlet.fractions + gd['source_lon'][a:b] = outlet.lon_source + gd['source_lat'][a:b] = outlet.lat_source + gd['source_x_ind'][a:b] = outlet.x_source + gd['source_y_ind'][a:b] = outlet.y_source + gd['source_decomp_ind'][a:b] = outlet.cell_id_source + gd['source_time_offset'][a:b] = outlet.offset + gd['source2outlet_ind'][a:b] = i + # -------------------------------------------------------- # + + # -------------------------------------------------------- # + # outlet specific inputs + gd['outlet_lon'][i] = outlet.lon + gd['outlet_lat'][i] = outlet.lat + gd['outlet_x_ind'][i] = outlet.domx + gd['outlet_y_ind'][i] = outlet.domy + gd['outlet_decomp_ind'][i] = outlet.cell_id + gd['outlet_number'][i] = i + gd['outlet_name'][i] = outlet.name + gd['outlet_upstream_gridcells'][i] = outlet.upstream_gridcells + gd['outlet_upstream_area'][i] = outlet.upstream_area + # -------------------------------------------------------- # + + # -------------------------------------------------------- # + # update src counter + a = b + # -------------------------------------------------------- # + # ---------------------------------------------------------------- # + + return gd +# -------------------------------------------------------------------- # diff --git a/rvic/core/.ipynb_checkpoints/share-checkpoint.py b/rvic/core/.ipynb_checkpoints/share-checkpoint.py new file mode 100644 index 0000000..a5378a4 --- /dev/null +++ b/rvic/core/.ipynb_checkpoints/share-checkpoint.py @@ -0,0 +1,345 @@ +# -*- coding: utf-8 -*- +''' +share.py +''' +import sys +import socket +import string +import time as time_mod +from .pycompat import OrderedDict, iteritems +from netCDF4 import default_fillvals +from getpass import getuser + + +# ----------------------- CONSTANTS --------------------------------- # +EARTHRADIUS = 6.37122e6 # meters +WATERDENSITY = 1000. + +# area +METERSPERKM = 1000. +METERSPERMILE = 1609.34 +METERS2PERACRE = 4046.856 + +# time +# reference time +REFERENCE_STRING = '0001-1-1 0:0:0' +REFERENCE_DATE = 10101 # i.e. REFERENCE_STRING +REFERENCE_TIME = 0 # i.e. REFERENCE_STRING +TIMEUNITS = 'days since ' + REFERENCE_STRING # do not change (MUST BE DAYS)! +TIMESTAMPFORM = '%Y-%m-%d-%H' +CALENDAR = 'noleap' +HOURSPERDAY = 24. +SECSPERHOUR = 3600. +MINSPERHOUR = 60. +MINSPERDAY = HOURSPERDAY * MINSPERHOUR +SECSPERDAY = HOURSPERDAY * SECSPERHOUR +MONTHSPERYEAR = 12 + +# length +MMPERMETER = 1000. +CMPERMETER = 100. + +# precision +PRECISION = 1.0e-30 +NC_DOUBLE = 'f8' +NC_FLOAT = 'f4' +NC_INT = 'i4' +NC_CHAR = 'S1' +MAX_NC_CHARS = 256 + +# fill values +FILLVALUE_F = default_fillvals[NC_DOUBLE] +FILLVALUE_I = default_fillvals[NC_INT] + +# filenames +RPOINTER = 'rpointer' + +# tracers +RVIC_TRACERS = ('LIQ',) # Before changing, update history module + +# Calendar key number for linking with CESM +CALENDAR_KEYS = {0: ['None'], + 1: ['noleap', '365_day'], + 2: ['gregorian', 'standard'], + 3: ['proleptic_gregorian'], + 4: ['all_leap', '366_day'], + 5: ['360_day'], + 6: ['julian']} + +VALID_CHARS = '-_. %s%s' % (string.ascii_letters, string.digits) + + +# ----------------------- NETCDF VARIABLES --------------------------------- # +class NcGlobals(object): + def __init__(self, + title=None, + casename=None, + casestr=None, + history='Created: {}'.format(time_mod.ctime(time_mod.time())), + institution='University of Washington', + source=sys.argv[0], + references='Based on the initial model of Lohmann, et al., ' + '1996, Tellus, 48(A), 708-721', + comment='Output from the RVIC Streamflow Routing Model.', + Conventions='CF-1.6', + RvicPourPointsFile=None, + RvicFdrFile=None, + RvicUHFile=None, + RvicDomainFile=None, + version=None, + hostname=None, + username=None): + + self.atts = OrderedDict() + + if title is not None: + self.atts['title'] = title + + if comment is not None: + self.atts['comment'] = comment + + if Conventions is not None: + self.atts['Conventions'] = Conventions + + if history is not None: + self.atts['history'] = history + + if source is not None: + self.atts['source'] = source + + if institution is not None: + self.atts['institution'] = institution + + if hostname is not None: + self.atts['hostname'] = hostname + else: + self.atts['hostname'] = socket.gethostname() + + if username is not None: + self.atts['username'] = username + else: + self.atts['username'] = getuser() + + if casename is not None: + self.atts['casename'] = casename + + if casestr is not None: + self.atts['casestr'] = casestr + + if references is not None: + self.atts['references'] = references + + if version is not None: + self.atts['version'] = version + else: + try: + from rvic import version + self.atts['version'] = version.short_version + except ImportError: + self.atts['version'] = 'unknown' + + if RvicPourPointsFile is not None: + self.atts['RvicPourPointsFile'] = RvicPourPointsFile + + if RvicUHFile is not None: + self.atts['RvicUHFile'] = RvicUHFile + + if RvicFdrFile is not None: + self.atts['RvicFdrFile'] = RvicFdrFile + + if RvicDomainFile is not None: + self.atts['RvicDomainFile'] = RvicDomainFile + + def update(self): + self.atts['history'] = \ + 'Created: {0}'.format(time_mod.ctime(time_mod.time())) + + +# Coordinate Variables +time = dict(long_name='time', + units=TIMEUNITS) + +time_bnds = dict() + +timesteps = dict(long_name='Series of timesteps', + nits='unitless') + +lon = dict(long_name='longitude', + nits='degrees_east') + +lat = dict(long_name='latitude', + units='degrees_north') + +xc = dict(long_name='longitude', + units='degrees_east') + +yc = dict(long_name='latitude', + units='degrees_north') + +# Data Variables +fraction = dict(long_name='fraction of grid cell that is active', + units='unitless') + +unit_hydrograph = dict(long_name='Unit Hydrograph', + units='unitless') + +avg_velocity = dict(long_name='Flow Velocity Parameter', + units='m s-1') + +avg_diffusion = dict(long_name='Diffusion Parameter', + units='m2 s-1') + +global_basin_id = dict(long_name='Global Basin ID from RvicFdrFile', + units='unitless') + +full_time_length = dict(long_name='Length of original unit hydrograph', + units='timesteps') + +subset_length = dict(long_name='Shortened length of the unit hydrograph', + units='timesteps') + +unit_hydrograph_dt = dict(long_name='Unit hydrograph timestep', + units='seconds') + +outlet_x_ind = dict(long_name='x grid coordinate of outlet grid cell', + units='unitless') + +outlet_y_ind = dict(long_name='y grid coordinate of outlet grid cell', + units='unitless') + +outlet_lon = dict(long_name='Longitude coordinate of outlet grid cell', + units='degrees_east') + +outlet_lat = dict(long_name='Latitude coordinate of outlet grid cell', + units='degrees_north') + +outlet_decomp_ind = dict(long_name='1d grid location of outlet grid cell', + units='unitless') + +outlet_number = dict(long_name='outlet number', + units='unitless') + +outlet_mask = dict(long_name='type of outlet point', + units='0-ocean, 1-land, 2-guage, 3-none') + +outlet_name = dict(long_name='Outlet guage name', + units='unitless') + +outlet_upstream_area = dict(long_name='Upstream catchment area contributing ' + 'to outlet', + units='m2') + +outlet_upstream_gridcells = dict(long_name='Number of upstream grid cells ' + 'contributing to outlet', + units='number of grid cells') + +source_x_ind = dict(long_name='x grid coordinate of source grid cell', + units='unitless') + +source_y_ind = dict(long_name='y grid coordinate of source grid cell', + units='unitless') + +source_lon = dict(long_name='Longitude coordinate of source grid cell', + units='degrees_east') + +source_lat = dict(long_name='Latitude coordinate of source grid cell', + units='degrees_north') + +source_decomp_ind = dict(long_name='1d grid location of source grid cell', + units='unitless') +source_time_offset = dict(long_name='Number of leading timesteps ommited', + units='timesteps') + +source2outlet_ind = dict(long_name='source to outlet index mapping', + units='unitless') + +ring = dict(long_name='Convolution Ring', + units='kg m-2 s-1') + +streamflow = dict(long_name='Streamflow at outlet grid cell', + units='kg m-2 s-1') + +storage = dict(long_name='Mass storage in stream upstream of outlet ' + 'grid cell', + units='kg m-2 s-1') + +# valid values http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.6/\ +# cf-conventions.html#calendar +timemgr_rst_type = dict(long_name='calendar type', + units='unitless', + flag_values='0, 1, 2, 3, 4, 5, 6', + flag_meanings='NONE, NO_LEAP_C, GREGORIAN, ' + 'PROLEPTIC_GREGORIAN, ALL_LEAP, ' + '360_DAY, JULIAN') + +timemgr_rst_step_sec = dict(long_name='seconds component of timestep size', + units='sec', + valid_range='0, 86400') + +timemgr_rst_start_ymd = dict(long_name='start date', + units='YYYYMMDD') + +timemgr_rst_start_tod = dict(long_name='start time of day', + units='sec', + valid_range='0, 86400') + +timemgr_rst_ref_ymd = dict(long_name='reference date', + units='YYYYMMDD') + +timemgr_rst_ref_tod = dict(long_name='reference time of day', + units='sec', + valid_range='0, 86400') + +timemgr_rst_curr_ymd = dict(long_name='current date', + units='YYYYMMDD') + +timemgr_rst_curr_tod = dict(long_name='current time of day', + units='sec', + valid_range='0, 86400') + +nhtfrq = dict(long_name='Frequency of history writes', + units='absolute value of negative is in hours, 0=monthly, ' + 'positive is time-steps', + comment='Namelist item') + +mfilt = dict(long_name='Number of history time samples on a file', + units='initless', + comment='Namelist item') + +ncprec = dict(long_name='Flag for data precision', + flag_values='1, 2', + flag_meanings='single-precision double-precision', + comment='Namelist item', + valid_range='1, 2') + +fincl = dict(long_name='Fieldnames to include', + comment='Namelist item') + +fexcl = dict(long_name='Fieldnames to exclude', + comment='Namelist item') + +nflds = dict(long_name='Number of fields on file', + units='unitless') + +ntimes = dict(long_name='Number of time steps on file', + units='time-step') +is_endhist = dict(long_name='End of history file', + flag_values='0, 1', + flag_meanings='FALSE TRUE', + comment='Namelist item', + valid_range='0, 1') + +begtime = dict(long_name='Beginning time', + units='time units') + +hpindex = dict(long_name='History pointer index', + units='units') + +avgflag = dict(long_name='Averaging flag', + units='A=Average, X=Maximum, M=Minimum, I=Instantaneous') + +name = dict(long_name='Fieldnames') + +long_name = dict(long_name='Long descriptive names for fields') + +units = dict(long_name='Units for each history field output') diff --git a/rvic/core/history.py b/rvic/core/history.py index df11bec..897d965 100644 --- a/rvic/core/history.py +++ b/rvic/core/history.py @@ -645,7 +645,7 @@ def __write_array(self): f.createDimension('outlets', self._num_outlets) nocoords = coords + ('nc_chars',) - char_names = stringtochar(self._outlet_name) + char_names = self._outlet_name.astype(str) f.createDimension(nocoords[1], char_names.shape[1]) # ------------------------------------------------------------ # diff --git a/rvic/core/param_file.py b/rvic/core/param_file.py index 4f9d652..6049266 100644 --- a/rvic/core/param_file.py +++ b/rvic/core/param_file.py @@ -49,7 +49,7 @@ def finish_params(outlets, dom_data, config_dict, directories): # subset (shorten time base) if options['SUBSET_DAYS'] and \ options['SUBSET_DAYS'] < routing['BASIN_FLOWDAYS']: - subset_length = int(options['SUBSET_DAYS'] * + subset_length = (options['SUBSET_DAYS'] * SECSPERDAY / routing['OUTPUT_INTERVAL']) outlets, full_time_length, \ before, after = subset(outlets, subset_length=subset_length) @@ -327,11 +327,10 @@ def subset(outlets, subset_length=None): outlet.offset = np.empty(outlet.unit_hydrograph.shape[1], dtype=np.int32) - out_uh = np.zeros((subset_length, outlet.unit_hydrograph.shape[1]), - dtype=np.float64) + out_uh = np.zeros((int(subset_length), int(outlet.unit_hydrograph.shape[1]))) - d_left = int(-1 * subset_length / 2) - d_right = int(subset_length / 2) + d_left = -1 * subset_length / 2 + d_right = subset_length / 2 for j in pyrange(outlet.unit_hydrograph.shape[1]): # find index position of maximum @@ -366,8 +365,8 @@ def subset(outlets, subset_length=None): outlet.offset[j] = left # clip and normalize - tot = outlet.unit_hydrograph[left:right, j].sum() - out_uh[:, j] = outlet.unit_hydrograph[left:right, j] / tot + tot = outlet.unit_hydrograph[int(left):int(right),int(j)].sum() + out_uh[:, j] = outlet.unit_hydrograph[int(left):int(right), int(j)] / tot outlet.unit_hydrograph = out_uh @@ -401,16 +400,16 @@ def group(outlets, subset_length): # ---------------------------------------------------------------- # # Source specific values - gd['unit_hydrograph'] = np.empty((subset_length, n_sources), + gd['unit_hydrograph'] = np.empty((int(subset_length), int(n_sources)), dtype=np.float64) - gd['frac_sources'] = np.empty(n_sources, dtype=np.float64) - gd['source_lon'] = np.empty(n_sources, dtype=np.float64) - gd['source_lat'] = np.empty(n_sources, dtype=np.float64) - gd['source_x_ind'] = np.empty(n_sources, dtype=np.int32) - gd['source_y_ind'] = np.empty(n_sources, dtype=np.int32) - gd['source_decomp_ind'] = np.empty(n_sources, dtype=np.int32) - gd['source_time_offset'] = np.empty(n_sources, dtype=np.int32) - gd['source2outlet_ind'] = np.empty(n_sources, dtype=np.int32) + gd['frac_sources'] = np.empty(int(n_sources), dtype=np.float64) + gd['source_lon'] = np.empty(int(n_sources), dtype=np.float64) + gd['source_lat'] = np.empty(int(n_sources), dtype=np.float64) + gd['source_x_ind'] = np.empty(int(n_sources), dtype=np.int32) + gd['source_y_ind'] = np.empty(int(n_sources), dtype=np.int32) + gd['source_decomp_ind'] = np.empty(int(n_sources), dtype=np.int32) + gd['source_time_offset'] = np.empty(int(n_sources), dtype=np.int32) + gd['source2outlet_ind'] = np.empty(int(n_sources), dtype=np.int32) # ---------------------------------------------------------------- # # ---------------------------------------------------------------- # diff --git a/rvic/core/share.py b/rvic/core/share.py index 14a58f6..466e89d 100644 --- a/rvic/core/share.py +++ b/rvic/core/share.py @@ -274,28 +274,28 @@ def update(self): timemgr_rst_step_sec = dict(long_name='seconds component of timestep size', units='sec', - range=[0, 86400]) + range='0, 86400') timemgr_rst_start_ymd = dict(long_name='start date', units='YYYYMMDD') timemgr_rst_start_tod = dict(long_name='start time of day', units='sec', - range=[0, 86400]) + range='0, 86400') timemgr_rst_ref_ymd = dict(long_name='reference date', units='YYYYMMDD') timemgr_rst_ref_tod = dict(long_name='reference time of day', units='sec', - range=[0, 86400]) + range='0, 86400') timemgr_rst_curr_ymd = dict(long_name='current date', units='YYYYMMDD') timemgr_rst_curr_tod = dict(long_name='current time of day', units='sec', - range=[0, 86400]) + range='0, 86400') nhtfrq = dict(long_name='Frequency of history writes', units='absolute value of negative is in hours, 0=monthly, ' @@ -310,7 +310,7 @@ def update(self): flag_values='1, 2', flag_meanings='single-precision double-precision', comment='Namelist item', - valid_range=[1, 2]) + valid_range='1, 2') fincl = dict(long_name='Fieldnames to include', comment='Namelist item') @@ -327,7 +327,7 @@ def update(self): flag_values='0, 1', flag_meanings='FALSE TRUE', comment='Namelist item', - valid_range=[0, 1]) + valid_range='0, 1') begtime = dict(long_name='Beginning time', units='time units') diff --git a/rvic/parameters.py b/rvic/parameters.py index e98c7aa..3266d1e 100644 --- a/rvic/parameters.py +++ b/rvic/parameters.py @@ -191,7 +191,7 @@ def gen_uh_init(config): if 'names' in pour_points: pour_points.fillna(inplace=True, value='unknown') for i, name in enumerate(pour_points.names): - pour_points.ix[i, 'names'] = strip_invalid_char(name) + pour_points.loc[name, 'names'] = strip_invalid_char(name) pour_points.drop_duplicates(inplace=True) pour_points.dropna() diff --git a/rvic/version.py b/rvic/version.py index e284ae3..cffe302 100644 --- a/rvic/version.py +++ b/rvic/version.py @@ -1,2 +1,2 @@ -version = '1.1.0' -short_version = '1.1.0' +version = '1.1.1' +short_version = '1.1.1' diff --git a/samples/configs/rvic.convolution.rasm.cfg b/samples/configs/rvic.convolution.rasm.cfg index ac2d35e..bda1c35 100644 --- a/samples/configs/rvic.convolution.rasm.cfg +++ b/samples/configs/rvic.convolution.rasm.cfg @@ -133,7 +133,7 @@ FILE_NAME: None [PARAM_FILE] #--rvic parameter file file (char) --> -FILE_NAME: ./samples/cases/sample_rasm_parameters/params/sample_rasm_parameters.rvic.prm.wr50a.20151024.nc +FILE_NAME: ./samples/flow_directions/pnw.RVIC.input_20140218.nc #-- ====================================== --# [INPUT_FORCINGS] diff --git a/samples/pour_points/.ipynb_checkpoints/columbia_sample_pour_points-checkpoint.csv b/samples/pour_points/.ipynb_checkpoints/columbia_sample_pour_points-checkpoint.csv new file mode 100644 index 0000000..c76fda6 --- /dev/null +++ b/samples/pour_points/.ipynb_checkpoints/columbia_sample_pour_points-checkpoint.csv @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ce13e2314839c6c2a560db0381b32df6e09300f35a0f4591dc256e0690f9610e +size 518 diff --git a/samples/pour_points/.ipynb_checkpoints/rasm_sample_pour_points-checkpoint.csv b/samples/pour_points/.ipynb_checkpoints/rasm_sample_pour_points-checkpoint.csv new file mode 100644 index 0000000..2acaaab --- /dev/null +++ b/samples/pour_points/.ipynb_checkpoints/rasm_sample_pour_points-checkpoint.csv @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:16458c8ced7dd66fb6116a81adb42ce2742cfc4b97b5d6f6788b90ce6668e735 +size 289 diff --git a/samples/pour_points/columbia_sample_pour_points.csv b/samples/pour_points/columbia_sample_pour_points.csv index c76fda6..2ed6d86 100644 --- a/samples/pour_points/columbia_sample_pour_points.csv +++ b/samples/pour_points/columbia_sample_pour_points.csv @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:ce13e2314839c6c2a560db0381b32df6e09300f35a0f4591dc256e0690f9610e -size 518 +oid sha256:0ee4f18f345d9e78ad5d988c6d264ac7e893e5db1981cf499b3c951402271a78 +size 547 diff --git a/samples/pour_points/rasm_sample_pour_points.csv b/samples/pour_points/rasm_sample_pour_points.csv index 2acaaab..4c4dfa0 100644 --- a/samples/pour_points/rasm_sample_pour_points.csv +++ b/samples/pour_points/rasm_sample_pour_points.csv @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:16458c8ced7dd66fb6116a81adb42ce2742cfc4b97b5d6f6788b90ce6668e735 -size 289 +oid sha256:1ca4a2db86df7f40db26610ac0ac0cf14974a113a779a33c2054151ef75c8e6e +size 79 diff --git a/scripts/.ipynb_checkpoints/rvic-checkpoint b/scripts/.ipynb_checkpoints/rvic-checkpoint new file mode 100644 index 0000000..69ec3ea --- /dev/null +++ b/scripts/.ipynb_checkpoints/rvic-checkpoint @@ -0,0 +1,62 @@ +#!/usr/bin/env python +"""RVIC command line interface""" +import argparse +import os.path + + +# -------------------------------------------------------------------- # +def main(): + """ + Get the script and path to the config_file + """ + # Parse arguments + parser = argparse.ArgumentParser(description='The RVIC streamflow routing ' + 'model is based on the original model of ' + 'Lohmann, et al., 1996, Tellus, 48(A), ' + '708-721') + parser.add_argument("script", type=str, + help="RVIC subprogram to run", + choices=['parameters', 'convolution', 'convert'], + default=None, nargs='?') + parser.add_argument("config_file", type=str, + help="Input configuration file", + default=None, nargs='?') + parser.add_argument("-np", "--numofproc", type=int, + help="Number of processors used to run job", default=1) + parser.add_argument("--version", action='store_true', + help="Return RVIC version string") + + args = parser.parse_args() + + if args.version: + from rvic import version + print(version.short_version) + return + + if (args.script and args.config_file): + if not os.path.isfile(args.config_file): + raise IOError('config_file: {0} does not ' + 'exist'.format(args.config_file)) + + if args.script == 'parameters': + from rvic.parameters import parameters + parameters(args.config_file, numofproc=args.numofproc) + else: + if args.numofproc > 1: + print('{0} processors were specified but script {1} only ' + 'excepts 1'.format(args.numofproc, args.script)) + if args.script == 'convolution': + from rvic.convolution import convolution + convolution(args.config_file) + if args.script == 'convert': + from rvic.convert import convert + convert(args.config_file) + else: + parser.print_help() + return +# -------------------------------------------------------------------- # + +# -------------------------------------------------------------------- # +if __name__ == "__main__": + main() +# -------------------------------------------------------------------- #