Lab: Analysis of UAS-derived DSMs and orthophotos

Analysis of UAS data processing results

Open in Colab

GIS based analysis of UAS data processing outputs

The steps below are for the GRASS GUI or command line. The notebook (badge above) runs the same analysis in Python in Google Colab, with the maps and statistics produced for you.

In this lab we analyze DSMs derived from UAS imagery in terms of their accuracy and the type of distortions they carry, depending on whether GCPs were used and on the processing software. Then we map terrain change between two UAS surveys (vegetation growth, erosion, and other impacts) and separate real change from errors and artifacts.

  • Including GCPs in the processing improves the results. Does it have an effect on the phenomenon you are investigating?
  • Agisoft Metashape is only one of the software solutions available for UAS data processing. The Lake Wheeler GRASS project holds DSMs of the same flight generated in Agisoft, Pix4D, and Trimble Business Center. Compare the differences between them.
  • In Assignment 2B you used given parameter values at each processing step. You can change modes, quality, values, and options to investigate whether they improve the results (Metashape manual).

Data

Getting started

Start GRASS in the Lake_Wheeler_NCspm project and create a new mapset for this assignment, for example analysis.

Every comparison in this lab produces a difference raster. Look at each one in two ways:

  • the map, with a color table that spreads the colors over the range where most of the values are
  • the statistics, from r.univar (in the GUI: Raster > Reports and statistics > Univariate raster statistics, check Extended statistics):
r.univar -e map=GCP_noGCP

The mean of a difference raster is the vertical shift between the two DSMs. The standard deviation and the range between the center and the edges of the area describe the bowl. Note these numbers for each comparison; the report asks for them.

To enter a color table, open r.colors (in the GUI: right click the layer, Set color table), go to the Define tab, and type the rules into enter values directly.

DSMs comparison - influence of GCPs

sample_DSM was generated with 4 GCPs; sample_DSM_noGCPs is the same imagery processed without them. Add the GCP_12_degrees vector map to the display to see the GCPs; the four inside the extent of sample_DSM (UAV 3, 4, 6, and 8) are the ones that were used.

Set the computational region to the DSM and compute the difference with raster map algebra:

g.region raster=sample_DSM
r.mapcalc expression="GCP_noGCP = sample_DSM - sample_DSM_noGCPs"
r.univar -e map=GCP_noGCP

The statistics show that the DSM without GCPs is about 36 m too low everywhere, with most of the values between 35 and 39 m. A color table that spreads its colors over that band shows the shape of the error instead of the shift:

0 red
35 orange
36 yellow
37 cyan
38 aqua
39 blue

Trees and buildings are reconstructed differently in the two DSMs and fall outside the band; they show as red or blue specks.

Important

What is the shape of the difference? Consider the distribution of the GCPs and the shape of the area.

Compare DSMs generated in different software

For the flight of 2015-06-20, the PERMANENT mapset holds DSMs with 11 GCPs generated in three software environments (Trimble Business Center, Agisoft, and Pix4D), and DSMs generated without GCPs in Agisoft and Pix4D.

Check the bowl effect for the Agisoft products:

g.region raster=2015_06_20_DSM_agi_11GCP
r.mapcalc expression="agi_GCP_agi_noGCP = 2015_06_20_DSM_agi_11GCP - 2015_06_20_DSM_agi_noGCP"
r.univar -e map=agi_GCP_agi_noGCP

Color table:

0 red
35 orange
36 yellow
37 cyan
38 aqua
42 blue

Check the bowl effect for the Pix4D products:

r.mapcalc expression="p4d_GCP_p4d_noGCP = 2015_06_20_pix4d_11GCP_dsm - 2015_06_20_DSM_pix4d_NoGCP"
r.univar -e map=p4d_GCP_p4d_noGCP

Color table:

50 red
69 orange
70 yellow
71 cyan
72 aqua
75 blue
Important

The two color tables cover different ranges, so do not compare the maps by color alone. Use the statistics. Which software has the larger vertical shift (the mean)? Which one has the more curved bowl (the standard deviation, and the range between the center and the edges)? These are two different questions and they do not have the same answer.

Now compare the three DSMs generated with GCPs against each other. The sign of each difference follows the name: agi_trimble is positive where the Agisoft DSM is higher than the Trimble DSM.

r.mapcalc expression="agi_trimble = 2015_06_20_DSM_agi_11GCP - 2015_06_20_DSM_Trimble_11GCP"
r.mapcalc expression="p4d_agi = 2015_06_20_pix4d_11GCP_dsm - 2015_06_20_DSM_agi_11GCP"
r.mapcalc expression="p4d_trimble = 2015_06_20_pix4d_11GCP_dsm - 2015_06_20_DSM_Trimble_11GCP"

Apply a color table from -1 m to 1 m to each of them, with white for differences within 10 cm:

