Assignment 5B

Analysis of multitemporal UAS data and its applications

Authors

Anna Petrasova

Corey White

Outline

Introduction to GRASS temporal framework
Analyze the time series of UAS DEMs
Estimate crop biomas and its temporal change
Evaluate impact of vegetation change on viewshed extent
Estimate volume of eroded and deposited soil

Data

You should have everything already for Mid Pines area, get anything missing from the Course logistics web page:

GRASS Project: LakeWheeler_NCspm
2013 lidar interpolated DEM and DSM: download, use r.unpack.
Fields polygon as packed raster use r.unpack in your mapset to import it
Text file with timeseries

Tools

GRASS GIS 8.4
GRASS Addon: r.patch.smooth

Install r.patch.smooth with the following command:

g.extension r.patch.smooth

Analyze time series of UAS DSMs

Create a new mapset and copy the UAS DSMs

Start GRASS GIS in Lake_Wheeler_NCspm Location and create a new mapset “assignment5b”. Open GRASS GIS in this mapset

g.copy raster=2015_03_18_DSM_agi_6GCP@PERMANENT,2015_03_18_DSM_agi_6GCP
g.copy raster=2015_06_20_DSM_agi_11GCP@PERMANENT,2015_06_20_DSM_agi_11GCP
g.copy raster=2015_07_27_DSM_agi_2GCP_4ID@PERMANENT,2015_07_27_DSM_agi_2GCP_4ID 
g.copy raster=2015_08_14_DSM_agi_3GCP_4ID@PERMANENT,2015_08_14_DSM_agi_3GCP_4ID 
g.copy raster=2015_09_04_DSM_agi_6GCP@PERMANENT,2015_09_04_DSM_agi_6GCP
g.copy raster=2015_09_21_DSM_agi_6GCP@PERMANENT,2015_09_21_DSM_agi_6GCP
g.copy raster=2015_10_06_DSM_agi_8GCPs@PERMANENT,2015_10_06_DSM_agi_8GCPs

Register the time series in temporal framework

First, create a space-time raster dataset (STRDS) with absolute temporal type:

t.create output=uas_dsm type=strds temporaltype=absolute semantictype=mean title="UAS_DSM" description="UAS-based 2015 DSMs processed by Agisoft"

Now, create a text file named “series.txt” with the following content in your working directory, or download the prepared file into your working directory.

These are the DSMs with their time stamps:

2015_03_18_DSM_agi_6GCP|2015-03-18
2015_06_20_DSM_agi_11GCP|2015-06-20
2015_07_27_DSM_agi_2GCP_4ID|2015-07-27
2015_08_14_DSM_agi_3GCP_4ID|2015-08-14
2015_09_04_DSM_agi_6GCP|2015-09-04
2015_09_21_DSM_agi_6GCP|2015-09-21
2015_10_06_DSM_agi_8GCPs|2015-10-06

t.register input=uas_dsm file=series.txt

Now check it was registered correctly:

t.info input=uas_dsm

See different ways we can list raster maps registered in the dataset:

t.rast.list input=uas_dsm columns=name,start_time where="start_time >= '2015-08-01'"
t.rast.list input=uas_dsm columns=name,start_time,min,max

Verify the listed raster maps are actually in the database:

g.list type=raster pattern="2015*"

You can also use the Data tab in Layer Manager.

We can set common color table to all maps in the dataset:

t.rast.colors input=uas_dsm color=elevation
d.rast 2015_06_20_DSM_agi_11GCP@assignment5b
d.rast 2015_10_06_DSM_agi_8GCPs@assignment5b

To get better idea about the dates and spatial extents, we can use Timeline tool accessible from Temporal - GUI tools - Timeline tool.

Select dataset uas_dsm, you can click on the drawn data points, and then also check 3D plot of spatio-temporal extents. Stretch the window as needed.

We can plot the elevation timeseries using Temporal Plot Tool.

Open the tool from Temporal - GUI tools - Temporal plot tool.
Select dataset uas_dsm and click on the arrow next to the X and Y coordinates field.
Then click somewhere on the Map Display and press Draw.

