LiDAR data
Different areas of Saguaro NP Rincon Mountain District were flown with airborne LiDAR between the years 2005, 2008, and 2011. In areas that were flown twice (or three times in the case of Cienega Creek) we can do difference maps to look for change over time.
All of the delivered LiDAR (*.LAS) data are projected in horizontal datum NAD83 (93 HARN), and vertical datum NAVD88. The Projection is Arizona State Plane (Zone: Central) and the units are in International Feet (ESPG:2223).
Digital Terrain Model Creation
In USFS FUSION (McGaughey 2013) I batch scripted the GridSurfaceCreate command to generate bare earth layers for all the PAG 2005, 2008, and 2011 LiDAR tiles (*.LAS) using classified bare earth returns (class = 2) only.
GridSurfaceCreate /median:3 /smooth:3 /class:2 J:\SAGU\DTM\%1.dtm 9.8425 f f 2 0 2 2 %1.las
Complete description of the GridSurfaceCreate parameters are given in the FUSION manual. I compared surfaces produced at 2 feet (Swetnam et al. 2013) re-sampled to 9.8425 feet or 3 meter resolution (using MergeDTM) and surfaces produced at 9.8425 feet directly by GridSurfaceCreate. I found the 3x3 nearest neighbor median and smooth switch caused blurring of erosion and deposition edges along stream channels. The 2ft merged to 9.8425ft surface did not have this issue because the 3x3 filter ran on a 6x6 foot area retained higher definition edges.
The second batch file called the process file calls a filelist.txt of the tile names and drive locations, e.g.
for /F "eol=; tokens=1* delims=,. " %%i in (J:\PAG2011\filelist.txt) do call process_tile %%i
In order to keep file sizes under 2 Gb DTMs were split north to south along township (36 sq mile areas) IDs , e.g. 10S, 11S, ..., 20S.
For the Rincon Basin the MergeDTM command on the 2ft DEM files generated for . The tiles were resized from 2ft to 3m pixels in the MergeDTM command line:
MergeDTM /cellsize:9.8425 /overlap:average /disk /verbose /precision:3 J:\SAGU\DEM\rincon_crk_3m.dtm I:\PAG2011\DTM\rinc_crk_08_11.txt
After the layer was generated the units were converted to meters using the Raster Calculator function in QGIS
The layer was then reprojected into UTM coordinates in QGIS using a GDAL Warp command:
gdalwarp -overwrite -s_srs EPSG:2223 -t_srs EPSG:2956 -tr 3 3 -r bilinear -of GTiff I:/SAGU/DEM/3m/rincon_crk_3m_stpln.tif J:/SAGU/DEM/3m/rincon_crk_3m_utm.tif
Intensity
Generated LiDAR derived intensity layers from the 2005, 2008, and 2011 LiDAR flights over Saguaro NP. Rescaled layers to 10m resolution with the r.resamp.stats command set to average.
The red laser LiDAR is in the range of 1064 nm, slightly longer than the NIR band for the NDVI. The LiDAR reflectance is also influenced by the type of surface the beam refracts off of. In general vegetation appears darker, having lower reflectivity, and bare rock appears brighter, having high reflectance.
The intensity is normalized to 8-bit 0-255 digital numbers (DN). The 2011 LiDAR had a characteristic threshold of 130-255 DN for vegetation, below 130 the brighter surfaces, particularly bare rock were readily identifiable from the aerial orthophotography.
The 2005 and 2008 LiDAR were much lower density than the 2011 flight and had large areas of Null (or no data) values; additionally the 2008 has several tiles where the intensity was not normalized correctly, making their values much darker than the rest of the flight area. Those areas were excluded from the workflow in later steps.
Post-analysis thoughts on intensity data are that they are not consistent enough, nor cover enough area to be of considerable advantage over the NAIP derived NDVI.
FUSION Gridmetrics
US Forest Service FUSION was used to create 10m gridded statistics of the LiDAR point clouds from all years (2005, 2008, 2011).
First, I generated the Gridmetrics using a list of 306 tiles from the three different flight years using a call batch command:
for /F "eol=; tokens=1* delims=,. " %%i in (J:\SAGU\FUSION\sagu_las_tiles.txt) do call gridmetrics_sagu %%i
and the execute batch command gridmetrics_sagu.bat:
GridMetrics /outlier:-3.2808,147.64 /ascii /minht:0 /raster:densitytotal,min,max,mean,mode,stddev,variance,cv,cover,abovemean,abovemode,skewness,kurtosis,AAD,p01,p05,p10,p20,p30,p40,p50,p60,p70,p80,p90,p95,iq,90m10,95m05,allcover,afcover,allcount,allabovemean,afabovemean,fcountmean,allcountmean,totalfirst,totalall J:\SAGU\DEM\3m\rincons.dtm 32.8084 32.8084 J:\SAGU\FUSION\GRIDMETRICS\10M\%~n1.csv %1.las
where the 3m DEM produced from the three LiDAR flights was used to normalize the height profile for vegetation. The GridMetrics produced 38 sets of elevation (height) metrics:
all_returns_all_metrics_elevation_aad.asc all_returns_all_metrics_elevation_all_above_mean_count.asc all_returns_all_metrics_elevation_all_above_mean_cover.asc all_returns_all_metrics_elevation_all_cover.asc all_returns_all_metrics_elevation_all_cover_count.asc all_returns_all_metrics_elevation_all_first_above_mean_cover.asc all_returns_all_metrics_elevation_all_first_cover.asc all_returns_all_metrics_elevation_count.asc all_returns_all_metrics_elevation_cover.asc all_returns_all_metrics_elevation_cover_above_mean.asc all_returns_all_metrics_elevation_cv.asc all_returns_all_metrics_elevation_density_count_total.asc all_returns_all_metrics_elevation_first_above_mean_count.asc all_returns_all_metrics_elevation_IQ.asc all_returns_all_metrics_elevation_kurtosis.asc all_returns_all_metrics_elevation_max.asc all_returns_all_metrics_elevation_mean.asc all_returns_all_metrics_elevation_min.asc all_returns_all_metrics_elevation_mode.asc all_returns_all_metrics_elevation_p01.asc all_returns_all_metrics_elevation_p05.asc all_returns_all_metrics_elevation_p10.asc all_returns_all_metrics_elevation_p20.asc all_returns_all_metrics_elevation_p30.asc all_returns_all_metrics_elevation_p40.asc all_returns_all_metrics_elevation_p50.asc all_returns_all_metrics_elevation_p60.asc all_returns_all_metrics_elevation_p70.asc all_returns_all_metrics_elevation_p80.asc all_returns_all_metrics_elevation_p90.asc all_returns_all_metrics_elevation_p90_minus_p10.asc all_returns_all_metrics_elevation_p95.asc all_returns_all_metrics_elevation_p95_minus_p05.asc all_returns_all_metrics_elevation_skewness.asc all_returns_all_metrics_elevation_stddev.asc all_returns_all_metrics_elevation_total_count.asc all_returns_all_metrics_elevation_total_first_count.asc all_returns_all_metrics_elevation_variance.asc
I then wrote a second command set which created 38 text files for each set of 306 tiles:
dir /b /s *all_returns_all_metrics_elevation_aad.asc >all_returns_all_metrics_elevation_aad.txt ... dir /b /s *all_returns_all_metrics_elevation_variance.asc >all_returns_all_metrics_elevation_variance.txt
I then ran the MergeRaster command from FUSION to mosaic all 306 tiles for each metric into a single file:
MergeRaster /overlap:max J:\SAGU\FUSION\GRIDMETRICS\10M\Merged\all_returns_all_metrics_elevation_aad.asc J:\SAGU\FUSION\GRIDMETRICS\10M\all_returns_all_metrics_elevation_aad.txt ... MergeRaster /overlap:max J:\SAGU\FUSION\GRIDMETRICS\10M\Merged\all_returns_all_metrics_elevation_variance.asc J:\SAGU\FUSION\GRIDMETRICS\10M\all_returns_all_metrics_elevation_variance.txt
Differenced elevation model (dDEM)
We are interested in modelling both pre- and post-fire hydrological impacts. As luck would have it, a meso-scale flood event that occurred over Southern Arizona in 2006 (Webb et al. 2008, Griffiths et al. 2009) happened between the 2005 and 2008 LiDAR flights conducted over the Tucson basin. The 2005 DEM (13-37 cm vertical uncertainty) was used as a pre-flood example and the 2008 DEM (18 cm vertical uncertainty) as the post-flood example. A DEM of difference, or dDEM was produced by subtracting the 2008 DEM from the 2005 DEM. There is significant evidence of post-flood impacts along Rincon Creek, Pantano Wash, and Tanque Verde Wash which are quantified by the dDEM that are also corroborated by aerial orthophotography from the same period of time.
The DEMs were produced in USFS FUSION using the DTM to DEM Creation information given above.