Meshing from a Digital Elevation Model (DEM) .adf file

Estimated time to run the notebook = 5min

# !! run preprocessor change the DEM shape !
# dtm_13 does not have the same shape anymore!

import os

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import pyCATHY.meshtools as mt
from pyCATHY import cathy_tools
from pyCATHY.importers import cathy_inputs as in_CT
from pyCATHY.importers import cathy_outputs as out_CT
from pyCATHY.plotters import cathy_plots as cplt

import rioxarray
import pyvista as pv

---

path2prj = "../SSHydro/"  # add your local path here
simu = cathy_tools.CATHY(dirName=path2prj,
                         prj_name="meshing_from_subcachment"
                         )

rootpath = os.path.join(simu.workdir + simu.project_name)


# Path to the directory containing the .adf file (not the file itself)
adf_folder = "../../data/dtmplot1/"

# Open the raster (typically named 'hdr.adf', but you only need the folder)
raster_DEM = rioxarray.open_rasterio(adf_folder, masked=True).isel(band=0)

# Create a mask of valid (non-NaN) data
valid_mask = ~np.isnan(raster_DEM)

# Apply the mask to get valid coordinates
valid_x = raster_DEM['x'].where(valid_mask.any(dim='y'), drop=True)
valid_y = raster_DEM['y'].where(valid_mask.any(dim='x'), drop=True)

# Get min and max valid coordinates
min_lon, max_lon = float(valid_x.min()), float(valid_x.max())
min_lat, max_lat = float(valid_y.min()), float(valid_y.max())

raster_DEM_masked = raster_DEM.where(
    (raster_DEM['x'] >= min_lon) & (raster_DEM['x'] <= max_lon), drop=True
).where(
    (raster_DEM['y'] >= min_lat) & (raster_DEM['y'] <= max_lat), drop=True
)

raster_DEM_masked = np.where(np.isnan(raster_DEM_masked), -9999, raster_DEM_masked)
np.shape(raster_DEM_masked)
np.shape(raster_DEM)

🏁 Initiate CATHY object

(161, 283)

fig, ax = plt.subplots(1)
img = ax.imshow(raster_DEM_masked)
plt.colorbar(img)

simu.show_input(prop="dem")

simu.update_prepo_inputs(
    DEM=raster_DEM_masked,
    delta_x=5,
    delta_y=5,
    ivert=1,
)

fig = plt.figure()
ax = plt.axes(projection="3d")
simu.show_input(prop="dem", ax=ax)
simu.create_mesh_vtk(verbose=True)

• 
• 
• 

🔄 Update hap.in file
🔄 Update dem_parameters file
🔄 Update dem_parameters file
🔄 Update hap.in file
🔄 Update dem_parameters file
🔄 Update dtm_13 file
🔄 update zone file
🔄 Update dem_parameters file
🔄 Update parm file
🔄 Update dem_parameters file
🍳 gfortran compilation
👟 Run preprocessor

 wbb...

 searching the dtm_13.val input file...
 assigned nodata value =  -9999.0000000000000

 number of processed cells =        1174

 ...wbb completed

 rn...
 csort I...
 ...completed

 depit...
 dem modifications =            0
 dem modifications =            0 (total)
 ...completed

 csort II...
 ...completed

 cca...

 contour curvature threshold value =    9.99999996E+11
 ...completed

 smean...
 mean (min,max) facet slope =  0.216937334 ( 0.012000000, 0.760631317)
 ...completed

 dsf...
 weight w_1 less than 1.0E-06: set w_1 = 0 and w_2 = 1
 weight w_1 less than 1.0E-06: set w_1 = 0 and w_2 = 1
 the drainage direction of the outlet cell (           2 ) is used
 ...completed

 hg...
 ...completed

 saving the data in the basin_b/basin_i files...

 ...rn completed

 mrbb...


 Select the header type:
 0) None
 1) ESRI ascii file
 2) GRASS ascii file
 (Ctrl C to exit)

 ->
 Select the nodata value:
 (Ctrl C to exit)

 ->
 Select the pointer system:
 1) HAP system
 2) Arc/Gis system
 (Ctrl C to exit)

 ->  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dem file

 min value = 0.567960E+03
 max value = 0.628820E+03
 number of cells =  1174
 mean value = 0.592434E+03

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 lakes_map file

 min value =     0
 max value =     0
 number of cells =  1174
 mean value =     0.000000

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 zone file

 min value =     1
 max value =     1
 number of cells =  1174
 mean value =     1.000000

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_w_1 file

 min value = 0.000000E+00
 max value = 0.100000E+01
 number of cells =  1174
 mean value = 0.558963E+00

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_w_2 file

 min value = 0.000000E+00
 max value = 0.100000E+01
 number of cells =  1174
 mean value = 0.441037E+00

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_p_outflow_1 file

 min value =     2
 max value =     6
 number of cells =  1174
 mean value =     2.272573

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_p_outflow_2 file

 min value =     0
 max value =     9
 number of cells =  1174
 mean value =     1.447189

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 A_inflow file

 min value = 0.000000000000E+00
 max value = 0.293249992721E+05
 number of cells =  1174
 mean value = 0.770825256348E+03

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_local_slope_1 file

 min value =-0.234000E+00
 max value = 0.750000E+00
 number of cells =  1174
 mean value = 0.189089E+00

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_local_slope_2 file

 min value =-0.189505E+00
 max value = 0.692965E+00
 number of cells =  1174
 mean value = 0.173164E+00

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_epl_1 file

 min value = 0.000000E+00
 max value = 0.500000E+01
 number of cells =  1174
 mean value = 0.469762E+01

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_epl_2 file

 min value = 0.000000E+00
 max value = 0.707107E+01
 number of cells =  1174
 mean value = 0.637239E+01

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_kSs1_sf_1 file

 min value = 0.000000E+00
 max value = 0.665000E+02
 number of cells =  1174
 mean value = 0.288060E+02

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_kSs1_sf_2 file

 min value = 0.000000E+00
 max value = 0.665000E+02
 number of cells =  1174
 mean value = 0.274720E+02

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_Ws1_sf file

 min value = 0.000000E+00
 max value = 0.952684E+01
 number of cells =  1174
 mean value = 0.141419E+01

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_Ws1_sf_2 file

 min value = 0.000000E+00
 max value = 0.667605E+01
 number of cells =  1174
 mean value = 0.126586E+01

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_b1_sf file

 min value = 0.000000E+00
 max value = 0.260000E+00
 number of cells =  1174
 mean value = 0.538160E-01

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_y1_sf file

 min value = 0.000000E+00
 max value = 0.000000E+00
 number of cells =  1174
 mean value = 0.000000E+00

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_hcID file

 min value =     0
 max value =     1
 number of cells =  1174
 mean value =     0.206985

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_q_output file

 min value =     0
 max value =     0
 number of cells =  1174
 mean value =     0.000000

 writing the output file...

 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 dtm_nrc file

 min value = 0.100000E+01
 max value = 0.100000E+02
 number of cells =  1174
 mean value = 0.813714E+01

 writing the output file...

 ...mrbb completed

 bb2shp...

 writing file river_net.shp

