Introduction

I am following the LiDAR project outline provided by Prof. Fouad. For this stage of my project the workflow I followed was informed mainly by the description of the project in the word doc as well as the lidR documentation. I first selected a area to study. Then I downloaded the relevant .las files. Next I read them into r, filtered for last returns and performed a validation step. Then I plotted the points to get an idea of what was going on. I used a TIN method to create a DTM as suggested by Prof. Fouad and following the instructions provided in the lidR documentation. My final step was to plot the newly created DTM.

For my study area I chose a section of the marsh in Jamaica Bay. Specifically, I chose the section of marsh that extends off the end of runway 4L/22R at JFK. Google maps listed the name of this area as JoCo Marsh so my variable names follow this convention. However, I’m not certain if that’s the name of this particular area.

Libraries

library(tidyverse)
library(janitor)
library(lidR)
library(rgl)
library(RCSF)
library(ggplot2)

Initial Steps

Here I set up my wd and read in a series of .las files I downloaded from https://orthos.dhses.ny.gov/# I chose the tiles that overlapped most of the marsh area. Some sections had mostly water so I excluded them. I use the “-last_only = TRUE” expression to read in only the last returns from the .las. I found out about this with the help of the documentation and an examination of the readLAS() function using ?readLAS. The reads for the 2017 data returned warnings for “withheld points” I’m not sure if this will be an issue going forward.

wd <- getwd()

# 2014 .las reads, these cover more area, inc only last return
joco_n_2014_las <- readLAS(file.path(wd, "Data", "18TXK015970.las"), filter = "-last_only = TRUE")
# joco_s_2014_las <- readLAS(file.path(wd, "Data", "18TXK015955.las"), filter = "-last_only = TRUE")

# 2017 .las reads, smaller than above, inc only last return
# joco_ne_2017_las <- readLAS(file.path(wd, "Data", "42165.las"), filter = "-last_only = TRUE")
# joco_nw_2017_las <- readLAS(file.path(wd, "Data", "40165.las"), filter = "-last_only = TRUE")
# joco_e_2017_las <- readLAS(file.path(wd, "Data", "42162.las"), filter = "-last_only = TRUE")
# joco_w_2017_las <- readLAS(file.path(wd, "Data", "40162.las"), filter = "-last_only = TRUE")
# joco_se_2017_las <- readLAS(file.path(wd, "Data", "42160.las"), filter = "-last_only = TRUE")
# joco_sw_2017_las <- readLAS(file.path(wd, "Data", "40160.las"), filter = "-last_only = TRUE")

Validation

For the validation step I followed the suggestion in the lidR documentation found at https://r-lidar.github.io/lidRbook/io.html#asprs-compliance. Note that I included all the read in objects but I comment them out because it was at this point that I decided to proceed with only the “joco_n_2014_las” object in order to examine a possible workflow in a simplified context.

Based on the validation it appears the .las has already been classified for ground points so I decided to skip the classification step but if there is a good way to double check this I would like to do that. Some other potential issues I noticed:

Duplicate points that could be removed. A missing crs which could impede mapping Large amount of negative outliers (ground below the water surface?) No normalization - is best practice to normalize if this shows up? Points with identical return number - not sure what this means

# validate the data (comment out as needed)
las_check(joco_n_2014_las)
## 
##  Checking the data
##   - Checking coordinates... ✓
##   - Checking coordinates type... ✓
##   - Checking coordinates range... ✓
##   - Checking coordinates quantization... ✓
##   - Checking attributes type... ✓
##   - Checking ReturnNumber validity... ✓
##   - Checking NumberOfReturns validity... ✓
##   - Checking ReturnNumber vs. NumberOfReturns... ✓
##   - Checking RGB validity... ✓
##   - Checking absence of NAs... ✓
##   - Checking duplicated points...
##     ⚠ 1535 points are duplicated and share XYZ coordinates with other points
##   - Checking degenerated ground points...
##     ⚠ There were 1 degenerated ground points. Some X Y Z coordinates were repeated
##   - Checking attribute population... ✓
##   - Checking gpstime incoherances
##     ✗ 3006158 pulses (points with the same gpstime) have points with identical ReturnNumber
##   - Checking flag attributes... ✓
##   - Checking user data attribute... ✓
##  Checking the header
##   - Checking header completeness... ✓
##   - Checking scale factor validity... ✓
##   - Checking point data format ID validity... ✓
##   - Checking extra bytes attributes validity... ✓
##   - Checking the bounding box validity... ✓
##   - Checking coordinate reference system... ✓
##  Checking header vs data adequacy
##   - Checking attributes vs. point format... ✓
##   - Checking header bbox vs. actual content... ✓
##   - Checking header number of points vs. actual content... ✓
##   - Checking header return number vs. actual content... ✓
##  Checking coordinate reference system...
##   - Checking if the CRS was understood by R...
##     ✗ EPSG code 32767 unknown
##  Checking preprocessing already done 
##   - Checking ground classification... yes
##   - Checking normalization... no
##   - Checking negative outliers...
##     ⚠ 2009641 points below 0
##   - Checking flightline classification... yes
##  Checking compression
##   - Checking attribute compression...
##    -  Synthetic_flag is compressed
##    -  Keypoint_flag is compressed
##    -  Withheld_flag is compressed
##    -  UserData is compressed
# las_check(joco_s_2014_las)
# 
# las_check(joco_ne_2017_las)
# las_check(joco_nw_2017_las)
# las_check(joco_e_2017_las)
# las_check(joco_w_2017_las)
# las_check(joco_se_2017_las)
# las_check(joco_sw_2017_las)

