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.