You should see a plot of elevation values. Try different places (road, field, building).

Find the area mapped by all surveys

We will find area mapped by all surveys. First we will set the extent to the maximum bounding box, use t.info to get the values:

t.info input=uas_dsm@assignment5b -g
g.region n=219942 s=219037 e=637439 w=636456 -pa res=0.3

Next, we use t.rast.series to derive raster representing the number of overlapping DSMs. By using -n flag we create masking raster representing the intersection of all DSMs, and then we change computational region to the extent of its non-null bounding box:

t.rast.series input=uas_dsm@assignment5b method=count output=count
t.rast.series input=uas_dsm@assignment5b method=count output=intersection -n
d.rast count
d.rast intersection
g.region zoom=intersection -p
r.mask raster=intersection

Note

How large is the area? (use r.report)

Temporal aggregation

Our space-time raster dataset contains two DSMs from September:

t.rast.list input=uas_dsm@assignment5b columns=name,start_time

We will temporally aggregate the dataset by month by averaging all DSMs within each month (in this case averaging the September DSMs):

t.rast.aggregate input=uas_dsm@assignment5b output=uas_dsm_aggr basename=uas_dsm_aggr suffix=time granularity="1 months" method=average
t.rast.list input=uas_dsm_aggr columns=name,start_time

Minimum elevation (core)

You should already have lidar DSM and DEMs if you haven’t done already, unzip it in your working directory and unpack:

r.unpack -o input=mid_pines_lidar2013_dsm.pack
r.unpack -o input=mid_pines_lidar2013_dem.pack

Note

Find the minimum elevation (core) - how close it is to lidar bare ground?

t.rast.series input=uas_dsm method=minimum output=minimum
r.mapcalc "diff_lidar_uas = mid_pines_lidar2013_dem - minimum"
r.colors map=diff_lidar_uas color=differences
d.rast diff_lidar_uas

Parts of the minimum raster are significantly below lidar ground, find out which DSM is causing it:

t.rast.series input=uas_dsm@assignment5b method=min_raster output=min_raster
d.rast min_raster

Display resulting raster min_raster and query it (select the layer in Layer Manager and use Query tool in Map Display) or use r.what:

r.what map=min_raster coordinates=636879,219432
r.what map=min_raster coordinates=637232,219651

Value 2 represents third (numbering starts with 0) raster in our time series. Since there is some problem with the third and fourth raster (2015_07_27_DSM_agi_2GCP_4ID, 2015_08_14_DSM_agi_3GCP_4ID@assignment5b), we will exclude them from the analysis:

t.rast.series input=uas_dsm@assignment5b method=minimum output=minimum where="start_time < '2015-07-01' or start_time > '2015-09-01'" --o
r.mapcalc "diff_lidar_uas = mid_pines_lidar2013_dem - minimum" --o
r.colors map=diff_lidar_uas color=differences

Apply now the color table dif_lidar_uav.txt we used for highlighting differences in the last assignment:

r.colors map=diff_lidar_uas rules=dif_lidar_uav.txt

Explain what you see.

Let’s look now which temporal datasets we have:

t.list columns=id,number_of_maps,start_time

Estimate crop biomass

We will use raster algebra on the time series to compute crop biomass. First we need to make sure that we take into account only fields (not forests).

You can use get the packed raster representation of fields (link is at the top).

We assume crop is anything higher than 0.3 and lower than 2 m (to exclude the building for example):

r.unpack -o fields.pack
t.rast.mapcalc inputs=uas_dsm@assignment5b expression="if (fields, if(uas_dsm - mid_pines_lidar2013_dem > 0.3 && uas_dsm - mid_pines_lidar2013_dem < 2, uas_dsm - mid_pines_lidar2013_dem, null()))" output=veg_uas_lidar basename=veg_uas_lidar
t.rast.univar input=veg_uas_lidar where="start_time < '2015-07-01' or start_time > '2015-09-01'"
d.rast veg_uas_lidar_2
r.report -n map=veg_uas_lidar_2 units=me nsteps=1
r.univar map=veg_uas_lidar_2

