open TELEMAC-MASCARET The mathematically superior suite of solvers

Dear forum,

I've been struggling with a rather shady topic about post-processing a flood map (for a risk-related map) at the office, mostly related to the spatial resolution of the flood-extent (generated using BK which is quite similar to Crayfish and post-telemac ones) and the provided High-res DEM which was used to generate the "BOTTOM" variable of the model.

The issue is that the flood-map extent has some areas that shouldn't be marked as flooded and some that should be (that is, the shape needs to be extended) according to the elevation in those areas, but not because of a numerical inaccuracy of TELEMAC, but because of the fact that the mesh generated , which is quite rather fine, will never have the same resolution as the original high-res DEM and the interpolated nodes have some degree of error too. I strongly believe this is a common problem in all hydraulic models so its understandable...isn't it?

Except that, in my country, everyone is used to see HEC-RAS "post-processed"-like flood extent map that tries to emulate a continuous flood extent and getting a really "detailed" output from a 1D model.

The HEC-RAS (RAS mapper in fact) infamous type of flood extent map:




The list of problems that I can think of that specific "correcting" process in a 1D model are many and I'm sure there are more.

Now, I do (following a proper evangelizing process of TELEMAC in the office) have to "correct" all those areas (i.e. extend or contract) according to the high-res DEM that was provided to match only the "detail" in the extent of the flood, that is, that it looks nicer and matches the high-res DEM (from which the BOTTOM variable comes from).

So, has anyone had this issue?. Did anybody find a good solution?

I don't want to re-invent the wheel, but if it needed I'm more than happy that hear your thoughts on a possible approach and code a solution based on them and post it here.

PD: Just wondering, how much different are the results of a fancy hec-ras flood extent map compared to a robust 2D model of Telemac in a "real-life" not the "laboratory" case...


Regards,

José Díaz.

Hello,

he best way would probably to interpolate the TELEMAC results (SELAFIN format) onto your DEM (other format). If you can work with python scripts, there is a function in ./scripts/python27/samplers/meshes.py called mapThinPlateSpline, which you can use as follows:


where npoin can be set to 30, alpha to 0.2 and where map is your TELEMAC result and slf.MESHX,slf.MESHY are your DEM (x,y) coordinate. As it happens, slf is also a SELAFIN, but you can adapt this to your own format parser.
Also, map.getVariablesAt(0,[0])[0] is the extracted values from the TELEMAC results for time 0 (the first time frame in the file) and variable 0 (the first variable in the file), both of which you can change of even loop over.

