GIS/MEA 584:
Mapping and Analysis Using UAS

Multispectral imagery processing and analysis

Outline:
  • process multispectral imagery in Agisoft Metashape
  • calculate Vegetation indices in GRASS GIS
This assignment is based on chapter 17 of the book “Fundamentals of capturing and processing drone imagery and data” (available through NCSU library: available through NCSU library)

Data:

Agisoft part: GRASS GIS part:

Software:

In the assignment we will be using GRASS GIS 7.8 and Agisoft Metashape Professional Course logistics webpage for links to software in case you don't have it installed already.

Agisoft Metashape part

Processing multispectral drone imagery

Add photos to chunks

Click Workflow->add folder-> RGB folder
In the Workspace pane right click on Workspace and select “New chunk” from the dropdown menu. Double click at the new chunk in the Workspace pane (it will appear bold when active) 
Click Workflow->add folder-> RED folder
Rename your Chunks by right-clicking on Chunk 1 > Rename. Change the name of the RGB chunk (Chunk 1) to “RGB” and the RED chunk (Chunk 2) to “RED”. Each chunk includes a sub-category Cameras, which contains a list of the loaded images. The chunks contain photos taken during the site survey as well as images captured during the ascent and the descent of the drone. You do not need the images from the ascent and descent

Camera calibration

Note: this step is only completed for the RED dataset: 
Click Tools > Camera calibration
to correct Black level. This setting adjusts the color contrast in the images. On the left side of the dialogue, choose the Sequoia Camera (if it is not already selected). Go to the main panel and choose the Bands tab. Change the Black level from 4819.9 to 0, and click OK.

Aligning photos

To exclude images from the ascent and descent, you can manually select them from the list using the image numbers below and then right-click Photos > Disable Camera. You may note that some additional photos are disabled. These photos have also been excluded from the alignment.
  • For RGB: select images 0005 to 0065 (ascent) and images 0180 to 0206 (descent).
  • For RED: select images 0004 to 0022 (ascent) and images 0180 to 0206 (descent). 
Click on Workflow > Align Photos
to align the photos. (you need to align the RED chunk and RGB chunk separately by activating each one of them)
  • Accuracy: RGB: Low, RED: Medium
  • Generic Preselection: Yes
  • Reference Preselection: Yes (drop-down at the right: source)
  • Key point limit: 40,000
  • Tie point limit: 4,000
  • Adaptive camera model fitting (Advanced Settings): Yes

Optimizing sparse point cloud

The result of the alignment process is a sparse point cloud, which contains tie points of the aligned photos. Tie points are locations or objects that are visually recognizable in the overlap area between two or more images. The tie points contain three different measures of accuracy: reprojection error, reconstruction uncertainty, and projection accuracy. Detailed descriptions of these three parameters are provided in the Agisoft user manual. For the purposes of this exercise, it is sufficient to know that Agisoft is able to calculate these three accuracy parameters, and you can optimize the sparse cloud by using the Gradual Selection tool to identify and delete points with high error according to the values of the accuracy parameters. You want to delete points with low Reconstruction Uncertainty values before moving on to the next steps. To do this: 
Click Model > Gradual Selection
to optimize the sparse cloud. In the dropdown menu, choose Reconstruction Uncertainty. You can decrease or increase the uncertainty value of the points by sliding the button from left to right. After clicking OK, you can zoom in and out of the point cloud to inspect the points and see which points have been highlighted as you change the level of reconstruction uncertainty. Set the levels as indicate below then click OK.:
  • RGB: 16.0
  • RED: 44.0
Click Delete on the keyboard to remove the selected points. The last step is to adjust the camera positions and correct the point cloud. 
Click Tools > Optimize Cameras.
Leave all settings as default and click OK.

Building a 3D Model (Mesh)

In this exercise, you will use the optimized sparse point cloud to build a 3D model (Mesh). A mesh is a 3D structure on which the photos are fitted to create an orthomosaic. You can imagine this process as putting paper mâché onto a wire model. The process works best when the surface of the mesh is smooth, without hard edges or gaps.
To develop a mesh, 
click Workflow > Build Mesh.
Choose the following settings for a maximum number of faces (geometric parts forming the mesh surface), and then click OK.
  • Source data: Sparse cloud
  • Surface type: Arbitrary (3D)
  • Face count: Custom > 0
  • Interpolation: Enabled (default)
  • Calculate vortex colors (Advanced settings): Yes
To create a smooth mesh, 
click Tools > Mesh > Smooth Mesh
and input the following values, then click OK.
  • Strength: 15
  • Fix borders (Advanced settings): Yes

Generating an Orthomosaic

You will now generate an orthomosaic using the mesh you created above. NOTE: This step works only if you have the licensed version of Agisoft, it in unavailable in the DEMO mode. If you are working in a DEMO mode, the outputs will be provided for you to download. 
Workflow > Build Orthomosaic
and input the below settings, then click OK.
  • Projection type: Geographic
  • Projection: WGS 84 / UTM 32N (EPSG::32632)
  • Surface: Mesh
  • Blending mode: Mosaic
  • Enable hole filling: Yes
  • Pixel size: leave default setting
To export your orthomosaic as a GeoTIFF file, 
click File > Export > Export Orthomosaic> Export JPEG/TIFF/PNG,
name the files as follows and save them.
  • RGB: “orthophoto_RGB”
  • RED: “orthophoto_RED”
Use the following settings:
  • Projection type: Geographic
  • Project: WGS 84 / UTM 32N (EPSG::32632)
  • Pixel size: leave default settings
  • Output file: Write big TIFF file

