In [2]:

    

# sandwich wind model
import numpy as np
from landlab import RasterModelGrid

from landlab.components.diffusion.diffusion import LinearDiffuser
from landlab.plot.imshow import imshow_node_grid
from landlab.io.netcdf import *
from landlab.io import read_esri_ascii

from pylab import show, figure
%matplotlib inline

mg,z = read_esri_ascii('./sandwich_zoom.asc',reshape=1,name='topographic__elevation')
mg.axis_name
print('DEM has ' + str(mg.number_of_node_rows) + ' rows, ' +
                  str(mg.number_of_node_columns) + ' columns, and cell size ' + str(mg.dx)) + ' m'
# I could not figure out how to convert -9999.0 to no data, but I think making the water
# close to zero is actuapp better
z[z<-1]=-.0099
print np.shape(z)
mg.set_closed_boundaries_at_grid_edges(True, True, True, True)
figure()
imshow_node_grid(mg, 'topographic__elevation', cmap='jet', grid_units=['m','m'])
show()


    

    




DEM has 317 rows, 351 columns, and cell size 1.05267790787e-05 m
(317, 351)






    




/home/csherwood/miniconda/envs/ioos/lib/python2.7/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str('face'):






    

In [13]:

    

# Add maps of scalar properties at the nodes:

# I did this two ways...both seemed to work, but the shapes of the field arrays are different.

#    Critical shear stress for mobilization
tcrit = 0.05*np.ones_like(z)
print np.shape(tcrit),tcrit
mg.add_field('node', 'tcrit', tcrit, units='Pa', copy=True, noclobber=False)
# will this plot correctly? - yes
figure()
imshow_node_grid(mg, 'tcrit', cmap='jet', grid_units=['m','m'])
show()

# or...
tcrit = mg.add_ones('node','tcrit')
print np.shape(tcrit),tcrit
# will this plot correctly? - yes
figure()
imshow_node_grid(mg, 'tcrit', cmap='jet', grid_units=['m','m'])
show()


    

    




(317, 351) [[ 0.05  0.05  0.05 ...,  0.05  0.05  0.05]
 [ 0.05  0.05  0.05 ...,  0.05  0.05  0.05]
 [ 0.05  0.05  0.05 ...,  0.05  0.05  0.05]
 ..., 
 [ 0.05  0.05  0.05 ...,  0.05  0.05  0.05]
 [ 0.05  0.05  0.05 ...,  0.05  0.05  0.05]
 [ 0.05  0.05  0.05 ...,  0.05  0.05  0.05]]






    













    




(111267,) [ 1.  1.  1. ...,  1.  1.  1.]






    

In [15]:

    

# Calculate the slopes from the DEM delz = (dzdx,dzdy)
# This seems to want z to be a 1D array.
# 1) why is it 2D? (I had to use reshape=1 when importing it)
# 2) Why does it seem ok to do imshow with either shape, but not .calculate_gradients...?
delz = mg.calculate_gradients_at_active_links(z)
print np.shape(delz),delz

# Make up a uniform vector field of wind direction U = (u,v) (at nodes or links?)

# Calculate wind stress (function of U and delz)

# Calculate flux vector Q = (qx,qy)

# Calculate flux divergence


    

    




---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-15-bdfe28e30bb9> in <module>()
      3 # 1) why is it 2D? (I had to use reshape=1 when importing it)
      4 # 2) Why does it seem ok to do imshow with either shape, but not .calculate_gradients...?
----> 5 delz = mg.calculate_gradients_at_active_links(z)
      6 print np.shape(delz),delz
      7 

/home/csherwood/miniconda/envs/ioos/lib/python2.7/site-packages/landlab/testing/decorators.pyc in _wrapped(self, *args, **kwds)
     10             six.print_('Entering: %s.%s' % (self.__class__.__name__,
     11                                             func.__name__))
---> 12         ans = func(self, *args, **kwds)
     13         if self._DEBUG_TRACK_METHODS:
     14             six.print_('Exiting: %s.%s' % (self.__class__.__name__,

/home/csherwood/miniconda/envs/ioos/lib/python2.7/site-packages/landlab/grid/raster.pyc in calculate_gradients_at_active_links(self, node_values, out)
   3508         """
   3509         diffs = gfuncs.calculate_diff_at_active_links(self, node_values,
-> 3510                                                       out=out)
   3511         return np.divide(diffs, self._dx, out=diffs)
   3512 

/home/csherwood/miniconda/envs/ioos/lib/python2.7/site-packages/landlab/grid/grid_funcs.pyc in calculate_diff_at_active_links(grid, node_values, out)
    137     if out is None:
    138         out = grid.empty(centering='active_link')
--> 139     return np.subtract(node_values[grid.activelink_tonode],
    140                        node_values[grid.activelink_fromnode], out=out)
    141 

IndexError: index 703 is out of bounds for axis 0 with size 317

In [ ]:

    

whos