Merge pull request #35 from jhamman/hotfix.1.0.2

Added the 'SEARCH_FOR_CHANNEL' option...

Author: Joe Hamman

rvic/core/aggregate.py

@@ -103,7 +103,7 @@ def make_agg_pairs(pour_points, dom_data, fdr_data, config_dict):
     log.info('NUMBER OF POINTS TO AGGREGATE TO: %i' % key_count)
     log.info('NUMBER OF POUR POINTS AGGREGATED: %i' % pp_count)
     log.info('EFFECIENCY OF: %.2f %%' % (100.*pp_count / num))
-    log.info('UNASSIGNED POUR POINTS: %i \n' % (num-pp_count))
+    log.info('UNASSIGNED POUR POINTS: %i' % (num-pp_count))
     log.info('-------------------------------------------------------------\n')
 
     # ---------------------------------------------------------------- #

rvic/core/param_file.py

@@ -124,14 +124,15 @@ def finish_params(outlets, dom_data, config_dict, directories):
     # fractions
     area = dom_data[domain['AREA_VAR']]
 
-    for source, outlet in enumerate(source2outlet_ind):
-        if outlet_y_ind.ndim == 0 or outlet_x_ind.ndim == 0:
+    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:
+    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],
@@ -412,7 +413,8 @@ def group(outlets, subset_length):
     a = 0
     for i, (cell_id, outlet) in enumerate(outlets.iteritems()):
         b = a + len(outlet.y_source)
-        log.debug('unit_hydrograph.shape %s', outlet.unit_hydrograph.shape)
+        log.debug('outlet.unit_hydrograph.shape %s',
+                  outlet.unit_hydrograph.shape)
         # -------------------------------------------------------- #
         # Point specific values
         gd['unit_hydrograph'][:, a:b] = outlet.unit_hydrograph

rvic/core/plots.py

@@ -13,7 +13,7 @@
     matplotlib_available = True
     try:
         from mpl_toolkits.basemap import Basemap
-        basemap_available = False  # Core dumping due to Geos Error
+        basemap_available = False
     except:
         basemap_available = False
 except:
