ESCI497Z:
UAS for Environmental Research and Monitoring
North
Fork Nooksack: Channel Elevation Change using LIDAR and UAS SFM
Objective:
Floodplain channel dynamics and channel migration have significant impacts on
salmon habitat. We will use both LIDAR data and UAS data to quantify floodplain
dynamics for a portion of the North Fork of the Nooksack River.
You
will be preparing a short (3-5 pages of text) write up describing your flight
plan.
PROCEDURES:
Step
1: Get the data. All of the data for this lab are located
in: J:/saldata/esci-497/N_Fork_Nooksack. This folder
is huge. Nevertheless, I’d suggest copying the ENTIRE
THING to your own folder on C:/temp. (At
the end of the day, I’d strongly urge you to delete
this stuff from C:/temp to make space for others.
Step
2: Review metadata (processing reports)
-LIDAR:
We
will use a LIDAR dataset from 2005. Take a look at the
processing report “2309 Final Report.pdf” for information about this data set.
-UAS: I did
flights over two days to cover this area. I have already processed all of the
images using Agisoft following the same general procedure that you used last
week from the heron colony. On Aug 30, 2016, I flew from our field off of Rutsatz Rd.. Take
a look at the processing report (Aug30_2016.pdf). These fights
covered the west half of the study area and included nearly 2400 images. In the
report, you can probably see that, in order to get sufficient overlap, I cam up with two different flight plans; one with flight
lines oriented at 360 degrees (due north) and the other at 300 degrees. I was
hoping that this would provide sufficient overlap to align the images over the
areas with tall trees. As you will see in working with the data, I still ended
up with big gaps; sections over forested areas where I
was unable to align the images. On Sept. 1, I flew from another site on the
north side of the river and about 2 km NE of our field along Rutsatz Rd.. This dataset
(described in Sept1_2016_flts1_5.pdf) includes over 4000 images. For these
flights, I increased the sidelap from 60% to 80%. I
was unable to reduce the flight speed to increase frontlap
so I simply flew most of the area twice; I repeated the first two legs of the flightplan. This is the portion of the
flight that mostly focuses on the south side of the river. This improved
the overlap somewhat and I did not have as many gaps.
The total area covered by the 2 days of flying covered
about 5 kilometers of river.
Step
3: View the data in ArcMap
We will mostly be working in ArcMap today. As I
mentioned above, I’ve done the Agisoft processing for
you. Doing so (with the ~6400 images) took days of processing. As you know, the
two most time consuming steps in Agisoft involved the photo alignment and the
generation of the dense point cloud. I typically start one of these tasks as I
left my office in the evening and the job was usually complete by the next
morning.
Start by opening ArcMap. Add all of your data layers
by clicking on the Add data
button. From the Add Data dialog box,
you may need to click on the Connect to a folder
button to navigate to your workspace on C:\temp. I am working in
C:\temp\wallin|\n_fork_nooksack. Do not try to work off of
the J-drive or off of your U-drive. You will not have enough space!!!!!!!! You should
see seven .tif files. Go ahead and add all of them to
your ArcMap project. These files include:
UAS:
all
with UTM projection
Aug30_2016_ortho_big_utm.tif (ortho of west half of
study area; 4 cm pixel size)
Sept1_2016_flts1_5_ortho_utm.tif (ortho of east half
of study area; 4 cm pixel size)
Sept1_2016_flts1_5_dem_utm.tif (DEM including both E
and W half of study area; 16 cm pixel size)
LIDAR:
all
with 2 meter pixel size and UTM projection
1st
return (top of canopy)
Nook2005_dem_2m.tif
Nook2005_dem2m_shade.tif (shaded relief)
Last
return (bare earth)
Nook2005_dtm_2m.tif
Nook200r_dtm2m_shade.tif (shaded relief)
We do not have an ortho image of the area at the time
of the LIDAR flight but you can look at Google Earth to find something close.
And, the shaded relief DEM (Nook2005_dem2m_shade) does help you visualize where
the trees are at this time.
Click the various layers on and off to get a sense of
what has changed over this 11 year period. You will note a great deal of
channel migration and changes in riparian vegetation.
Step
4: Elevation change
Click on the Search
button. In the search box, type in “Raster calculator.” Select Raster Calculator (Spatial Analyst).
This is probably the second item in the list. This will bring up the Raster Calculator window:
From here, things are pretty simple. Double-click on
the Sept1_2016_flts1_5_dem_utm.tif”, then click on the minus sign “-“, then
double-click on the Nook2005_dem_2m.tif layer. Then specify a location for your
output raster. I would strongly suggest creating a new folder in your workspace for
your output. When you wrap up today, you will want to save all of your output
but probably NOT all of the huge files that you are starting with. I created a
“nook_out” folder in my workspace, like this: C:\temp\wallin|\n_fork_nooksack\nook_out
Note that Arc limits you to filenames of no more than
13 characters:
After specifying your output raster, click OK in the Raster Calculator window.
It will take a few seconds to complete this
calculation. The resulting output file will have a 2 meter pixel size and it
will be added to your Table of Contents in Arc. My result looks like this:
This grey-scale image is rather underwhelming and hard
to interpret. So, right-click on your difference image name and to Properties-Symbology. Click on Classified. Say “Yes” to calculating a
histogram. Select a color ramp that pleases you, then go to Classify and select some reasonable
class breaks. The default here is 5 classes but you can increase this if you would
like. Then click on OK-Apply. Here
is one possible way to symbolize the elevation changes:
In this case, the warm colors represent elevation LOSS
between 2016 and 2005 (2005 elevations higher) and the cool colors represent
elevation GAIN over this interval. All elevations are in meters.
Reality-check:
Before diving too far into interpretation, let’s do a
reality-check. Take a look at the Ortho image from 2016. What feature(s) in the
study area would you expect NOT to change in elevation over this 11 year
interval?
Think real hard……..
How about the roads? Yes, some of the roads might have
been resurfaced but this would only result in a change of a centimeter or two.
And yes, there might be a few places where there might be some
slumping/subsidence but this is likely to occur in only a handful of places.
So, let’s take a look. Zoom in on a section of road
(with the 2016 ortho layer activated) and use the Identify
tool to get the elevation difference for this location. Record this value by entering this
into a spreadsheet. Zoom out and
move to another part of the study area and try again. Repeat this for a few
dozen locations that are well distributed around the entire study area and then
calculate the mean, maximum and minimum differences. From this, you should see
a systematic bias. One of the layers seems to be a bit higher than the other.
You can use your mean difference from your road points for bias correction (below).
But first, which elevation layer is “wrong?”
GCP
data:
One the things that you should have copied along with
all of the other data is a folder called GCP_data.
This includes a shape file called “point.ge.” Add this to your Arc project.
This will add about 14 points, scattered across the study area. For the ones
that are located on roads, calculate the difference between the GPS elevation
(GNSS height) and the elevation on the 2016 and 2005 DEM elevations. Based on
these results, which layer seems to most closely align with the GPS elevations?
Bias
correction:
Go back to the Raster
calculator “correct” your 2016-2005 difference layer by adding or
subtracting the mean difference from your road points. After doing so, redo
your symbology.
Interpretation:
Zoom in on some of the big losses and gains and see if
you can figure out what is causing them. To help you interpret this, toggle
between the 20016-2005 elevation difference, the 2005 shaded relief, and the
2016 ortho images.
What are the drivers of these big (+ or – 10-20 m)
losses and gains?
What is the driver of the more modest (+ or – 1-2 m)
losses and gains?
Note that the maximum elevation loss is 100 m and the
maximum elevation gain is 52 m. Do these max/min values seem
reasonable/realistic? Are these values “real” or due to some anomaly in the
data?
Frequency
Distribution of Change in Channel Migration Zone:
Some of the area in our 2016-2005 difference image is
“upland” that is not within the channel migration zone. We would
like to quantify the amount and percentage of the area within the channel migration
zone has undergone different types of change. You can visualize the channel
migration zone by viewing the LIDAR-based bare earth shaded relief layer
(nook2005_dtm2m_shade.tif). Using this, I manually created a polygon that
approximately delineates the channel migration zone. Add the shapefile
channel.shp to your project (click the button
again and point to “channel” in the same folder with all of your other input
layers).
We now want to clip our 2016-2005 difference image
using the “channel” polygon. To do this, open the Arc Toolbox and
go to Data Management
tools-Raster-Raster Processing-Clip
From the Clip
dialog box, your Input raster is
your 2016-2005 difference layer, the Output
extent is the channel shapefile. It is critically important that you check
the use Input Features for Clipping Geometry box! This is what forces
your input layer to be clipped to the exact boundary of our channel shapefile.
Specify an Output Raster Dataset
and, as before, I’d strongly urge you to put this in the folder with your
2016-2005 difference layer. Click OK to
run the clipping procedure. The output will be added to your project.
Note that this output layer consists of real numbers
(decimal values) defining elevation change within the channel migration zone.
To simplify this, we’d like to convert these to integer values.
To do this, in the ArcToolbox, go to Spatial
Analyst tools-Math-Int
In the Int
dialog, simply specify your input layer and a name for the output layer (again,
put this in the folder with your other output files). You can then export the
attribute table, bring it into Excel and generate a frequency distribution.
Report:
As with last week’s lab, I would like you to prepare a
brief report based on your work today.