EDA

In this step I just plotted the points in “joco_n_2014_las” as expected they were very similar because they are all last returns. By manipulating the 3d plot you will be able to make out the area of the marsh as well as a section of runway 4L/22R.

# initial EDA
# plot(joco_n_2014_las)
# rglwidget()

DTM Creation

For this section I chose the TIN method because that was Prof. Fouad’s go to method. The lidR documentation provided a straightforward way to implement this using a built in algorithm. Once I created the TIN I visualized it with a plot.

# DTM creation
joco_n_2014_tin <- rasterize_terrain(joco_n_2014_las, res = 1, algorithm = tin())
## Delaunay rasterization[==============================--------------------] 61% (2 threads)
Delaunay rasterization[===============================-------------------] 62% (2 threads)
Delaunay rasterization[===============================-------------------] 63% (2 threads)
Delaunay rasterization[================================------------------] 64% (2 threads)
Delaunay rasterization[================================------------------] 65% (2 threads)
Delaunay rasterization[=================================-----------------] 66% (2 threads)
Delaunay rasterization[=================================-----------------] 67% (2 threads)
Delaunay rasterization[==================================----------------] 68% (2 threads)
Delaunay rasterization[==================================----------------] 69% (2 threads)
Delaunay rasterization[===================================---------------] 70% (2 threads)
Delaunay rasterization[===================================---------------] 71% (2 threads)
Delaunay rasterization[====================================--------------] 72% (2 threads)
Delaunay rasterization[====================================--------------] 73% (2 threads)
Delaunay rasterization[=====================================-------------] 74% (2 threads)
Delaunay rasterization[=====================================-------------] 75% (2 threads)
Delaunay rasterization[======================================------------] 76% (2 threads)
Delaunay rasterization[======================================------------] 77% (2 threads)
Delaunay rasterization[=======================================-----------] 78% (2 threads)
Delaunay rasterization[=======================================-----------] 79% (2 threads)
Delaunay rasterization[========================================----------] 80% (2 threads)
Delaunay rasterization[========================================----------] 81% (2 threads)
Delaunay rasterization[=========================================---------] 82% (2 threads)
Delaunay rasterization[=========================================---------] 83% (2 threads)
Delaunay rasterization[==========================================--------] 84% (2 threads)
Delaunay rasterization[==========================================--------] 85% (2 threads)
Delaunay rasterization[===========================================-------] 86% (2 threads)
Delaunay rasterization[===========================================-------] 87% (2 threads)
Delaunay rasterization[============================================------] 88% (2 threads)
Delaunay rasterization[============================================------] 89% (2 threads)
Delaunay rasterization[=============================================-----] 90% (2 threads)
Delaunay rasterization[=============================================-----] 91% (2 threads)
Delaunay rasterization[==============================================----] 92% (2 threads)
Delaunay rasterization[==============================================----] 93% (2 threads)
Delaunay rasterization[===============================================---] 94% (2 threads)
Delaunay rasterization[===============================================---] 95% (2 threads)
Delaunay rasterization[================================================--] 96% (2 threads)
Delaunay rasterization[================================================--] 97% (2 threads)
Delaunay rasterization[=================================================-] 98% (2 threads)
Delaunay rasterization[=================================================-] 99% (2 threads)
Delaunay rasterization[==================================================] 100% (2 threads)
# plot_dtm3d(joco_n_2014_tin, bg = "white")
# rglwidget()

Next Steps

Going forward, the project requires a DEM difference to be created. First off, what I created was called a DTM in the lidR documentation and I’m not sure this is the same as a DEM so I will need to confirm this before moving foward. The documentation suggests a normalization step follow the creation of the DTM so I will likely perform this. I also want to better understand the warnings that cropped up in the validation. I moved forward regardless in the interests of time but for the next stage I want to double check to make sure they are not impacting the quality of my analysis.

Assuming everything works out I will be moving forward with a difference using the corresponding tiles from 2017 that overlap with my 2014 tile. Once I have proven the process I will continue with the entire study area. Then I will create a map. At this point I would like to use the interactive map options we explored with Prof. Sun last section but this will depend on what time allows.