Geospatial simulations of dynamic processes: water flow
- to update the 2013 DEM to 2015 state of terrain, perform smooth fusion of lidar and UAS data
- run simulation of surface water flow using path sampling technique on the following DSMs: lidar DEM, UAS DSM, fused lidar+UAS
- optional: assess impact of changes in topography and landuse on water flow pattern and dynamics
- visualize your results using animations
- Use the data from previous assignment, see the Assignments web page
- sparse point clouds for the targeted area from 5 flights in 2015: March (03), June (06), September (09), October 6th (10) and October 29th (11)
- orthophotos generated from the flights Mid Pines
Smooth fusionOne possible approach:
- Use binning to create a 0.2m resolution DEM from lidar for entire area and a DSM from UAV for the region area to be updated. (r.in.lidar), the rasters will have holes but that is OK
- use map algebra to replace all cells in UAV DSM that fall more than 0.5m below lidar DEM with nulls - these are errors that we have observed in the area with trees with points over 10m below the surface
- use map algebra to replace the lidar data in the fields (use the fields mask) by UAV data in the binned DEM: you can patch them exactly at the boundary or leave a narrow band of no-data
- create smooth DEM based on the combined raster with holes to points (select one of the options
- reinterpolate the raster with holes with r.fill.nulls
- export the raster as points or convert it to points and interpolate the smooth raster using v.surf.rst - r.out.xyz may be the easiest, r.to.vect may work too
- randomly sample the combined DEM with output as 3D vector points and interpolate the smooth raster using v.surf.rst
Surface water flow modeling
- Run r.watershed to identify watershed boundaries to make sure you run the simulations in area with at least one complete watershed.
- Select subregion with a watershed and run path sampling based simulation of water flow (r.sim.water) using the lidar DEM, UAV DSM and the fused DEM (lidar updated with UAS data).
- Create animation for the results, discuss differences.
- Select additional UAS DSM from our time series and run water flow simulation, evaluate for change in flow patterns.
Workflow in GRASS GIS 7
Smooth fusionImport UAV data as points, set region and interpolate:
Note: Importing the point clouds is time consuming, but using flags -t (do not create attribute table), -b (do not build topology) and -r (limit import to the current region) make the procedure much faster.
We will create a field mask based on the difference between the UAV surface and lidar DEM. Link to lidar DEM is provided on Assignments page. We will assume that if the difference is larger than 30 cm, then it represents vegetation. Set appropriate color table for differences we used in one of the previous assignments. We will use r.grow module to dilate the mask.
g.region n=219657 s=219320 e=637082 w=636731 res=0.3 -p v.in.lidar -trbo input=2015_06_20_Points_11GCP_Ncspm.las output=uav_06_all v.surf.rst npmin=100 input=uav_06_all elevation=uav_06_interp tension=30 smooth=1
The mask now represents areas which will be replaced by lidar data. We will fuse interpolated UAV and lidar raster data using raster algebra weighted average expression. We compute weights as distance from masked cells. We evaluate the results with shaded relief.
r.mapcalc "diff_uav_06_lidar = uav_06_interp - mid_pines_lidar2013_dem" r.colors map=diff_uav_06_lidar rules=dif_lidar_uav.txt r.mapcalc "diff_uav_06_lidar_mask = if (uav_06_interp - mid_pines_lidar2013_dem > 0.3, 1, null())" r.grow input=diff_uav_06_lidar_mask output=diff_uav_06_lidar_mask_grow radius=3.01 new=1
Now we will merge our result into larger area covered by lidar using the same technique:
r.grow.distance -m input=diff_uav_06_lidar_mask_grow distance=distance_mask r.mapcalc merged_lidar_uav = if(distance_mask > 10, uav_06_interp, (1 - distance_mask/10) * mid_pines_lidar2013_dem + (distance_mask/10) * uav_06_interp)" r.relief input=merged_lidar_uav output=merged_lidar_uav_relief zscale=100
g.region n=n-10 s=s+10 e=e-10 w=w+10 r.mapcalc "rectangle = 1" g.region raster=mid_pines_lidar2013_dem r.grow.distance -m input=rectangle distance=distance_rectangle r.mapcalc "merged_lidar_uav_full = if(distance_rectangle > 10, mid_pines_lidar2013_dem, (1 - distance_rectangle/10) * merged_lidar_uav + (distance_rectangle/10) * mid_pines_lidar2013_dem)" r.relief input=merged_lidar_uav_full output=merged_lidar_uav_full_relief zscale=100
Water flow simulationThe aim of the processing is to investigate the flow pattern simulated on different DSMs. The lidar 2013 and the timeseries of UAS-based DSMs will be used as input for simulations using r.sim.water. First, we examine the watershed boundaries to see which area to include:
r.watershed elevation=merged_lidar_uav_full threshold=50000 accumulation=merged_flowacc basin=merged_basin r.to.vect -t input=merged_basin output=merged_basin type=area
We will compute water flow on our fused DEM within our selected region. We need first order spatial derivatives which can be derived either using r.slope.aspect or during interpolation phase with v.surf.rst.
The simulation writes water depth and discharge raster maps so we can look at first results soon although the simulation takes a lot of time. Depending on your processing power, you can make the region smaller.
g.region n=219727 s=219319 e=637191 w=636701 res=0.3 -p r.slope.aspect elevation=merged_lidar_uav_full dx=merged_dx dy=merged_dy
r.sim.water -t elevation=merged_lidar_uav_full dx=merged_dx dy=merged_dy rain_value=30 man_value=0.15 depth=merged_depth discharge=merged_disch nwalkers=100000 niterations=20 output_step=1 hmax=0.2 halpha=8.0 hbeta=1.0
Animation and visualisation the resultsAnimation tool allows to create animated gifs from a series of raster maps. Load all the *_depth maps to the animation tool. Export the animation as GIF file. You can use ImageMagick to add a pause in the end of the simulation:
convert animation.gif \( +clone -set delay 300 \) animation_with_pause.gif
The water flow pattern is better visible in 2d when it's draped over a shaded relief.
You can also use shaded maps in the animation tool.
d.shade shade=merged_lidar_uav_full_relief color=merged_depth.20 brighten=30
Validation using orthophotos (optional)In order to determine which DSM (lidar based or UAS based) gives more realistic simulation results, the depth maps can be compared with orthophotos. The flights in March, June and both of the October flights were executed after intense rainfall and some of the remaining water (puddles on the road) is still visible.
The orthophotos are available in the .pack format. Use r.unpack for import.
Display the depth simulation rasters after 20 min (the last in each time series) for each simulation. Change the transparency to see the orthophoto below.