NCSU GIS 714:
Geospatial Computing and Simulations

Data simulation

Resources:

For animating in GRASS GIS 7, see the video instructions or instructions for the Spatio-temporal data handling and visualization in GRASS GIS workshop

Start GRASS GIS

Start GRASS - click on GRASS icon or type
grass

In startup pannel set GRASS GIS Database Directory to path to datasets, for example on MS Windows, C:\Users\myname\grassdata.
For GRASS Location select nc_spm_08_grass7 (North Carolina, State Plane, meters) and for GRASS Mapset create a new mapset (called e.g. HW2_data_simulation).
Click Start GRASS session.

Change working directory to where you will keep your external files (such as color ramps, time series layers, etc.):
Settings > GRASS working environment > Change working directory > select/create any directory
or type cd (stands for change directory) into the GUI Console and hit Enter: 

cd

Now you can use the commands from the assignment requiring external text files without the need to specify the full path to the file.

Simple deterministic surfaces

Create surfaces defined by mathematical functions. First, set region to a small rural watershed area at 1m resolution using a predefined region: 

g.region rural_1m -p

Using map algebra generate a tilted plane with 20% slope, slightly rotated, with 50m offset. Cut off the plane at the crossection with elev_lid792_1m surface, but let it protrude 2m to make it visible above the surface (see the equation in the lecture): 

r.mapcalc "tiltplane = 0.02*row()+0.2*col()+50"
r.mapcalc "tiltpl_under = if(tiltplane < elev_lid792_1m + 2,tiltplane,null())"

Generate a plane that crosses the surface along the region diagonal with 10% tilt and protrudes 5m through the elev_lid792_1m surface. Keep in mind that you will need to add an offset so that the plane intersects the elevation surface within the given region - experiment with some offesets and visualize the plane with elev_lid792_1m in 3D. Provide the equation in mathematical and map algebra form.

In the lecture we have shown mathematical surface with hills and valleys generated using the following function 

r.mapcalc "waves = sin(0.4*row())+0.3*cos(col())+80"

Using map algebra generate surface with hills, depressions or valleys that is different from the one shown in the lecture, but is close enough to the elev_lid792_1m surface that you can visualize them together in 3D.

Stochastic methods

Investiagte impact of noise with different distributions on water flow pattern, watershed and stream delineation using random surfaces. Read the manual pages for the relevant commands to understand the methods and parameters.

First we explore different types of random surfaces: uniform with values within the interval [-3,3], Gaussian distribution around the mean=0 with sigma=0.3, and spatially dependent gaussian distribution with distance 10m. Then we add these noise surfaces to our elevation surface: 

g.region rural_1m -p
r.surf.random out=random_neg3_3 min=-3 max=3
r.colors random_neg3_3 co=bcyr
d.legend -v -s -d random_neg3_3 at=5,50,7,10
r.univar random_neg3_3
r.mapcalc "elev_lidnoise = elev_lid792_1m + random_neg3_3/10"
r.colors elev_lidnoise co=elevation

r.surf.gauss out=gauss_m0_s03 mean=0 sigma=0.3
r.colors gauss_m0_s03 co=differences
r.univar gauss_m0_s03
r.mapcalc "elev_lidnoise_gauss = elev_lid792_1m + gauss_m0_s03"
r.colors elev_lidnoise_gauss co=elevation

r.random.surface output=gauss_spatial_d10 distance=10
r.univar gauss_spatial_d10
r.mapcalc "gauss_surfspat = (gauss_spatial_d10 - 134.441) * (0.3/254)"
r.mapcalc "elev_lidnoise_gspat = elev_lid792_1m + gauss_surfspat"
r.colors elev_lidnoise_gspat co=elevation
Can you explain the constants used in the mapcalc command for gauss_surfspat? Where do these constants come from and why are they needed?

Visualy compare flow accumulation petterns and watersheds derived from elev_lid792_1m, elev_lidnoise, elev_lidnoise_gauss, elev_lidnois_gspat: 

r.watershed -a elev_lid792_1m thresh=5000 accum=accum_5K drain=draindir_5K basin=basin_5K
r.watershed -a elev_lidnoise thresh=15000 accum=accum_15K_noise drain=draindir_15K_noise basin=basin_15K_noise
r.watershed -a elev_lidnoise_gspat thresh=10000 accum=accum_10K_noisegspat drain=draindir_10K_noisegspat basin=basin_10K_noisegspat
r.to.vect basin_5K out=basin_5K type=area
r.to.vect basin_10K_noisegspat out=basin_10K_noisegspat type=area
r.to.vect basin_15K_noise out=basin_15K_noise type=area
What can you say about the spatial pattern of flow accumulation derived from these three DEMs?

Optional: For quantitative comparison, you can extract selected basins and compare their area using r.report, compute distances between the boundaries, extract streams using r.stream.extract and compare their characteristics using the r.stream* add-ons.

Random fractal surfaces

