fastscape-lem/fastscape

How to add my DEM data to simulate the evolution process?

Opened this issue · 2 comments

Sorry, I didn't see any example of adding a DEM download from other websites as a simulated object. I don't know how to add my DEM to the program if the fastscape supports doing that. I appreciate that someone could help with that. Thanks.

For example, I have a DEM in the file ("E:/PhD/TopoToolbox/TTLEMStudy/ttlem/topotoolbox-master/DEMdata/WulongshanForTTLEM30.tif"), and I try to use it but into failure. I want to find some examples of this problem, but nothing.

You can first load your DEM file as an xarray DataArray using rioxarray.

The default models in fastscape like basic_model have an "init_topography" process that computes the initial topography (e.g., flat surface with random perturbations, escarpment, etc., see initial conditions).

If you already have an initial topography, you can just drop this process from the model:

from fastscape.models import basic_model


model = basic_model.drop_processes("init_topography")

The newly created model will then expose topography__elevation directly as a model input:

>>> model
<xsimlab.Model (15 processes, 12 inputs)>
grid
    shape           [in] ('shape_yx',) nb. of grid nodes in (y, x)
    length          [in] ('shape_yx',) total grid length in (y, x)
boundary
    status          [in] () or ('border',) node status at borders
fs_context
uplift
    rate            [in] () or ('y', 'x') uplift rate
tectonics
surf2erode
diffusion
    diffusivity     [in] () or ('y', 'x') diffusivity (transport co...
init_erosion
flow
drainage
spl
    k_coef          [in] () or ('y', 'x') bedrock channel incision ...
    area_exp        [in] drainage area exponent
    slope_exp       [in] slope exponent
    tol_rel         [in] relative tolerance (Gauss-Siedel convergence)
    tol_abs         [in] absolute tolerance (Gauss-Siedel convergence)
    max_iter        [in] max nb. of iterations (Gauss-Siedel conver...
erosion
vmotion
topography
    elevation       [inout] ('y', 'x') surface topography elevation
terrain

You can then assign your loaded DEM as simulation input:

import rioxarray
import xsimlab as xs


dem = rioxarray.open_rasterio("dem.tif")

ds_in = xs.create_setup(
    model=model,
    input_vars={
        ...,
        "topography__elevation": dem,
    },
    ...
)

# check_dims="transpose" is only useful if rioxarray does not load the
# DEM with x and y coordinates in the same order
ds_out = ds_in.xsimlab.run(model=model, check_dims="transpose")

Thanks a lot! I will try.