This Weeks Lab
We are going to create two DEM layers this week using Lidar Data. We have been playing with different elevation models to utilize in combination with other remotes sensed data – usually for remote area. Today we are focusing on urban remote sensing (well as urban as Prince George is).
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, so as before – we will load data into QGIS for a quick peek.
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, and then make use of Open Street Map using the OpenLayers plugin (web pulldown menu) to add some back ground data.
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).
1. DEM data for the comparison – getting some raw LAS files
Open up a terminal 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 – but 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:
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 withing the 201 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
We are pros at this now – no need for explanation from Scott – right? 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 that is the way we go. 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 centre. 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?
Burning Vectors in for the GeoWall and Terrain Bender Assignment
You can use the data from the Quesnel Lake data set to test your vector burning skills. You should have a Landsat clipped layer and a vector lakes layer clipped to the same extents from that exercise.
You should scale the RGB composite you want to use in the GeoWall by using the scale function – creating new 8 bit channels. That way when you export to a tiff file later on – the data is already stretched.
Run POLY2RAS to rasterize the lakes layer you created into a channel in the scaled file.
EASI Modelling – of course
Run a script that represents your channel placements to fit the process below
if %lakeschannel = “the value of the ras2poly for the lakes” then
%channel holding the blue RGB values = 242endif
Do the same for the green and red channels using values 16 and 13 (for instance)
Look at your new combination in Focus. Then export as tiff with the channels in the right RGB order.