Use Python shell (tab in Layer Manager) to compute the biomass volume for June (paste each line into Python shell and press Enter):

sum = gs.parse_command('r.univar', map='veg_uas_lidar_2', flags='g')['sum']
float(sum) * gs.region()['nsres'] * gs.region()['ewres']

Remove the mask before we proceed in further analysis:

r.mask -r

Evaluate impact of vegetation change on viewshed extent

We would like to set up a webcam to monitor this area. First we will compute a cumulative viewshed of this area to find places from which we will have likely good visibility. We pick several places distributed all over the area from which we compute viewsheds on a lidar DSM with observer height of 3 meters (webcam height).

g.region raster=mid_pines_lidar2013_dsm res=1 -pa
r.viewshed input=mid_pines_lidar2013_dsm output=vshed_1 coordinates=637087,219376 observer_elevation=3
r.viewshed input=mid_pines_lidar2013_dsm output=vshed_2 coordinates=637004,219653 observer_elevation=3
r.viewshed input=mid_pines_lidar2013_dsm output=vshed_3 coordinates=636782,219374 observer_elevation=3
r.viewshed input=mid_pines_lidar2013_dsm output=vshed_4 coordinates=636839,219150 observer_elevation=3
r.viewshed input=mid_pines_lidar2013_dsm output=vshed_5 coordinates=636665,219527 observer_elevation=3
r.viewshed input=mid_pines_lidar2013_dsm output=vshed_6 coordinates=636863,219503 observer_elevation=3

Now we compute the cumulative viewshed as the number of places (from the coordinates above) from which a cell is visible. This gives us better idea which places are visually more prominent and therefore likely to have good view of the fields.

r.series input=vshed_1,vshed_2,vshed_3,vshed_4,vshed_5,vshed_6 output=cum_viewshed method=count
d.rast cum_viewshed

Display the result with legend to see which places are potentially good for placing a webcam.

Since we need the webcam to capture the field during entire year we can run t.rast.series with method=maximum on the time series of UAS DSMs to derive the DSM series envelope surface and use it to analyze the viewshed. Of course, this assumes that the crops will be the same the following year. We exclude the problematic DSMs as we did while deriving the core surface:

g.region raster=mid_pines_lidar2013_dsm
t.rast.series input=uas_dsm@assignment5b method=maximum output=maximum where="start_time < '2015-07-01' or start_time > '2015-09-01'"

We need to update the lidar-based DSM with the UAS envelope data and analyze the viewshed while taking into account the corn growing in the field. First we need to make sure that we take into account entire field to create a lidar-based DSM updated in the fields using UAS data.

r.mapcalc "maximum_clip = if (fields, maximum)"
r.patch.smooth input_a=maximum_clip input_b=mid_pines_lidar2013_dsm output=lidar_maximum_dsm smooth_dist=10
r.relief lidar_maximum_dsm out=lidar_maximum_dsm_relief zscale=5
d.rast lidar_maximum_dsm_relief

Look at the relief map around the building - where do people park their cars?

Now we selected a location based on the cumulative viewshed and look at the difference between the viewshed based on lidar and the envelope surface:

r.viewshed input=mid_pines_lidar2013_dsm output=viewshed_lidar_dsm observer_elevation=3 coordinates=636917,219223
r.viewshed input=lidar_maximum_dsm output=viewshed_maximum observer_elevation=3 coordinates=636917,219223

r.colors viewshed_lidar_dsm co=reds
r.colors viewshed_maximum co=greens

To compare the viewsheds display the resulting maps with transparency.

Note

Discuss the result. Can you find a better place for the webcam? Would placing the webcam higher help?

Estimate volume of a building

We need to first extract the mask of the building:

g.region n=219388.5 s=219366.9 w=637100.1 e=637138.5 res=0.3 -pa
r.mapcalc "building = 2015_06_20_DSM_agi_11GCP - mid_pines_lidar2013_dem"
r.mapcalc "building_mask = if(building > 2, 1, null())"

Now we compute the volume

r.volume input=building clump=building_mask

Note

What is the volume of the building?