GRASS GIS part

Generating spectral indices from UAS-derived orthomosaics

The RGB sensor on the Parrot Sequoia is equipped with a rolling shutter while the multispectral sensors use a global shutter. These differences can impact the registration of the different bands, and cause geometric distortions. Additionally, the four spectral sensors are slightly offset from each other on the Parrot Sequoia, which leads to band misregistration and can affect the calculation of spectral indices by preventing the bands from being combined into an orthomosaic. To achieve optimal results, it is necessary to perform an extra step to align the bands when preprocessing spectral orthomosaics.

Start Grass GIS Create a new location based on EPSG: 32632 (do not apply transformation)
Open the PERMANENT mapset in this location

Change the current working directory to the directory where you downloaded the LAS files using cd command and path or in case you work in command line in GRASS GUI just type cd and press enter and select the directory using a dialog. 

cd ~/Downloads
set the region 
g.region n=5612389.15971435 s=5612164.15046685 e=466168.04885871 w=465909.24655539 res=0.1 -p
Import the RGB orthophoto into GRASS 
r.in.gdal input=ortho_RGB.tif output=orthophoto_RGB
Import the multispectral layers 
r.in.gdal input=ortho_GRE.tif output=orthophoto_GRE
r.in.gdal input=ortho_RED.tif output=orthophoto_RED
r.in.gdal input=ortho_NIR.tif output=orthophoto_NIR
r.in.gdal input=ortho_REG.tif output=orthophoto_REG

Working with the RGB Orthomosaic

You will use the orthomosaics to compute vegetation indices. You will first work with the RGB orthomosaic before moving on to multiple single-band orthomosaics in the next part. You have downloaded the orthophoto with three bands combined (red, green and blue). We need to separate the bands to use them in the indices calculation 
g.region n=5612389.15971435 s=5612164.15046685 e=466168.04885871 w=465909.24655539 res=0.1 -p
r.rgb input=orthophoto_RGB@PERMANENT red=ortho.red green=ortho.green blue=ortho.blue
calculate VARI: Visible Atmospherically Resistant Index 
i.vi output=VARI viname=vari red=ortho.red@PERMANENT green=ortho.green@PERMANENT blue=ortho.blue@PERMANENT
r.colors -e map=VARI@PERMANENT color=bgyr 
d.legend -f -s -d rast=VARI
r.univar VARI
Calculate RI: Redness Index 
r.mapcalc expression="RI = 1.0 * (pow ( ortho.red@PERMANENT, 2 ) )/( ortho.blue@PERMANENT *pow( ortho.green@PERMANENT ,2))"
r.colors -e map=RI@PERMANENT color=bcyr
d.legend -f -s -d rast=RI
r.univar RI
Calculate BI: Brightness Index 
r.mapcalc expression="BI = 1.0 * ( sqrt(pow( ortho.red@PERMANENT ,2)+pow( ortho.green@PERMANENT ,2)+pow( ortho.blue@PERMANENT ,2)) )/3"
d.legend -f -s -d rast=BI
r.univar BI
r.colors -e map=BI@PERMANENT color=byg

Working with multispectral orthomosaic

Calculate NDVI: Normalized Difference Vegetation Index 
i.vi output=NDVI viname=ndvi red=orthophoto_RED@PERMANENT nir=orthophoto_NIR@PERMANENT
d.legend -f -s -d rast=NDVI
or 
r.mapcalc expression="NDVI2 = 1.0 * ( orthophoto_NIR@PERMANENT - orthophoto_RED@PERMANENT )/( orthophoto_NIR@PERMANENT + orthophoto_RED@PERMANENT )" 
r.colors map=NDVI2@PERMANENT color=ndvi
d.legend -f -s -d rast=NDVI2

r.univar NDVI
r.univar NDVI2
Calculate IPVI: Infrared Percentage Vegetation Index 
i.vi output=IPVI viname=ipvi red=orthophoto_RED@PERMANENT nir=orthophoto_NIR@PERMANENT
r.colors map=IPVI@PERMANENT color=bcyr
d.legend -f -s -d rast=IPVI
r.univar IPVI
Calculate PVI: Perpendicular Vegetation Index 
i.vi output=PVI viname=pvi red=orthophoto_RED@PERMANENT nir=orthophoto_NIR@PERMANENT
r.colors map=PVI@PERMANENT color=bcyr
d.legend -f -s -d rast=PVI
r.univar PVI
Calculate SAVI: Soil-adjusted Vegetation Index 
i.vi output=SAVI viname=savi red=orthophoto_RED@PERMANENT nir=orthophoto_NIR@PERMANENT
r.colors map=SAVI@PERMANENT color=byr
d.legend -f -s -d rast=SAVI
r.univar SAVI
Calculate DVI: Difference Vegetation Index 
i.vi output=DVI viname=dvi red=orthophoto_RED@PERMANENT nir=orthophoto_NIR@PERMANENT
r.colors map=DVI@PERMANENT color=byr
d.legend -f -s -d rast=DVI
r.univar DVI
Calculate NDRE: Normalized Difference Red Edge Index 
r.mapcalc expression="NDRE = 1.0 * ( orthophoto_NIR@PERMANENT - orthophoto_REG@PERMANENT )/( orthophoto_NIR@PERMANENT + orthophoto_REG@PERMANENT )" 
r.colors -e map=NDRE@PERMANENT color=bgyr 
d.legend -f -s -d rast=NDRE
r.univar NDRE
Refer to lecture, supplemental materials and i.vi help page to interpret the results of the calculation of above indices.