Assignment 5B

Analysis of multitemporal UAS data

Published

October 28, 2026

Open in Colab

Task

In this assignment you analyze the seven UAS DSMs of Mid Pines flown between March and October 2015 as a time series. You register them as a space-time raster dataset in the GRASS temporal framework, find the area covered by every survey, and aggregate the series by month. You derive the core surface, compare it with the lidar bare ground to find the two distorted surveys, and recompute the core without them. You then apply the series: crop height and biomass from map algebra on the series, a webcam viewshed on the lidar DSM and on the envelope DSM, and the volume of a building.

Keep in mind that the seven flights were processed with different numbers of GCPs, from two to eleven. The surveys with few GCPs are the ones the core surface exposes. Everything in Parts 2 to 4 is computed inside the intersection of the surveys and under a mask, so the statistics refer to that area only.

Work through the steps below in the GRASS GUI or command line, or open the notebook in Colab (badge above), which runs the same steps in Python and produces the maps and statistics for you. The Timeline Tool and the Temporal Plot Tool in Part 1 are GUI only; the notebook plots the same things. Either way, the report is the same.

Data

Lake Wheeler GRASS project: Lake_Wheeler_NCspm (EPSG:3358; you have it from the earlier assignments; the seven 2015 DSMs are in PERMANENT)
2013 lidar DEM and DSM of Mid Pines (zip with GRASS raster packs; you have it from Assignments 4B and 5A, restore with r.unpack)
Fields raster: fields.pack (GRASS raster pack of the crop fields, 1 inside the fields, null elsewhere)
Series file: series.txt (the DSM names and their dates, for t.register)
Color rules for difference maps: dif_lidar_uav.txt (the same file as in Assignment 5A)

Software

GRASS 8.5
The addon r.patch.smooth, which you installed in Assignment 5A; if not:

g.extension extension=r.patch.smooth

Part 1: Register the time series

Create a mapset and copy the DSMs

Start GRASS in the Lake_Wheeler_NCspm project, create a new mapset multitemporal, and switch to it. Put the data files in your working directory.

g.mapset -c mapset=multitemporal

The temporal database lives in the mapset, and maps from other mapsets cannot be registered in it, so copy the seven DSMs from PERMANENT:

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

The names record the flight date and the number of GCPs used in processing.

Create and register the dataset

Create a space-time raster dataset (STRDS) with absolute time:

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

The file series.txt lists each DSM with its date, one per line, separated by a pipe:

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
t.info input=uas_dsm

t.rast.list lists the registered maps; the where clause is SQL on the time columns, and columns can include statistics of each map:

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

The maps are still ordinary rasters in the mapset, which g.list type=raster pattern="2015*" or the Data tab in the Layer Manager shows. Registration only records their time in the database.

Look at the series

Set one color table for all maps in the dataset and display two of them:

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

The Timeline Tool (Temporal > GUI tools > Timeline tool) shows when each map was taken. Select the dataset uas_dsm, click on the points, and check 3D plot of spatio-temporal extents to see that the flights covered different areas.

The Temporal Plot Tool (Temporal > GUI tools > Temporal plot tool) plots the elevation of every map at one location. Select the dataset uas_dsm, click the arrow next to the X and Y coordinates field, click a point in the Map Display, and press Draw. Try a point on the road, in a field, and on the building.

Part 2: Intersection and aggregation

Area mapped by all surveys

Per-cell statistics over the series are only meaningful where every survey has data. Set the region to the union of the extents, which t.info -g reports, then count the maps with data in each cell with t.rast.series. With the -n flag the count is null wherever any map is null, which gives the intersection:

t.info input=uas_dsm -g
g.region n=219942 s=219037 e=637439 w=636456 res=0.3 -pa
t.rast.series input=uas_dsm method=count output=count
t.rast.series input=uas_dsm method=count output=intersection -n
d.rast map=count
d.rast map=intersection

Shrink the region to the intersection and set it as the mask for the rest of Parts 2 to 4:

g.region zoom=intersection -p
r.mask raster=intersection
r.report -n map=intersection units=h

18.65 hectares

Every raster command from here on reads and writes only inside the mask, until r.mask -r at the end of Part 4.