Finally, you may apply a filter or a mask on VARSOR, knowing that if the DEM value is higher than the (interpolated) water level, then you keep your DEM value (or leave those value as 0 (or NaN if your format supports those).

Hope this help,
Sébastien.

Dear Sébastien,

I'm sorry for this late reply, I couldn't find the function "mapThin..." in the meshes module or even in the other files in the python scripts folder (v7p1r1 release), is this perhaps something you added to the trunk?

Thank you for your time!

Regards,

José D.

EDIT:

Gave myself the answer without checking first, it *was* on the trunk version of the python scripts and not in the release one.

I'll let you know how this goes with that specific function!

Thank you!

Dear Sébastien and if somehow anyone is related to this...

I'm starting to touch an issue with Telemac's free surface variable output, which a near the (invisible in Freesurface variable)"boundary" of wet and dry elements seems to behave a little bit unexpected.








the attachment shows a test cross section where the free surface (interpolated to a 2m spaced points) "bends" to match the elevation of a node (which is the edge between "dry" nodes and "wet" ones) even when I explicitly extract all nodes which WaterDepth is higher or equal to 10cm (threshold)

This isn't really an issue when low flows are simulated and extracting the flood map with bluekenue or other tool (base flow as the pictures) but it really is when analyzing a Dam break study: The flood map generated here includes those really high nodes making the maximum contour level of water (relating to the DEM) skyrocket (sometimes I expect a 36-37m level of water and get 45-60m of elevation in the edges of the flood map polygon generated)

I've been reading also about "THRESHOLD FOR RECEDING PROCEDURE..:" (set as high as 15cm) and that doesn't solve the problem, and also re-meshing is imposible (its a 700k nodes on a really long river whose computation times are at limit of acceptable for sensibility cases)

Essentially, this "polish" procedure of flood map has become more and more important as time passes, and I refuse to believe I'm the only one that has happened this for a large scale model because is near imposible to model such large domain with a really-fine mesh...

So,before panicking, the routine that was thinking to implement to extract a proper flood map is:

1. Read the values of the results: Water depth (WD), free surface (FS)
2. first filter: Discards nodes of FS that has a water depth below a threshold (lets say 5cm),discarding both "land nodes" and small "ponds" getting as results "only water bodies nodes"
3. second filter:Discard the nodes that has a high slope(?) or gradient maybe(?) or are outside (outliers, both high and low) a percentile range (say 20-80 % of surrounding data of near "N" points ) which are most probably the nodes pointed in the above img?.
4. With the resulted FS nodes apply an interpolation (thinPlateSplines as you suggested) to every cell (DEM) that lies within the mesh
5. Third filter: remove (make NoData or NaN) or mask every value that mathes the following condition: Nodes<HighResDEM+tol
6. Write the result output to two raster: FS Values (for debug time) and a boleean one (water=1 and no water=0)
7. Convert the boleean raster to polygon Shapefile (and remove the polygon value of 0)

Is there a way to mask nodes in the SELAFIN class? (I'm thinkin of using numpy masked arrays but dont want to step outside of the selafin class too much).

But the big issue here exactly having trouble with the 3rd step of that procedure! dont know how to implement that second filter (identify and remove those nasty nodes) before the interpolation of the free surface...

is that a decent procedure? would you do it differently?

Regards,

José D.

Hi José

It's a well known problem.
There is no automatic solution to deal with this result, just potential solution or ideas to try to solve this problem...
The main problem is link to the mesh structure which not allow you to have a precise location of the 0 water depth.

Personally, I started to develop a tool which works on Shapefile to post-process the results. The idea is to identify elements where at least one node is dry and from the wet node value of free surface, try to cut the element at the location where the free surface = bottom level.
This solution works quite well where the bottom gradients are positive but doesn't work when you propagate water wave on a dry flat bottom...

Your solution is also an option. This gives me another idea:

• extract a submesh with only positive water depth
• from this mesh and the global mesh, build a second mesh with only the first dry node linked to the first extracted mesh
• gives to the dry nodes the free surface value of the nearest wet node
• compute the intersection of this free surface plane with the DEM taken from the lidar

Hope this helps

Regards

Hi Christophe and all,

Your "idea" is what I have been using on ArcGIS based on a scripted toolbox. It works fine but this method has some limitations since certain dry areas of the mesh can be mapped as flooded if the local bathymetry appears to be lower than the "nearest" wet node. Basically, you need to clean the flooded area in a last step.

But this is the best method I found so far. It can however be improved if it is performed only on "dry" nodes located at a certain maximal distance from the "wet" nodes (buffer around the flooded mesh). Did not try it yet though.

Best regards
PL

Dear Pl and Christophe,

Never thought of building a second mesh...that could be a viable approach indeed.

On the other hand, could you please post the script based in the arcgis toolbox? (maybe I could get a license here to test it.. ).
Ithink I could mimic the procedure using Gdal or Np/scipy to share it or even improve it...

Regards,

José D.

Hi,

Actually, in my method I don't build a second mesh - I just deal with nodes shapefiles given as input and containing max water depth and max water levels by modifying the field "max water level" on dry nodes with the water level from the nearest wet node. So you need to create this shapefile in a first step in Blue Kenue or fudaa (the script does not use the telemac scripts for reading directly selafin files).

My script was written for ArcGIS 9.3 and does not work with version 10.x because ArcGIS changed their python environment when launching version 10. I still have some fixes to include before I can share it.

Christophe, maybe can this be included in your QGIS based pluggin?

Best regards
PL

Hi
Of course we could include such script in the plugin in the future.
This kind of post treatment is also on our TODO list but needs time to be developed in QGIS...
If some scripts are available, we will be happy to integrate them on our plugin

Regards