In [1]:
# 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 matplotlib import pyplot as plt
%matplotlib inline
# try to get z
mg,z = read_esri_ascii('./sandwich_bluff_utm.asc',reshape=0,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)
plt.figure()
imshow_node_grid(mg, 'topographic__elevation', cmap='jet', grid_units=['m','m'])
plt.show()
In [2]:
# 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)
xsect = mg.node_vector_to_raster(z, flip_vertically=True)[:,5]
ycoords = mg.node_vector_to_raster(mg.node_y, flip_vertically=True)[:,5]
print np.shape(xsect),xsect
plt.figure()
plt.plot(ycoords, xsect)
plt.show()
In [3]:
# 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.flat),delz
plt.figure() plt.pcolormesh(delz, cmap='RdBu', vmin=0, vmax=5) plt.show()
It seems like this should work too...
...but no. How do we look at maps of gradients?
In [4]:
xdelz = mg.link_vector_to_raster(delz, flip_vertically=True)[:,5]
print np.shape(xdelz),xdelz
In [ ]:
# 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
for i in range(25):
... g = mg.calculate_gradients_at_active_links(z)
... qs = -kd*g
... dqsdx = mg.calculate_flux_divergence_at_nodes(qs)
... dzdt = -dqsdx
... z[interior_nodes] += dzdt[interior_nodes]*dt
Make up a uniform vector field of wind direction U=(u,v) at nodes
In [6]:
# Can I make shorter names?
uxn = 'land_surface_air_flow__x_component_of_velocity'
vxn = 'land_surface_air_flow__y_component_of_velocity'
u = mg.add_zeros('node', uxn)
v = mg.add_zeros('node', vxn)
print('Shape of u {shape}'.format(shape=u.shape))
# Yes
u_at_link = mg.map_mean_of_link_nodes_to_link('land_surface_air_flow__x_component_of_velocity')
v_at_link = mg.map_mean_of_link_nodes_to_link('land_surface_air_flow__y_component_of_velocity')
u_at_link = u_at_link[mg.active_links]
v_at_link = v_at_link[mg.active_links]
print('Shape of u_at_link {shape}'.format(shape=u_at_link.shape))
print('Shape of v_at_link {shape}'.format(shape=v_at_link.shape))
print('Shape of delz {shape}'.format(shape=delz.shape))
In [ ]:
# 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