Python FlowRouter.route_flow方法代码示例

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]
def test_degenerate_drainage():
    """
    This "hourglass" configuration should be one of the hardest to correctly
    re-route.
    """
    mg = RasterModelGrid(9,5)
    z_init = mg.node_x.copy()*0.0001 + 1.
    lake_pits = np.array([7, 11, 12, 13, 17, 27, 31, 32, 33, 37])
    z_init[lake_pits] = -1.
    z_init[22] = 0.  # the common spill pt for both lakes
    z_init[21] = 0.1  # an adverse bump in the spillway
    z_init[20] = -0.2  # the spillway
    z = mg.add_field('node', 'topographic__elevation', z_init)

    fr = FlowRouter(mg)
    lf = DepressionFinderAndRouter(mg)
    fr.route_flow()
    lf.map_depressions()

    correct_A = np.array([ 0.,   0.,   0.,   0.,   0.,
                           0.,   1.,   1.,   1.,   0.,
                           0.,   4.,   1.,   3.,   0.,
                           0.,   1.,  10.,   1.,   0.,
                          21.,  21.,   1.,   1.,   0.,
                           0.,   1.,   9.,   1.,   0.,
                           0.,   4.,   1.,   3.,   0.,
                           0.,   1.,   1.,   1.,   0.,
                           0.,   0.,   0.,   0.,   0.])
    
    thelake = np.concatenate((lake_pits, [22])).sort()

    assert_array_almost_equal(mg.at_node['drainage_area'], correct_A)

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]
def test_three_pits():
    """
    A test to ensure the component correctly handles cases where there are
    multiple pits.
    """
    mg = RasterModelGrid(10,10,1.)
    z = mg.add_field('node', 'topographic__elevation', mg.node_x.copy())
    # a sloping plane
    #np.random.seed(seed=0)
    #z += np.random.rand(100)/10000.
    # punch some holes
    z[33] = 1.
    z[43] = 1.
    z[37] = 4.
    z[74:76] = 1.
    fr = FlowRouter(mg)
    lf = DepressionFinderAndRouter(mg)
    fr.route_flow()
    lf.map_depressions()
    
    flow_sinks_target = np.zeros(100, dtype=bool)
    flow_sinks_target[mg.boundary_nodes] = True
    # no internal sinks now:
    assert_array_equal(mg.at_node['flow_sinks'], flow_sinks_target)
    
    # test conservation of mass:
    assert_almost_equal(mg.at_node['drainage_area'
                                       ].reshape((10,10))[1:-1,1].sum(), 8.**2)
    # ^all the core nodes
    
    # test the actual flow field:
    nA = np.array([  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                     8.,   8.,   7.,   6.,   5.,   4.,   3.,   2.,   1.,   0.,
                     2.,   2.,   1.,   1.,   2.,   1.,   1.,   1.,   1.,   0.,
                    26.,  26.,  25.,  15.,  11.,  10.,   9.,   8.,   1.,   0.,
                     2.,   2.,   1.,   9.,   2.,   1.,   1.,   1.,   1.,   0.,
                     2.,   2.,   1.,   1.,   5.,   4.,   3.,   2.,   1.,   0.,
                     2.,   2.,   1.,   1.,   1.,   1.,   3.,   2.,   1.,   0.,
                    20.,  20.,  19.,  18.,  17.,  12.,   3.,   2.,   1.,   0.,
                     2.,   2.,   1.,   1.,   1.,   1.,   3.,   2.,   1.,   0.,
                     0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.])
    assert_array_equal(mg.at_node['drainage_area'], nA)
    
    #test a couple more properties:
    lc = np.empty(100, dtype=int)
    lc.fill(XX)
    lc[33] = 33
    lc[43] = 33
    lc[37] = 37
    lc[74:76] = 74
    assert_array_equal(lf.lake_map, lc)
    assert_array_equal(lf.lake_codes, [33, 37, 74])
    assert_equal(lf.number_of_lakes, 3)
    assert_array_almost_equal(lf.lake_areas, [2., 1., 2.])
    assert_array_almost_equal(lf.lake_volumes, [2., 2., 4.])

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]
def test_internal_closed():
    """
    Test closed nodes in the core of the grid.
    """
    fr = FlowRouter(mg)
    fr.route_flow()
    assert_array_almost_equal(A_target, mg.at_node['drainage_area'])
    assert_array_equal(frcvr_target, mg.at_node['flow_receiver'])
    assert_array_equal(links2rcvr_target, mg.at_node['links_to_flow_receiver'])
    assert_array_almost_equal(A_target, mg.at_node['water__volume_flux'])
    assert_array_almost_equal(steepest_target,
                              mg.at_node['topographic__steepest_slope'])

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]
def test_irreg_topo():
    """
    Tests D4 routing on a toy irregular topo.
    """
    fr = FlowRouter(mg)
    fr.route_flow(method='D4')
    assert_array_equal(A_target_D4, mg.at_node['drainage_area'])
    assert_array_equal(frcvr_target_D4, mg.at_node['flow_receiver'])
    assert_array_equal(upids_target_D4, mg.at_node['upstream_node_order'])
    assert_array_equal(links2rcvr_target_D4,
                       mg.at_node['links_to_flow_receiver'])
    assert_array_almost_equal(steepest_target_D4,
                              mg.at_node['topographic__steepest_slope'])

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]
def test_accumulate_D8():
    """
    Test accumulation works for D8 in a simple scenario
    """
    fr = FlowRouter(mg)
    fr.route_flow()
    assert_array_equal(A_target, mg.at_node['drainage_area'])
    assert_array_equal(frcvr_target, mg.at_node['flow_receiver'])
    assert_array_equal(upids_target, mg.at_node['upstream_node_order'])
    assert_array_equal(links2rcvr_target, mg.at_node['links_to_flow_receiver'])
    assert_array_equal(A_target, mg.at_node['water__volume_flux'])
    assert_array_equal(steepest_target,
                       mg.at_node['topographic__steepest_slope'])

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]
def test_variable_Qin():
    """
    Tests a variable Qin field.
    """
    Qin_local = np.zeros(25, dtype=float)
    Qin_local[13] = 2.
    mg.add_field('node', 'water__volume_flux_in',
                 Qin_local, units='m**3/s')
    fr = FlowRouter(mg)
    fr.route_flow()
    Qout_local = np.zeros_like(Qin_local)
    Qout_local[10:14] = 200.
    assert_array_equal(Qout_local, mg.at_node['water__volume_flux'])
    assert_array_equal(A_target, mg.at_node['drainage_area'])

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]

