GIS595/MEA792: UAV/lidar Data Analytics

This is an unmaintained course material, please see current material at:

Fusion of uas and lidar data

Outline:

analyze the differences between the lidar and uas DEM/DSM
fuse lidar and uas DSM for viewshed analysis: compare simple patching with map algebra with reinterpolation from maximum elevation points

Data:

you should already have both the lidar and uas DSMs and points clouds.
Fields polygon as packed raster use r.unpack in your mapset to import it

Tools:

Your GRASS GIS 7 installation should also include libLAS library which is used by GRASS modules v.in.lidar and r.in.lidar (standalone GRASS GIS for MS Windows, OSGeo4W and Ubuntu packages contain libLAS).
libLAS installation should include command line tools lasinfo and las2txt.

Analyze the difference between the lidar and UAS DSM

First we need to evaluate the DSMs to find out which are suitable for fusion

compare errors at GCPs (lidar, agisoft, trimble or pix4d in june)
compute the difference between lidar DSM and selected uas: raster, profile/difference along a road and roof
use above to identify potential systematic error (you will see approx constant shift in profiles and most difs have the same sign, use mean or median as syst. error estimate, re-do the profile to verify the correction), use Mitasova et al. 2009 for reference
identify distortions - use corrected dsm and lid-uav diff, look for unexpected pattern, verify with profile. Put into report and exclude distorted areas from study area. Why are the distortions / errors low along the road but high in the fields? What would you have to do if you needed to correct the distorted area (add gcps using stable features in lidar and reprocess)
now you have two DSMs: corrected lidar DSM and (corrected) low-distortion subset of UAS DSM: create a new DSM for the entire area using at least two approaches and compare how the two approaches work for viewshed mapping

Suggested color table for differences, used in the workflow as dif_lidar_uav.txt:

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

Compute the differences between the lidar and uav DSMs

g.region rast=mid_pines_lidar2013_dsm -p
r.mapcalc "diff_lidardsm_agi_june = mid_pines_lidar2013_dsm - 2015_06_20_DSM_agi_11GCP"
r.colors diff_lidardsm_agi_june rules=dif_lidar_uav.txt
r.mapcalc "diff_lidardsm_pix4d_june = mid_pines_lidar2013_dsm - 2015_06_20_pix4d_11GCP_dsm"
r.colors diff_lidardsm_pix4d_june rules=dif_lidar_uav.txt

Evaluate possible systematic shift along stable features using the profile tool or the following commands

r.profile -g input=diff_lidardsm_agi_june output=road.txt coordinates=637208,219491,637059,219734
v.in.ascii input=road.txt output=profile_road separator=space columns="x double,y double,profile double,diff double"
v.univar map=profile_road column=diff

Because the shift in lidar is small we will ignore it in our first set of fusion examples, but we will use the analysis above top select the region with minimal distortion in the UAS DSM and compute clipped DSM

g.region n=219720.9 s=219257.1 w=636762.9 e=637190.4 res=0.3 -p
r.mapcalc "2015_06_20_DSM_agi_11GCP_cl = 2015_06_20_DSM_agi_11GCP"

Patch UAS DSM with lidar DSM to cover the entire Mid Pines study area and use shaded relief to check for patch edges

g.region rast=mid_pines_lidar2013_dsm -p
r.patch 2015_06_20_DSM_agi_11GCP_cl,mid_pines_lidar2013_dsm out=lid_uas_patch
r.relief lid_uas_patch out=lid_uas_patchsh

We have visible discontinuity along the patch edges so we will use map algebra to do the patching with averaged overlap. First we define a region for 6m (20 cell) overlap (you can do this also interactively using mouse)

g.region rast=2015_06_20_DSM_agi_11GCP_cl -p
g.region s=s-6 n=n+6 e=e+6 w=w-6 -p

Then, within this region, we compute a raster with lidar+uas DSM average in the overlap and UAS elsehwere

r.mapcalc "lid_uas_overlap_avg = if(isnull(2015_06_20_DSM_agi_11GCP_cl),(2015_06_20_DSM_agi_11GCP+mid_pines_lidar2013_dsm)/2.,2015_06_20_DSM_agi_11GCP_cl)"

We can now merge it with lidar for entire study area and check the patch edges using shaded relief

g.region rast=mid_pines_lidar2013_dsm
r.mapcalc "lidar_uas_mergedavg = if(isnull(lid_uas_overlap_avg), mid_pines_lidar2013_dsm, lid_uas_overlap_avg)"
r.relief lidar_uas_mergedavg out=lidar_uas_mergedavgsh

We now have smoother transition from UAS to lidar, but you can still identify two edges when viewing the patched DSM in 3D view (nviz). Therefore we will use weighted average to smooth out the overlap with the following workflow (try to explain it in your report):

g.region rast=2015_06_20_DSM_agi_11GCP_cl -p
g.region s=s+10 n=n-10 e=e-10 w=w+10 -p
r.mapcalc "dummy = 1"
g.region rast=2015_06_20_DSM_agi_11GCP_cl -p
r.grow.distance -m input=dummy distance=distance
g.region rast=mid_pines_lidar2013_dsm
r.mapcalc "patched_smooth = if(isnull(distance), mid_pines_lidar2013_dsm, (1 - distance/10.) * 2015_06_20_DSM_agi_11GCP_cl +(distance/10) * mid_pines_lidar2013_dsm)"

We will create a smooth transition when modeling the flow using a different method, in this assignment we will explore the viewshed analysis.

We would like to set up a webcam to monitor this area during the growing season so we need to update the lidar-based DSM with the UAS 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. You can use digitizer to create a polygon or get the packed raster representation of fields from the link above.

r.unpack -o fields.pack
r.mapcalc "uas_june_clpol = 2015_06_20_DSM_agi_11GCP * fields" 
r.patch uas_june_clpol,mid_pines_lidar2013_dsm out=lid_uas_patch_polygon 
r.relief lid_uas_patch_polygon out=lid_uas_patch_polygonsh

We compute viewshed from the point 637100,219360 using lidar only, uas only, and patched lidar+uas,

r.viewshed input=mid_pines_lidar2013_dsm output=viewshed_lidar coordinates=637100,219360
r.viewshed input=2015_06_20_DSM_agi_11GCP output=viewshed_uas coordinates=637100,219360
r.viewshed input=lid_uas_patch_polygon output=viewshed_liduaspatchpoly coordinates=637100,219360 --o
r.colors viewshed_lidar co=reds
r.colors viewshed_uas co=greens
r.colors viewshed_liduaspatchpoly co=blues

To compare the viewsheds display the resulting maps with transparency. Discuss the result.

Note: To find the location and height of the webcam that will capture the entire field during entire year we can run r.series on the time series of UAS DSMs to derive the DSM envelope surface and use it to analyze the viewshed. Of course, this assumes that the crops will be the same the following year.