@@ -32,9 +32,9 @@ def uhs(data, title, case_id, plot_dir):
     Plot diagnostic plot showing all unit hydrographs
     """
     today = date.today().strftime('%Y%m%d')
-    file_name = "{}_{}_{}.png".format(title.lower().replace(" ", "_"),
-                                      case_id.lower().replace(" ", "_"),
-                                      today)
+    file_name = "{0}_{1}_{2}.png".format(title.lower().replace(" ", "_"),
+                                         case_id.lower().replace(" ", "_"),
+                                         today)
     pfname = os.path.join(plot_dir, file_name)
 
     fig = plt.figure()
@@ -56,9 +56,9 @@ def _fractions_grid(data, dom_x, dom_y, title, case_id, plot_dir):
     # ---------------------------------------------------------------- #
     # Plot Fractions
     today = date.today().strftime('%Y%m%d')
-    file_name = "{}_{}_{}.png".format(title.lower().replace(" ", "_"),
-                                      case_id.lower().replace(" ", "_"),
-                                      today)
+    file_name = "{0}_{1}_{2}.png".format(title.lower().replace(" ", "_"),
+                                         case_id.lower().replace(" ", "_"),
+                                         today)
     pfname = os.path.join(plot_dir, file_name)
 
     mask = data <= 0.0
@@ -75,7 +75,7 @@ def _fractions_grid(data, dom_x, dom_y, title, case_id, plot_dir):
     plt.title(title)
     plt.xlabel('x')
     plt.ylabel('y')
-    plt.gca().invert_yaxis()
+    # plt.gca().invert_yaxis()
     fig.savefig(pfname)
     # ---------------------------------------------------------------- #
     return pfname
@@ -90,9 +90,9 @@ def _fractions_map(data, dom_x, dom_y, title, case_id, plot_dir):
     # ---------------------------------------------------------------- #
     # Plot Fractions
     today = date.today().strftime('%Y%m%d')
-    file_name = "{}_{}_{}.png".format(title.lower().replace(" ", "_"),
-                                      case_id.lower().replace(" ", "_"),
-                                      today)
+    file_name = "{0}_{1}_{2}.png".format(title.lower().replace(" ", "_"),
+                                         case_id.lower().replace(" ", "_"),
+                                         today)
     pfname = os.path.join(plot_dir, file_name)
 
     fig = plt.figure(figsize=(8, 8))
@@ -109,7 +109,6 @@ def _fractions_map(data, dom_x, dom_y, title, case_id, plot_dir):
     # define projection
     midx = int(dom_x.shape[1]/2)
     midy = int(dom_x.shape[0]/2)
-    print midx, midy
 
     projection = {'projection': 'stere',
                   'lon_0': dom_x[-1, midx],
@@ -136,7 +135,6 @@ def _fractions_map(data, dom_x, dom_y, title, case_id, plot_dir):
     m.drawmeridians(meridians, labels=[0, 0, 0, 1])
 
     x, y = m(dom_x, dom_y)  # compute map proj coordinates.
-    print 'passed 1'
     cs = m.pcolormesh(x, y, data, cmap=cmap)
     m.colorbar(cs, location='right', pad="5%")
     plt.title(title)

rvic/core/utilities.py

@@ -50,6 +50,40 @@ def latlon2yx(plats, plons, glats, glons):
 # -------------------------------------------------------------------- #
 
 
+# -------------------------------------------------------------------- #
+# Search neighboring grid cells for channel
+def search_for_channel(source_area, routys, routxs, search=2, tol=10):
+    """Search neighboring grid cells for channel"""
+
+    log.debug('serching for channel')
+
+    new_ys = np.empty_like(routys)
+    new_xs = np.empty_like(routxs)
+
+    for i, (y, x) in enumerate(zip(routys, routxs)):
+        area0 = source_area[y, x]
+
+        search_area = source_area[y-search:y+search+1, x-search:x+search+1]
+
+        if np.any(search_area > area0*tol):
+            sy, sx = np.unravel_index(search_area.argmax(), search_area.shape)
+
+            new_ys[i] = y + sy - search
+            new_xs[i] = x + sx - search
+
+            log.debug('Moving pour point to channel y: '
+                      '{0}->{1}, x: {2}->{3}'.format(y, new_ys[i],
+                                                     x, new_xs[i]))
+            log.debug('Source Area has increased from {0}'
+                      ' to {1}'.format(area0, source_area[new_ys[i], new_xs[i]]))
+        else:
+            new_ys[i] = y
+            new_xs[i] = x
+    return new_ys, new_xs
+
+# -------------------------------------------------------------------- #
+
+
 # -------------------------------------------------------------------- #
 # Write rpointer file
 def write_rpointer(restart_dir, restart_file, timestamp):

rvic/core/variables.py

@@ -27,7 +27,7 @@ class Point(object):
     '''
 
     def __init__(self, lat=-9999.9, lon=-9999.9, domx=-9999, domy=-9999,
-                 routx=-9999, routy=-9999, name='', cell_id=-9999):
+                 routx=-9999, routy=-9999, name=None, cell_id=-9999):
         '''Defines x and y variables'''
         self.lat = lat
         self.lon = lon
@@ -76,9 +76,9 @@ def __init__(self, param_file, case_name, calendar, out_dir, file_format,
         f = Dataset(param_file, 'r')
         self.n_sources = len(f.dimensions['sources'])
         self.n_outlets = len(f.dimensions['outlets'])
-        self.subset_length = f.variables['subset_length'][0]
-        self.full_time_length = f.variables['full_time_length'][0]
-        self.unit_hydrograph_dt = f.variables['unit_hydrograph_dt'][0]
+        self.subset_length = int(f.variables['subset_length'][:])
+        self.full_time_length = int(f.variables['full_time_length'][:])
+        self.unit_hydrograph_dt = f.variables['unit_hydrograph_dt'][:]
         self.source_lon = f.variables['source_lon'][:]
         self.source_lat = f.variables['source_lat'][:]
         self.source_x_ind = f.variables['source_x_ind'][:]

rvic/parameters.py

@@ -11,6 +11,7 @@
 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
@@ -160,6 +161,7 @@ def gen_uh_init(config_file):
     fdr_lon = config_dict['ROUTING']['LONGITUDE_VAR']
     fdr_vel = config_dict['ROUTING']['VELOCITY']
     fdr_dif = config_dict['ROUTING']['DIFFUSION']
+    fdr_area = config_dict['ROUTING']['SOURCE_AREA_VAR']
     try:
         fdr_data, fdr_vatts, _ = read_netcdf(fdr_file)
         fdr_shape = fdr_data[fdr_var].shape
@@ -169,10 +171,19 @@ def gen_uh_init(config_file):
         if fdr_data[fdr_lat][-1] > fdr_data[fdr_lat][0]:
             log.debug('Flow Direction inputs came in upside down, flipping '
                       'everything now.')
-            var_list = fdr_data.keys()
-            var_list.remove(fdr_lon)
-            for var in var_list:
-                fdr_data[var] = np.flipud(fdr_data[var])
+
+            remove_vars = []
+
+            for var, data in fdr_data.iteritems():
+                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]
         # ---------------------------------------------------------------- #
 
         # ---------------------------------------------------------------- #
@@ -250,6 +261,14 @@ def gen_uh_init(config_file):
                                        glats=fdr_data[fdr_lat],
                                        glons=fdr_data[fdr_lon])
 
+        if options['SEARCH_FOR_CHANNEL']:
+            routys, routxs = search_for_channel(fdr_data[fdr_area],
+                                                routys, routxs, 
+                                                tol=10, search=2)
+            # 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,
@@ -259,7 +278,7 @@ def gen_uh_init(config_file):
         for i in xrange(len(lats)):
             if 'names' in pour_points.keys():
                 name = pour_points['names'].values[i]
-                name = name.replace("'", "").replace(" ", "_")
+                name = name.replace("'", '').replace(" ", "_")
             else:
                 # fill name filed with p-outlet_num
                 name = 'p-{0}'.format(i)