Monthly aggregation

The dataset has two surveys in September. t.rast.aggregate combines the maps that fall in each interval of a chosen granularity into one map, here the monthly average, and registers the results as a new dataset with interval time stamps:

t.rast.aggregate input=uas_dsm 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,end_time

Six monthly maps result; only the September one averages two surveys.

Part 3: Core surface and distorted surveys

Minimum elevation

Restore the lidar surfaces from Assignment 5A, then compute the core surface, the minimum elevation of the series at each cell, and subtract it from the lidar bare ground:

r.unpack -o input=mid_pines_lidar2013_dsm.pack
r.unpack -o input=mid_pines_lidar2013_dem.pack
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
r.univar map=diff_lidar_uas
d.rast map=diff_lidar_uas

mean: 0.87, minimum: -4.35, maximum: 15.34

A DSM cannot lie below the ground it was flown over, so a core surface well below the lidar DEM points at a distorted survey. Over much of the area the core is close to the lidar ground, but in places it lies meters below it.

Which survey is responsible

method=min_raster writes, for each cell, the index of the map that supplied the minimum, counted from 0 in the order of the time stamps:

t.rast.series input=uas_dsm method=min_raster output=min_raster
d.rast map=min_raster
r.what map=min_raster coordinates=636879,219432
r.what map=min_raster coordinates=637232,219651
r.report -n map=min_raster units=p nsteps=7

636879|219432||2 637232|219651||3

Index 2 is the third map, the July 27 survey, and index 3 the fourth, August 14. Together they supply the minimum in two thirds of the cells. These two flights were processed with only two and three GCPs.

Recompute the core without them

Exclude the two surveys with a where clause on the time stamps and recompute the core and the difference:

t.rast.series input=uas_dsm method=minimum output=minimum where="start_time < '2015-07-01' or start_time > '2015-09-01'" --overwrite
r.mapcalc "diff_lidar_uas = mid_pines_lidar2013_dem - minimum" --overwrite
r.univar map=diff_lidar_uas
r.colors map=diff_lidar_uas rules=dif_lidar_uav.txt
d.rast map=diff_lidar_uas

mean: 0.08

The mean difference drops from 0.87 m to 0.08 m. The color rules file is the one from Assignment 5A: grey within 10 cm, color beyond. What remains colored, and why?

t.list shows the datasets in the mapset:

t.list columns=id,number_of_maps,start_time

Part 4: Crop height and biomass

Map algebra on a series applies one expression to every map and registers the results as a new dataset. We compute the height of the UAS surface above the lidar bare ground, keep it only inside the crop fields and between 0.3 and 2 m, to drop the soil noise below and the building and trees above, and set the rest to null. The fields raster comes from the pack:

r.unpack -o input=fields.pack
t.rast.mapcalc inputs=uas_dsm 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.list input=veg_uas_lidar columns=name,start_time

In the expression, uas_dsm stands for each registered map in turn. The outputs are named veg_uas_lidar_1 to veg_uas_lidar_7, numbered from 1 in time order, so veg_uas_lidar_2 is June 20. Summarize the crop height per survey, leaving out the two distorted ones:

t.rast.univar input=veg_uas_lidar where="start_time < '2015-07-01' or start_time > '2015-09-01'"

The June survey has by far the most cells with crop and a mean height of 0.8 m. Look at it and compute its crop area and volume:

d.rast map=veg_uas_lidar_2
r.report -n map=veg_uas_lidar_2 units=me
r.univar map=veg_uas_lidar_2

The volume of crop above ground is the sum of the heights times the cell area. In the Python tab of the Layer Manager, run:

import grass.script as gs
stats = gs.parse_command("r.univar", map="veg_uas_lidar_2", flags="g")
region = gs.region()
float(stats["sum"]) * region["nsres"] * region["ewres"]

49,605 cubic meters over 62,232 square meters of crop

Remove the mask before Part 5, which works on the whole study area:

r.mask -r

Part 5: Viewsheds and the envelope DSM

Cumulative viewshed

We want to place a webcam that watches the fields. First compute a cumulative viewshed from six candidate points on the lidar DSM, with the camera 3 m above the ground, at 1 m resolution:

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