First we generate two series of surfaces with dimension 2.9 to 2.01. The series shows creation of a fractal surface. Register the series in temporal framework, assign the rasters a common color table and create an animation (feel free to select larger n or a different dimension). Prepare text file fractal_dem290_series.txt with the names and "years" in your series 
fractal_dem_d290.1|1001
fractal_dem_d290.2|1002
fractal_dem_d290.3|1003
fractal_dem_d290.4|1004
fractal_dem_d290.5|1005
fractal_dem_d290.6|1006
fractal_dem_d290.7|1007
fractal_dem_d290.8|1008
fractal_dem_d290.9|1009
fractal_dem_d290.10|1010
and a similar one for the fractal_dem201_series.txt. 
fractal_dem_d201.1|1001
fractal_dem_d201.2|1002
fractal_dem_d201.3|1003
fractal_dem_d201.4|1004
fractal_dem_d201.5|1005
fractal_dem_d201.6|1006
fractal_dem_d201.7|1007
fractal_dem_d201.8|1008
fractal_dem_d201.9|1009
fractal_dem_d201.10|1010
Then generate the series, register them and visualize, change color table, or run simple analytics. Run watershed analysis on the final surfaces. 
g.region rural_1m -p
r.surf.fractal out=fractal_dem_d290 dim=2.90 n=10
r.info fractal_dem_d290
t.create output=fractal_dem_290 type=strds temporaltype=relative title="Fractal DEMs with d=2.90" description="Generated data"
t.register input=fractal_dem_290 type=raster file=fractal_dem290_series.txt unit=year
t.rast.list in=fractal_dem_290
t.rast.colors in=fractal_dem_290 co=bcyr
t.rast.univar in=fractal_dem_290

r.surf.fractal out=fractal_dem_d201 dim=2.01 n=10
r.info fractal_dem_d201
t.create output=fractal_dem_201 type=strds temporaltype=relative title="Fractal DEMs with d=2.01" description="Generated data"
t.register input=fractal_dem_201 type=raster file=fractal_dem201_series.txt unit=year
t.rast.colors in=fractal_dem_201 co=bcyr

r.watershed -a fractal_dem_d290 thresh=15000 accum=f290_accum_15K drain=f290_draindir_15K basin=f290_basin_15K
r.watershed -a fractal_dem_d201 thresh=15000 accum=f201_accum_15K drain=f201_draindir_15K basin=f201_basin_15K

Rescale the fractal surface to [-0.3,0.3] and add to our DEM as a simulated noise. The water flow pattern may represent flow over rough surface, where the roughness is represented in the DEM (e.g. as a grass cover) rather than a constant value (e.g. Mannings coefficient). 

r.univar fractal_dem_d290
r.mapcalc "elev_noisefractal = elev_lid792_1m +0.6*(fractal_dem_d290-50.9649)/3754.7"
r.colors elev_noisefractal co=elevation
r.univar elev_noisefractal
r.watershed -a elev_noisefractal thresh=15000 accum=_accum_15K_noisefrac drain=draindir_15K_noisefrac basin=basin_15K_noisefrac
Can you explain the constants used in the mapcalc command for elev_noisefractal? Where do these constants come from and why are they needed?

Compute elevation surface from the lidar bare ground data with different tension parameter and compare the impact on flow modeling. Also compare the flow pattern from simulated noise with the noise from lidar data preserved in the surface interpolated with higher tension parameter: 

g.region raster=elev_lid792_1m -p
v.surf.rst elev_lid792_bepts elev=elev_lid_default npmin=120 
v.surf.rst elev_lid792_bepts elev=elev_lid_t120 npmin=120 ten=120
r.watershed -a elevation=elev_lid_default threshold=5000 accumulation=accum_5K_liddef drainage=draindir_5K_liddef basin=basin_5K_liddef
r.watershed -a elevation=elev_lid_t120 threshold=5000 accumulation=accum_5K_lidt120 drainage="draindir_5K_lidt120 basin=basin_5K_lidt120
Which of the surfaces with simulated noise produce the flow pattern closest to the surface with noise from lidar point cloud?

Random point data

In this example we analyze how errors in point position can influence viewshed analysis. First, install r.viewshed.cva addon: 

g.extension r.viewshed.cva
Generate random points and perturb them. For both points set compute a cumulative viewshed and compare them: 
g.region raster=elev_lid792_1m
v.random output=random_view npoints=10 seed=4
r.viewshed.cva input=elev_lid792_1m vector=random_view output=cumulative
v.perturb input=random_view output=random_view_perturb parameters=20 seed=1
r.viewshed.cva input=elev_lid792_1m vector=random_view_perturb output=cumulative_perturb
r.report map=cumulative units=p
r.report map=cumulative_perturb units=p
What is the impact of error in point position on total viewshed area? What is result for the perturbation with normal distribution? (read the manual page to leaern how to run it).

Read Starek, M.J., et al., 2020, Viewshed simulation and optimization for digital terrain modelling with terrestrial laser scanning, International Journal of Re-mote Sensing, 41:16, 6409-6426 to learn how simulated annealing can be used to find optimal locations for a terrestrial laser scanner to cover the most area.

Additional commands you can use to simulate point data in raster or vector representation: r.random, r.random.cells, r.sample.category, v.perturb, v.random, v.kcv

Generating patches

Create several landscape cover scenarios, for example, to explore the following scenario - if we need to convert 30% of forest to development, is it better to create a single compact area or many smaller areas? Read the manual pages for r.pi to understand what the commands perform and feel free to change the examples below to get more realistics results. 
g.extension r.pi
g.region raster=landclass96 -p
r.pi.nlm.circ output=circ1 landcover=30 count=20 size=526,474
g.region raster=elevation -p
r.pi.nlm.circ output=circ2 landcover=30 count=500 size=1500,1350
r.pi.nlm input=landclass96 output=nlm.3 nullval=1,5,6,7 landcover=10 sharp=1.0
r.report circ1 unit=p
r.report circ2 unit=p
r.report nlm.3 unit=p