Note: The following floating-point exceptions are signalling:
IEEE_UNDERFLOW_FLAG IEEE_DENORMAL

🔄 Update parm file
🛠  Recompile src files [7s]
🍳 gfortran compilation [16s]
❌ Linker error:
/usr/bin/ld: cannot find -llapack: No such file or directory
/usr/bin/ld: cannot find -lblas: No such file or directory
collect2: error: ld returned 1 exit status

👟 Run processor
b''

meshfile = rootpath + “/vtk/” + simu.project_name + “.vtk”

# mesh2plot = pv.read(meshfile)

# pl = pv.Plotter(off_screen=True)

# pl.add_mesh(mesh2plot)
# pl.show_bounds()

# pl.show()
# mesh2plot.plot(show_edges=True,
#                show_axes=True,
#                show_bounds=True
#                )

simu.run_preprocessor(verbose=False)
# simu.run_processor(IPRT1=3,verbose=True)

# # simu.read_inputs('atmbc')
simu.update_parm(TIMPRTi=[1800,7200])

grid3d = simu.read_outputs('grid3d')

t_atmbc = [0,86400,86400*4]
v_atmbc_t_rain = 1e-7
netValue = np.zeros(len(t_atmbc))
netValue[0] = v_atmbc_t_rain
# fig, ax = plt.subplots()
# ax.imshow(v_atmbc)


#% Create an empty dataframe of SPP and set default SPP properties
df_SPP_map = simu.init_soil_SPP_map_df(nzones=1,nstr=15)
SPP_map = simu.set_SOIL_defaults(SPP_map_default=True)


SPP_map['PERMX'] = 6.88e-4
SPP_map['PERMY'] = 6.88e-4
SPP_map['PERMZ'] = 6.88e-4

#% Update soil file
simu.update_soil(SPP_map=SPP_map)


simu.update_ic(INDP=0,
               IPOND=0,
               pressure_head_ini=-15
                  )


simu.update_atmbc(
                    HSPATM=1,
                    IETO=1,
                    time=t_atmbc,
                    netValue=netValue
                  )

simu.update_parm(
                        IPRT=4,
                        VTKF=2, # dont write vtk files
                        )

🍳 gfortran compilation
👟 Run preprocessor
🔄 Update parm file
🔄 Update soil
homogeneous soil
🔄 Update ic
Single value detected for ic ==> assumming it homogeneous
🔄 Update atmbc
🔄 Update parm file
🔄 Update parm file

int(grid3d['nnod'])


# simu.run_processor(IPRT1=2,
#                     DTMIN=1e-2,
#                     DTMAX=1e3,
#                     DELTAT=1e2,
#                     TRAFLAG=0,
#                     verbose=True
#                     )

1271

# # pl = pv.Plotter(notebook=True)
# pl = pv.Plotter(off_screen=True)



# cplt.show_vtk(unit="pressure",
#               timeStep=1,
#               path=simu.workdir + "/meshing_from_subcachment/vtk/",
#               ax=pl,
#               )
# # pl.show()
# image = pl.screenshot('test.png')