Each viewshed holds the vertical angle to the visible cells and null elsewhere. Count, for every cell, from how many of the six points it is visible:

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

Display the result with a legend. Cells visible from many points are the prominent places; a camera there would see much of the area. About half of the area is visible from none of the six points.

Envelope DSM

The camera has to see the fields through the whole season, so the test surface is the envelope of the series, the maximum elevation at each cell, fused into the lidar DSM inside the fields. Compute the envelope without the two distorted surveys at the native resolution, clip it to the fields, and fuse it with r.patch.smooth as in Assignment 5A:

g.region raster=mid_pines_lidar2013_dsm -p
t.rast.series input=uas_dsm method=maximum output=maximum where="start_time < '2015-07-01' or start_time > '2015-09-01'"
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 input=lidar_maximum_dsm output=lidar_maximum_dsm_relief zscale=5
d.rast map=lidar_maximum_dsm_relief

This assumes next year’s crop looks like this year’s. Look at the relief around the building: the envelope keeps everything that was ever there, including the parked cars.

Compare the viewsheds

Pick a camera location from the cumulative viewshed, here a point with a wide view, and compute its viewshed on the lidar DSM and on the envelope DSM at 1 m resolution:

g.region raster=mid_pines_lidar2013_dsm res=1 -pa
r.viewshed input=mid_pines_lidar2013_dsm output=viewshed_lidar_dsm coordinates=636917,219223 observer_elevation=3
r.viewshed input=lidar_maximum_dsm output=viewshed_maximum coordinates=636917,219223 observer_elevation=3
r.colors map=viewshed_lidar_dsm color=reds
r.colors map=viewshed_maximum color=greens
r.report -n map=viewshed_lidar_dsm units=h
r.report -n map=viewshed_maximum units=h

viewshed_lidar_dsm: 11.0 hectares viewshed_maximum: 6.4 hectares

Display both with transparency to see where the crop closes the view. Would a higher camera or another candidate point do better?

Part 6: Volume of a building

The difference between a DSM and the bare ground over the footprint of a feature gives its volume. Set the region to the building, compute the height of the June DSM above the lidar DEM, and take everything more than 2 m high as the footprint:

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())"
r.volume input=building clump=building_mask

Total Volume = 897.44

r.volume sums the height over the clump and multiplies by the cell area: 2,473 cells, 223 square meters, an average height of 4 m, and a volume of about 900 cubic meters.

Report

Prepare a report (3 to 5 pages including figures):

Series and intersection: the registered dataset (t.info output or the Timeline Tool), the count and intersection maps, the intersection area, and a temporal plot at a field point and a road point with a sentence on what each curve shows.
Aggregation: the list of monthly maps and what the September aggregate represents.
Core surface: the lidar minus core map before and after excluding the distorted surveys, with their means, the min_raster map, and how you identified the two surveys and why they are distorted.
Crop height and biomass: the June crop height map, the t.rast.univar table for the valid surveys, the June crop area and volume, and a paragraph on how the crop changes through the season.
Viewsheds: the cumulative viewshed, the two viewsheds at your camera point with their areas, and a discussion of how the crop changes the view and where you would put the camera.
Building volume: the volume, how the footprint was defined, and what limits the accuracy of the estimate.

Submission

Upload to Moodle, Assignment 5B:

Your report as a PDF.
Optional: your completed notebook (.ipynb) if you worked in Colab.

Grading

Series and intersection: the dataset is registered with correct dates, the intersection is derived and used as the mask, and the temporal plots are read correctly.
Aggregation: the monthly dataset is produced and explained.
Core surface: the two distorted surveys are identified from the core and min_raster, and the improvement after excluding them is quantified.
Crop height and biomass: the height series is derived with map algebra, the June area and volume are computed, and the seasonal change is interpreted.
Viewsheds: the cumulative and the two single viewsheds are computed and the effect of the crop on visibility is explained.
Building volume: the volume is computed and its limitations are stated.
Interpretation: conclusions tie back to the lecture: space-time dataset, intersection, aggregation, core and envelope, distorted surveys.
On-time submission.