This Weeks Lab
- DEM Differencing
- Ortho Rectification of Rapideye imagery
We are going to create two DEM layers this week using Lidar Data. We have been playing with different elevation models in combination with other remotes sensed data – usually for remote area (i.e. Topographic Effect). Today we are focusing on urban remote sensing – sorry back in Prince George.
We have data from 2010 and 2014 to work with. Most of you have been at UNBC for that time, so (as you are all very in tune with your city) you should recall changes that occurred downtown.
The data are in the /home/data directory (not in the labs directory) and for a change – we will start by using QGIS for a quick peek at data.
Open QGIS and build a shape file index for the images in /home/data/pg_city_2010_lidar/ProjectAOI/Images (raster –> miscellaneous –> create index
Now do the same for /home/data/pg_lidar_ortho_2014/ortho2014
Make sure the project is set to the UTM 10 projection (should be by default if the shape files you created are loaded into the project), and then make use of Open Street Map using the OpenLayers plugin
- Plugins –> Manage and Install Plugins
- Setting tab in Plugins menu –> click on all three radio buttons
- All tab –> Search openlayers –> click to select OpeLayers Plugin –> install (bottom right) –> close plugin menu
- Web menu –> OpenLayers –> Open Street map
We want to work with the sheets that cover the George street section of downtown.
We do not actually have imagery for 2010 (although there are images for 2009 in the 2009 directory of /home/labs/geog457/dem_differencing), but we can see how the imagery (index grids) overlap. We will be using the LAS layers for 093G09712(3). You can load the appropriate image tile into QGIS by adding raster layer into QGIS (layer –> load raster layers).
1. DEM data for the comparison – getting some LAS files
LAS file definition from ASPRS – https://www.asprs.org/committee-general/laser-las-file-format-exchange-activities.html
Open up a terminal on the Linux workstation and get prepared for some typing! We can see from the vector index layers that the data more or less share common boundaries. We need to create and area that is common for both – and we mean the exact ground is covered. Scott beat his head on this for a few hours trying to align data without LASTools and using only PCI - more hair lost for no reason…
Question: What limitation does VDEMINT have in regards to setting up an exact resolution and extent file for comparison?
Finding our extents
Create a folder to work in (I created one called dem-differencing).
Issue the command
We get data that include the following lines:
Minimum and Maximum Attributes (min,max)
Min X, Y, Z: 516416.83, 5973815.73, 563.91
Max X, Y, Z: 518063.84, 5975213.19, 687.65
Bounding Box: 516416.83, 5973815.73, 518063.84, 5975213.19
Time: 0.000000, 0.000000
Return Number: 1, 7
Return Count: 1, 7
Flightline Edge: 0, 1
Intensity: 9983, 65535
Scan Direction Flag: 0, 0
Scan Angle Rank: 0, 0
Classification: 1, 9
Point Source Id: 49, 60
User Data: 0, 0
Minimum Color (RGB): 0 0 0
Maximum Color (RGB): 0 0 0
We need to create an extent to work with by subsetting the 2010 data to match the above data set.
Issue the command (use you own user name – not Scott’s of course):
las2las -i /home/data/pg_city_2010_lidar/ProjectAOI/Raw/NonClassified/LAS/093g09712_raw.las -o /home/emmons/geog457/dem_difference/2010_las_subset.las -e “516420.00 5973820.00 518000.00 5975210.00″
Do the same for the 2014 data using the exact extents.
Why do we subset the 2014 – could we not just use the data extents as the LAS data is within the 2010 data?
View the LAS data
wine /home/lastools/bin/lasview “the path to your 2010 data”
Do the same for the 2014 data
What is different between the datasets?
2. Getting the data into PCI
We need to launch the EASI command interface from the Geomatica tool set (not focus). Start it up…
We see we need to set some parameters
points = “/path to your subsetted file for 2010.las data
filv =”/path to a new pix file for the subsetted 2010 las data
If it looks good
Now do the same for the 2014 data creating a second pix file for the subsetted 2014 las data.
3. Create DEM a for each year from points in the two pix file vector segments
Use the algorithm called VDEMINT. We are pros at using these algorithms – no need for explanation from Scott – right? For an explanation of how we use VDEMINT – check out the creation of DEMs in PCI lab:
We want the resolution to be 1 metre and use 1024 as the memory setting.
Check out the resolution and extents to make sure they will both line up – you should look at them in focus (as if you could resist!)
4. View the DEMs superimposed and make one file for modelling
Layer -> add -> RGB and load channel 1 (2010) in red and channel 2 (2014) in both green and blue. Red areas will show where the 2014 DEM is falsely low, and blue-green areas where it is falsely high.
One file fits all
Either use one of the existing DEM pix files, or create a new one (probably easier to manage). Add both channels into the new file.
What could we do to reduce the size of these files by creating new channels and moving the two DEM layers into it?
5. Generate the difference channel and display in pseudocolour
You can do this in two ways: using either operation ARI or tools-> EASI modeling. For the former, you would direct the result to your file and it will generate the new channel, or if you use EASI modeling, you need to first create an empty 32 bit (or 16 bit if you managed the data differently) channel
I (Scott) prefer using EASI – I cannot recall even using ARI, so guess that means you are on your own to experiment. Subtract one channel from the other, lets try 2010 – 2014.
6. View the histogram for the created difference channel to understand the values created.
What do you see? True changes are likely near the center. You can restrict the histogram by using sensible values that you know could occur in the real world – 100 metre changes are not realistic.
Start a new project and open the file directly into Focus (as a RGB and look at the data). Pretty easy to see the changes.
Load the difference channel into Focus as a pseudo colour channel.
Edit the colours, by:
Right-click on the layer in the files tree. Select edit pseudocolours
Select range-based instead of single value (tabs) and play with the setting to reveal the data you are looking for. He is an example Scott worked up:
You can turn off some of the classes to see the changes in the landscape. You can load the 2014 and 2009 imagery to see the changes under the PCT layer.
7. Calculating change – Questions
You can do volume and area changes in the scene, but we are really looking for what happened downtown. Well … what happened.
8. Bitmap creation for volume calculation
We will use the reliability layer to create a bitmap which identifies all pixels with ‘0’ reliability (= clouds or other effects), and thus to help isolate areas of real change. We assume that non-zero values are meaningful. We might also exclude the values of 104 and 108.
You can add a new empty bitmap – in the files tree on left, right click and add one bitmap, note its number, then in easi-modelling:
If %diff channel > “some realistic number” and %x < “another realistic number” then
%%y = 1
(where x is the reliability channel number, and y is the new bitmap)
You can also add a new channel and do the same thing, but modelling under a bitmap works much better.
Display the result visually to check -
7. Volume (task VLM)
We can calculate the volume of change, but would need to exclude the erroneous DEM values from cars, trees, clouds etc
Question – explain how you decided upon your base value?
Rectification for Rapideye imagery
We are going to have a look at GeoRectification of data from the high resolution sensor RapidEye. We are lucky with this lab in that we can rely on improvements in sensor technologies to rectify the data at hand. Do you recall the methods we employed in Geog 300 to rectify an aerial photo? Perhaps we should practice this technique using PCI?!
Roger has created a PCI image from the Rapid Eye data we received for Prince George. You can see the original downloads as well as the pix file in the /home/labs/geog457/lab4-rectification/rapid_eye_pg_area directory.
Points to think about:
- What is the file format the files come in
- What was the process that the files were loaded into the pix file – we should all be pros at this by now
- What are the other non raster layers (segments in the file).
Look at the original data
Load the pix file in the /home/labs/geog457/lab4-rectification/rapid_eye_pg_area directory using the file and north to the top options. Now load the pg boundary shape file and the DEM raster layer in the /home/labs/geog457/lab4-rectification directory. How do things line up?
Now load the subset of the Rapid Eye data Scott created from the DEM extents into the project (pg2013re_city_bound file). How does this layer look? What is going on here?
As well as the in being much easier than what we did in Geog300, we are going to include a DEM to OrthoRectify our imagery (to compensate for changes in relief). You need to open up OrthoEngine and move through the steps of creating a project, loading data and completing the orthocorrection by:
- Create a project (new project in your lab4-rectification directory)
- Pick the proper math model – you have to check out the options under the radio buttons
- Set the projection to UTM Zone 10 NAD 83 (pick the most appropriate projection from the list)
- Data input – we only have one uncorrected image to load – load the pg2013re_city_bound.pix file
- We will skip all the other steps (as we are not using any GCPs or any other layers)
- Select Ortho generation –> schedule a generation –> =make sure you pick a DEM (select the channel when it loads)
Your panel should look something like the image below (click it):
Then run the ortho rectification and look at the results in Focus.
Why is everything so bizarre looking? Lets discuss the steps that led to the output.
Now Try it again but use the eye_ball pix file Scott created. Lets discuss what happened this time.
Look Through the Data you Created and run through some of the tasks Roger was thinking we should look – see the following notes:
- 321 composite (natural colour)
- 532 (infra-red)
- 543 (NIR -Red Edge – Red)
Check the DNs in the display
- most features seem to have DNs that either all increase or all decrease through the 3 bands
- e.g. water – low DN decreasing 3-4-5
- forest: high DNs increasing 3-4-5 Red edge alone in grayscale
- which features have the highest values ? (soccer fields?)
- what is correlation between 5 and 4, tools -> scatterplot
- 4 and 3 (how is red edge correlated with NIR and Red respectively … which is it closer to ?
- View histograms: 1. visible bands are single-peak RE and NIR are bi-modal
- more pronounced in NIR (why)?
- Note minimum value in each band – it progressively gets lower. WHY ?
- It is because of haze – haze increases the DNs in blue, less so in green and so on as wavelength increases Band ratios
Ratios – in raster calculator run these ratios:
- in display only, and to 32 bit display to see the decimal values 5/4
- 5/3 (traditional vegetation ratio)
- 4/3 (RE/Red) they all seem to be trimodal – what are the 3 peaks ? (Roger thinks they are water, unvegetated e.g. urban, and vegetated – remember the 5, 4,3 composite/stretch)