-40 red
-1 orange
-0.5 yellow
-0.1 grey
0 white
0.1 grey
0.5 cyan
1 aqua
35 blue

The Trimble DSM covers a smaller area than the other two, so the maps that include it are empty (null) outside its extent. That is expected.

Important

Compare the maps. How do they relate to each other? What patterns (artifacts) can you recognize? Where are the differences: along tree lines and buildings, at the edges of the block, in stripes along the flight lines?

Terrain analysis of the DSMs

Terrain analysis shows the bowl more clearly than the elevation does: a dome becomes a slow gradient in slope and a radial pattern in aspect, on top of the real terraces and tree lines.

Run r.slope.aspect and r.relief on the sample DSM with and without GCPs and compare the results side by side:

g.region raster=sample_DSM
r.slope.aspect elevation=sample_DSM slope=sample_DSM_slope aspect=sample_DSM_aspect pcurvature=sample_DSM_curv
r.relief input=sample_DSM output=sample_DSM_relief
r.slope.aspect elevation=sample_DSM_noGCPs slope=sample_DSM_noGCP_slope aspect=sample_DSM_noGCP_aspect pcurvature=sample_DSM_noGCP_curv
r.relief input=sample_DSM_noGCPs output=sample_DSM_noGCP_relief

You can also install GRASS addons with g.extension:

and apply them to the sample data or to your own data. The region is still set to sample_DSM.

Install and run r.local.relief:

g.extension extension=r.local.relief
r.local.relief input=sample_DSM output=local_relief_sample_DSM
r.local.relief input=sample_DSM_noGCPs output=local_relief_sample_DSM_noGCPs

Install and run r.shaded.pca:

g.extension extension=r.shaded.pca
r.shaded.pca input=sample_DSM output=shaded_pca_sample_DSM
r.shaded.pca input=sample_DSM_noGCPs output=shaded_pca_sample_DSM_noGCPs

Install and run r.skyview:

g.extension extension=r.skyview
r.skyview input=sample_DSM output=skyview_sample_DSM
r.skyview input=sample_DSM_noGCPs output=skyview_sample_DSM_noGCPs

Detect terrain change using two UAS surveys

The PERMANENT mapset also holds the DSM of the flight of 2015-10-06, processed in Agisoft with 8 GCPs. Subtract it from the June DSM:

g.region raster=2015_06_20_DSM_agi_11GCP
r.mapcalc expression="diff_jun_oct_agis = 2015_06_20_DSM_agi_11GCP - 2015_10_06_DSM_agi_8GCPs"
r.univar -e map=diff_jun_oct_agis

The difference is June minus October. Positive values mean the surface was higher in June: a loss of elevation by October, for example a harvested crop. Negative values mean the surface was higher in October: a gain, for example a crop that grew over the summer. Keep this sign in mind when you read the map.

Apply the same -1 m to 1 m color table as above:

-40 red
-1 orange
-0.5 yellow
-0.1 grey
0 white
0.1 grey
0.5 cyan
1 aqua
35 blue

The roads did not change between June and October, so the difference on them is the error of the comparison. Note how large it is; anything smaller than that elsewhere cannot be called change.

Classify and measure the change

Call differences beyond half a meter loss or gain and everything in between stable, then let r.report sum the area of each class in hectares. Half a meter is well above the error on the roads; change the threshold if your reading of the roads suggests a different value, and say so in the report.

r.mapcalc expression="change_class = if(diff_jun_oct_agis > 0.5, 1, if(diff_jun_oct_agis < -0.5, 3, 2))"
r.category map=change_class separator=: rules=- <<EOF
1:loss (June higher)
2:stable
3:gain (October higher)
EOF
r.colors map=change_class rules=- <<EOF
1 orange
2 white
3 blue
EOF
r.report -n map=change_class units=h,p

In the GUI, enter the category labels and the color rules through the Define tabs of r.category and r.colors, and run r.report from Raster > Reports and statistics > Report basic information.

Orthophoto for context

To interpret the change you need to know what was on the fields. The November 2015 orthophoto was flown after the harvest; it does not match either DSM date, but it shows the field boundaries, the roads, and the tree lines. Change the working directory to the folder you downloaded it to (Settings > GRASS working environment > Change working directory) and restore it with r.unpack:

r.unpack ortho_2015_11_24.pack -o

Display the orthophoto with the change map on top of it. The change map is opaque, so set its opacity in the Layer Manager (right click the layer, Change opacity level) to about 50 percent, or display only the loss and gain cells:

r.mapcalc expression="change_only = if(change_class == 2, null(), change_class)"
r.colors map=change_only raster=change_class
Important

Which of the loss and gain patches follow field boundaries (crops), which follow tree lines (canopy growth or leaf-off differences), and which follow flight lines or the edges of the block (processing artifacts)? What is the difference on the roads, and what does it say about the accuracy of a change map from two UAS surveys?