#.........这里部分代码省略.........
        apply_slope : bool
            Whether to leave the filled surface flat (default), or apply a
            gradient downwards through all lake nodes towards the outlet.
            A test is performed to ensure applying this slope will not alter
            the drainage structure at the edge of the filled region (i.e.,
            that we are not accidentally reversing the flow direction far
            from the outlet.)

        Return fields
        -------------
        'topographic__elevation' : the updated elevations
        'sediment_fill__depth' : the depth of sediment added at each node
        """
        self.original_elev = self._elev.copy()
        # We need this, as we'll have to do ALL this again if we manage
        # to jack the elevs too high in one of the "subsidiary" lakes.
        # We're going to implement the lake_mapper component to do the heavy
        # lifting here, then delete its fields. This means we first need to
        # test if these fields already exist, in which case, we should *not*
        # delete them!
        existing_fields = {}
        spurious_fields = set()
        set_of_outputs = self._lf.output_var_names | self._fr.output_var_names
        try:
            set_of_outputs.remove(self.topo_field_name)
        except KeyError:
            pass
        for field in set_of_outputs:
            try:
                existing_fields[field] = self._grid.at_node[field].copy()
            except FieldError:  # not there; good!
                spurious_fields.add(field)

        self._fr.route_flow(method=self._routing)
        self._lf.map_depressions(pits=self._grid.at_node['flow_sinks'],
                                 reroute_flow=True)
        # add the depression depths to get up to flat:
        self._elev += self._grid.at_node['depression__depth']
        # if apply_slope is none, we're now done! But if not...
        if apply_slope:
            # new way of doing this - use the upstream structure! Should be
            # both more general and more efficient
            for (outlet_node, lake_code) in zip(self._lf.lake_outlets,
                                                self._lf.lake_codes):
                lake_nodes = np.where(self._lf.lake_map == lake_code)[0]
                lake_perim = self.get_lake_ext_margin(lake_nodes)
                perim_elevs = self._elev[lake_perim]
                out_elev = self._elev[outlet_node]
                lowest_elev_perim = perim_elevs[perim_elevs != out_elev].min()
                # note we exclude the outlet node
                elev_increment = ((lowest_elev_perim-self._elev[outlet_node]) /
                                  (lake_nodes.size + 2.))
                assert elev_increment > 0.
                all_ordering = self._grid.at_node['upstream_node_order']
                upstream_order_bool = np.in1d(all_ordering, lake_nodes,
                                              assume_unique=True)
                lake_upstream_order = all_ordering[upstream_order_bool]
                argsort_lake = np.argsort(lake_upstream_order)
                elevs_to_add = (np.arange(lake_nodes.size, dtype=float) +
                                1.) * elev_increment
                sorted_elevs_to_add = elevs_to_add[argsort_lake]
                self._elev[lake_nodes] += sorted_elevs_to_add
        # now put back any fields that were present initially, and wipe the
        # rest:
        for delete_me in spurious_fields:
            self._grid.delete_field('node', delete_me)

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]
print 'fixed vals in grid: ', mg.fixed_value_node_properties['values']

# Display a message
print 'Running ...' 

#instantiate the components:
fr = FlowRouter(mg)
sde = SedDepEroder(mg, input_file)
vid = VideoPlotter(mg, data_centering='node')

time_on = time()
#perform the loops:
for i in xrange(nt):
    #print 'loop ', i
    mg.at_node['topographic_elevation'][mg.core_nodes] += uplift_per_step
    mg = fr.route_flow(grid=mg)
    #mg.calculate_gradient_across_cell_faces(mg.at_node['topographic_elevation'])
    #neighbor_slopes = mg.calculate_gradient_along_node_links(mg.at_node['topographic_elevation'])
    #mean_slope = np.mean(np.fabs(neighbor_slopes),axis=1)
    #max_slope = np.max(np.fabs(neighbor_slopes),axis=1)
    #mg,_,capacity_out = tl.erode(mg,dt,slopes_at_nodes='steepest_slope')
    #mg,_,capacity_out = tl.erode(mg,dt,slopes_at_nodes=max_slope)
    mg_copy = deepcopy(mg)
    mg,_ = sde.erode(mg,dt)
    #print sde.iterations_in_dt
    #print 'capacity ', np.amax(capacity_out[mg.core_nodes])
    #print 'rel sed ', np.nanmax(sed_in[mg.core_nodes]/capacity_out[mg.core_nodes])
    if i%100 == 0:
        print 'loop ', i
        print 'max_slope', np.amax(mg.at_node['steepest_slope'][mg.core_nodes])
        pylab.figure("long_profiles")

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]
time_on = time.time()

#instantiate the components:
fr = FlowRouter(mg)
sp = StreamPowerEroder(mg, './drive_sp_params.txt')
#load the Fastscape module too, to allow direct comparison
fsp = Fsc(mg, './drive_sp_params.txt')

#perform the loop:
elapsed_time = 0. #total time in simulation
while elapsed_time < time_to_run:
    print(elapsed_time)
    if elapsed_time+dt>time_to_run:
        print("Short step!")
        dt = time_to_run - elapsed_time
    mg = fr.route_flow()
    #print 'Area: ', numpy.max(mg.at_node['drainage_area'])
    #mg = fsp.erode(mg)
    mg,_,_ = sp.erode(mg, dt, node_drainage_areas='drainage_area', slopes_at_nodes='topographic__steepest_slope', K_if_used='K_values')
    #add uplift
    mg.at_node['topographic__elevation'][mg.core_nodes] += uplift*dt
    elapsed_time += dt

time_off = time.time()
print('Elapsed time: ', time_off-time_on)

#Finalize and plot
elev = mg['node']['topographic__elevation']
elev_r = mg.node_vector_to_raster(elev)

# Clear previous plots

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]
grid.set_inactive_boundaries(True, False, True, True)
z = grid.add_zeros('node', 'Land_Surface__Elevation')
z[6] = 4.5
z[7] = 3.
z[8] = 1.
z[11] = 4.
z[12] = 2.8
z[13] = 2.


# Get array of interior (active) node IDs
interior_nodes = grid.get_active_cell_node_ids()

# Route flow
flow_router = FlowRouter(grid)
grid = flow_router.route_flow()

for i in range(grid.number_of_nodes):
    print(i, grid.node_x[i], grid.node_y[i], z[i], grid.node_status[i], \
          r[i], a[i], q[i], ss[i], rl[i])

# Let's take a look for debugging
#print 'node  receiver  flow_link'
#for i in interior_nodes:
#    print i, r[i], rl[i]

# Calculate lengths of flow links
flow_link_length = ones(size(z))
flow_link_length[interior_nodes] = grid.link_length[rl[interior_nodes]] #DEJH suspects a node ordering bug here - rl is not in ID order, but interior_nodes is
print('fll:', flow_link_length)

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]

#.........这里部分代码省略.........
    
        #interior_nodes are the nodes on which you will be calculating flow 
        self.interior_nodes = grid.get_active_cell_node_ids()

        if total_t==None:
            total_t = self.rain_duration
        if rainrate==None:
            rainrate = self.rainfall_rate
        if rainduration==None:
            rainduration = total_t
            
        elapsed_time = 0
        
        #below is for calculating water surface slope at interior nodes
        w_slope=np.zeros(self.interior_nodes.size)

        t = [] #time array for plotting
        t.append(0.) #initialize array
        self.dqds= grid.zeros(centering='node')
                
        while elapsed_time < total_t:
            
            # Calculate time-step size for this iteration (Bates et al., eq 14)
            dtmax = self.alpha*grid.dx/np.sqrt(self.g*np.amax(self.h))
        
            # Take the smaller of delt or calculated time-step
            dt = min(dtmax, total_t)
        
            # Calculate the effective flow depth at active links. Bates et al. 2010
            # recommend using the difference between the highest water-surface
            # and the highest bed elevation between each pair of cells.
            zmax = grid.max_of_link_end_node_values(z) #array of length of active links
            w = self.h+z   # water-surface height, array of length num nodes
            wmax = grid.max_of_link_end_node_values(w) #array of length of active links
            hflow = wmax - zmax #array of length of active links
        
            # Calculate water-surface slopes: across links, but water heights are 
            #defined at nodes.
            water_surface_slope = grid.calculate_gradients_at_active_links(w)
       
      
            # Calculate the unit discharges (Bates et al., eq 11)
            self.q = (self.q-g*hflow*dtmax*water_surface_slope)/ \
                (1.+g*hflow*dtmax*0.06*0.06*abs(self.q)/(hflow**ten_thirds))                
            #NOTES:    
            # q is calculated at links
            # water_surface_slope is at links
            # hflow is at links, but w is at nodes

            # Calculate water-flux divergence at nodes
            self.dqds = grid.calculate_flux_divergence_at_nodes(self.q,self.dqds)
            
            # Calculate rate of change of water depth
            self.dhdt = rainrate-self.dqds

            
            # Second time-step limiter (experimental): make sure you don't allow
            # water-depth to go negative
            if np.amin(self.dhdt) < 0.:
                shallowing_locations = np.where(self.dhdt<0.)
                time_to_drain = -self.h[shallowing_locations]/self.dhdt[shallowing_locations]
                dtmax2 = self.alpha*np.amin(time_to_drain)
                dt = np.min([dtmax, dtmax2, total_t])
        
            # Update the water-depth field
            self.h[self.interior_nodes] = self.h[self.interior_nodes] + self.dhdt[self.interior_nodes]*dt


            # Let's calculate shear stress at the nodes.  
            # First get water height at the nodes.
            # Then calculate the maximum gradient in water surface elevations (S).
            # Then you can calculate shear stress! (rho g h S)       
            # h (water depth) is at nodes
            
            w = self.h+z   # water-surface height, array of length num nodes
            

            #Below if is for limiting shear stress calculations to only times when
            #q surpasses a threshold (in this case q should be in m^3/sec)
            #Note that this threshold is HARDWIRED below (on right of >)

            #self.slopes_at_node, garbage = grid.calculate_steepest_descent_on_nodes(w, water_surface_slope)
            fr=FlowRouter(grid)
            r, a, q, ss, s, d = fr.route_flow(w)
            self.slopes_at_node = ss
            if self.q.any()*grid.dx> 0.2:
                self.tau = self.rho*self.g*self.slopes_at_node*self.h

                #for i in range(24899): 
                #    #w_slope[i],garbage=grid.calculate_gradient_across_cell_faces(grid, w, self.interior_nodes)    #grid.calculate_max_gradient_across_node_d4(w,self.interior_nodes[i])
                #    self.tau[self.interior_nodes[i]]=self.rho*self.g*self.slopes_at_node[self.interior_nodes[i]]*self.h[self.interior_nodes[i]]
                #    #self.tau[i] = tau[i]

            self.tau[np.where(self.tau<0)] = 0       

            #self.tau[interior_nodes] = self.tau[interior_nodes] + self.dtds[interior_nodes]*dt
            #tau[np.where(tau<0)] = 0
            # Update model run time
            elapsed_time += dt
            print "elapsed time", elapsed_time

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]
              5., 3., 2., 3., 5.,
              5., 4., 4., 4., 5.,
              5., 5., 5., 5., 5.])
mg['node']['topographic__elevation'] = z

print('Running ...')

# instantiate the components:
fr = FlowRouter(mg)
sp = StreamPowerEroder(mg, './drive_sp_params_discharge.txt')
# load the Fastscape module too, to allow direct comparison
fsp = Fsc(mg, './drive_sp_params_discharge.txt')

# perform the loop (once!)
for i in xrange(1):
    fr.route_flow(method='D8')
    my_Q = mg.at_node['water__volume_flux']*1.
    sp.erode(mg, dt, node_drainage_areas='drainage_area',
             slopes_at_nodes='topographic__steepest_slope',
             Q_if_used=my_Q)
    # no uplift

# print the stream power that was calculated:
print('stream power values:')
print(mg.at_node['stream_power_erosion'])

# Finalize and plot
elev = mg['node']['topographic__elevation']
elev_r = mg.node_vector_to_raster(elev)

# Clear previous plots

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]
def test_composite_pits():
    """
    A test to ensure the component correctly handles cases where there are
    multiple pits, inset into each other.
    """
    mg = RasterModelGrid(10, 10, 1.)
    z = mg.add_field('node', 'topographic__elevation', mg.node_x.copy())
    # a sloping plane
    #np.random.seed(seed=0)
    #z += np.random.rand(100)/10000.
    # punch one big hole
    z.reshape((10,10))[3:8,3:8] = 0.
    # dig a couple of inset holes
    z[57] = -1.
    z[44] = -2.
    z[54] = -10.
    fr = FlowRouter(mg)
    lf = DepressionFinderAndRouter(mg)
    fr.route_flow()
    lf.map_depressions()
    
    flow_sinks_target = np.zeros(100, dtype=bool)
    flow_sinks_target[mg.boundary_nodes] = True
    # no internal sinks now:
    assert_array_equal(mg.at_node['flow_sinks'], flow_sinks_target)
    
    # test conservation of mass:
    assert_almost_equal(mg.at_node['drainage_area'
                                       ].reshape((10,10))[1:-1,1].sum(), 8.**2)
    # ^all the core nodes
    
    # test the actual flow field:
    nA = np.array([  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                     8.,   8.,   7.,   6.,   5.,   4.,   3.,   2.,   1.,   0.,
                     1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.,   0.,
                     1.,   1.,   1.,   4.,   2.,   2.,   8.,   4.,   1.,   0.,
                     1.,   1.,   1.,   8.,   3.,  15.,   3.,   2.,   1.,   0.,
                     1.,   1.,   1.,  13.,  25.,   6.,   3.,   2.,   1.,   0.,
                     1.,   1.,   1.,  45.,   3.,   3.,   5.,   2.,   1.,   0.,
                    50.,  50.,  49.,   3.,   2.,   2.,   2.,   4.,   1.,   0.,
                     1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.,   0.,
                     0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.])
    assert_array_equal(mg.at_node['drainage_area'], nA)
    
    # the lake code map:
    lc = np.array([XX, XX, XX, XX, XX, XX, XX, XX, XX, XX,
                   XX, XX, XX, XX, XX, XX, XX, XX, XX, XX,
                   XX, XX, XX, XX, XX, XX, XX, XX, XX, XX,
                   XX, XX, XX, 57, 57, 57, 57, 57, XX, XX,
                   XX, XX, XX, 57, 57, 57, 57, 57, XX, XX,
                   XX, XX, XX, 57, 57, 57, 57, 57, XX, XX,
                   XX, XX, XX, 57, 57, 57, 57, 57, XX, XX,
                   XX, XX, XX, 57, 57, 57, 57, 57, XX, XX,
                   XX, XX, XX, XX, XX, XX, XX, XX, XX, XX,
                   XX, XX, XX, XX, XX, XX, XX, XX, XX, XX])
    
    #test the remaining properties:
    assert_equal(lf.lake_outlets.size, 1)
    assert_equal(lf.lake_outlets[0], 72)
    outlets_in_map = np.unique(lf.depression_outlet_map)
    assert_equal(outlets_in_map.size, 2)
    assert_equal(outlets_in_map[0], 72)
    assert_equal(lf.number_of_lakes, 1)
    assert_equal(lf.lake_codes[0], 57)
    assert_array_equal(lf.lake_map, lc)
    assert_almost_equal(lf.lake_areas[0], 25.)
    assert_almost_equal(lf.lake_volumes[0], 63.)

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]
#load the Fastscape module too, to allow direct comparison
fsp = Fsc(mg, input_file_string)

try:
    #raise NameError
    mg = copy.deepcopy(mg_mature)
except NameError:
    print('building a new grid...')
    out_interval = 50000.
    last_trunc = time_to_run #we use this to trigger taking an output plot
    #run to a steady state:
    #We're going to cheat by running Fastscape SP for the first part of the solution
    for (interval_duration, rainfall_rate) in precip.yield_storm_interstorm_duration_intensity():
        if rainfall_rate != 0.:
            mg.at_node['water__volume_flux_in'].fill(rainfall_rate)
            mg = fr.route_flow()
            #print 'Area: ', numpy.max(mg.at_node['drainage_area'])
            mg,_,_ = sp.erode(mg, interval_duration, Q_if_used='water__volume_flux', K_if_used='K_values')
        #add uplift
        mg.at_node['topographic__elevation'][mg.core_nodes] += uplift*interval_duration
        this_trunc = precip.elapsed_time//out_interval
        if this_trunc != last_trunc: #a new loop round
            print('made it to loop ', out_interval*this_trunc)
            last_trunc=this_trunc

    mg_mature = copy.deepcopy(mg)

else:
    #reinstantiate the components with the new grid
    fr = FlowRouter(mg)
    sp = StreamPowerEroder(mg, input_file_string)

# 需要导入模块: from landlab.components.flow_routing.route_flow_dn import FlowRouter [as 别名]
# 或者: from landlab.components.flow_routing.route_flow_dn.FlowRouter import route_flow [as 别名]
from landlab import VoronoiDelaunayGrid  # , RasterModelGrid
from landlab.components.flow_routing.route_flow_dn import FlowRouter
from landlab.components.stream_power.stream_power import StreamPowerEroder
from landlab.plot.imshow import imshow_node_grid
import numpy as np
from matplotlib.pyplot import figure, show

nnodes = 10000

x, y = np.random.rand(nnodes), np.random.rand(nnodes)
mg = VoronoiDelaunayGrid(x,y)
#mg = RasterModelGrid(4,5)

z = mg.add_field('node', 'topographic__elevation', np.random.rand(nnodes)/10000., copy=False)

fr = FlowRouter(mg)
spe = StreamPowerEroder(mg, 'drive_sp_params_voronoi.txt')

for i in xrange(100):
    z[mg.core_nodes] += 0.01
    fr.route_flow()
    spe.erode(mg, 1.)

imshow_node_grid(mg, 'topographic__elevation')

show()