Error in example raster calc step if only using a single .las/.laz file

[akamoske]

Error in example raster calc step if only using a single .las/.laz file

[akamoske]

Fixed it!!!!

[akamoske]

OK, major updates. I just did a test, where I uploaded a .laz file that I have on my computer (that I've used a 1000 times for testing and just tested again) to figshare. Published it. Downloaded it into R using your code. And I get the exact same errors!

ERROR: 'chunk with index 0 of 1 is corrupt' after 2745 of 5276415 points   
Warning: there are 1474 points flagged 'withheld'.   
Warning: there are 1378 points flagged 'synthetic'.   
Error: cannot allocate vector of size 68662.6 Gb   

[akamoske]

Based on this and a trying a ton of different ways of downloading data, I think there is something funky with the way the R is downloading files via Windows. Something is causing the file to get corrupted during the download. I did some research on the the "chunk with index 0 or 1 is corrupt" issue with LiDAR data and it has something to do with the file itself when getting compressed. I'm going to try a couple of other ways to download the data and see if anything changes.

[akamoske]

ahh...

I know what the problem is. I just need to add a qualifying statement so that if there is only one file then it does not try to merge the list together.

I'll fix it and push it!

Good eye!

[serbinsh]

Aha! yes of course such that if z>1 then it merges, otherwise it ends the function call, or whatever

[akamoske]

Yup! Just tested it with 2 .laz files and then again with a single .laz file. Works good on my end!

Does it work for you?

I just pushed it to the main branch.

[akamoske]

I just added this code:

#if there is more than 1 raster in the list lets do this
  if(length(lad.rstack.list) > 1){

    #now that all the raster stacks are complete lets merge them together as one big happy raster
    lad.ras <- base::do.call(raster::merge, lad.rstack.list)

  } else {

    #set the only object in the list to this variable
    lad.ras <- lad.rstack.list[[1]]

  }

[serbinsh]

@akamoske It worked!! Nice one.

By the way, I created this simple example script which automatically downloads NEON lidar data from the NEON tutorial to use in your code. It may be helpful for your example since this way you dont need to ship example data with the package but instead use "the cloud" to provide example data to end users:

## Get missing R packages dependencies (if needed)
list.of.packages <- c("uuid","rlas","devtools")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[, "Package"])]
if (length(new.packages)) {
  print("installing : ")
  print(new.packages)
  install.packages(new.packages, repos = "http://cran.rstudio.com/", dependencies = TRUE)
}

## install LiDARforestR from github
if (!("LiDARforestR" %in% installed.packages()[, "Package"])) {
  devtools::install_github("akamoske/LiDARforestR")
}


## Setup scratch folder to contain datasets
scratch_folder <- file.path("~/scratch/neon_data/")
if (! file.exists(scratch_folder)) dir.create(scratch_folder,recursive=TRUE)
setwd(file.path(scratch_folder))
getwd()

## Download NEON example .las file
download.file(url = "https://ndownloader.figshare.com/files/7024955",destfile = file.path(scratch_folder,"neon_lidar_example.las"),
              method = "auto")

## Load LiDARforestR package
library(LiDARforestR)

## Convert .laz or .las files into a list of voxelized lidar arrays
laz.data <- laz.to.array(laz.files.path = file.path(scratch_folder), 
                         voxel.resolution = 10, 
                         z.resolution = 1)

## Level each voxelized array in the list to mimic a canopy height model
level.canopy <- canopy.height.levelr(lidar.array.list = laz.data)

## Estimate LAD for each voxel in leveled array in the list 
lad.estimates <- machorn.lad(leveld.lidar.array.list = level.canopy, 
                             voxel.height = 1, 
                             beer.lambert.constant = NULL)

## Convert the list of LAD arrays into a single raster stack
lad.raster <- lad.array.to.raster.stack(lad.array.list = lad.estimates, 
                                        laz.array.list = laz.data, 
                                        epsg.code = 32611)

## Create a single LAI raster from the LAD raster stack
lai.raster <- raster::calc(lad.raster, fun = sum, na.rm = TRUE)

## Generate a quick raster plot of the resulting total canopy LAI values for each pixel
plot(lai.raster)

I ran this on my desktop mac at home and everything worked!

[akamoske]

Cool! I just tried to run it and got this error.

> ## Convert .laz or .las files into a list of voxelized lidar arrays 

laz.data <- laz.to.array(laz.files.path = file.path(scratch_folder), 
                                     voxel.resolution = 10,                        
                                      z.resolution = 1) 

"Processing: C:/Users/Aaron Kamoske/Documents/scratch/neon_data/neon_lidar_example.las" 

WARNING: invalid bounding box [ 256000 4.111e+06 1.72046e+270 / 257000 4.112e+06 510.25 ]    WARNING: only 122 bytes until point block when trying to read 10375 bytes into header.vlrs[0].data    WARNING: only 0 bytes until point block after reading 1 of 2 vlrs. skipping remaining vlrs   
Warning: there are 1667416 points flagged 'withheld'.    
Warning: there are 1977164 points flagged 'synthetic'.    
Error: cannot allocate vector of size 68719.2 Gb   
--
 
> | >
>

[akamoske]

I'm going to look into it a bit right now, because this would be a cool part to include

[serbinsh]

What in the world?!?! That is a whacky error! Perhaps mac/PC differences in how the script needs to be written or we need to make a tweak to be PC/mac compliant? Did it actually download the las file to the scratch folder?

[akamoske]

Yeah it downloaded it into the scratch folder

[akamoske]

It seems like it is a problem with the NEON .las file. I'm on Windows 10. It won't let me read the .las file into R

[akamoske]

Hold on...I think I have an easy fix!

[akamoske]

Nope, nevermind. It is something wrong with the NEON .las file on Windows 10

[akamoske]

Here are some stats on the .las file. There is something funky going on...

scan angles range from: 0 to 114
classifications from 0 to 24
Intensity from to 7424

Also there are some negative x , y, z coordinates in there such as these: 1757822 6625956 -1508102

[serbinsh]

What in the world! (confirming again why I HATE Windows). So is it a file type error? What is the size of the file on Windows? ~180 mbs?

Really strange.....wait, perhaps it has something to do wtih the "method" in the download.file() function? Perhaps it is corrupting the file? Is there perhaps a setting that needs to be set on Windows?

[akamoske]

Yeah, WIndows sucks.

The file is ~180 mbs. I'm trying to do a bit of digging with it to see. I'll look into the download.file() function

[serbinsh]

https://stat.ethz.ch/R-manual/R-devel/library/utils/html/download.file.html

maybe try changing the method to method="wininet" ? That should happen automatically but maybe its not and it is using an incompatible method? Or perhaps it is something about .las on windows? Perhaps that is why .laz is better?

OK good luck! I feel bad for sending you down a rabbit hole :(

[akamoske]

Yeah, I'm trying the "wininet" method right now. I've used .las files a bunch on windows. All .las is is an uncompressed .laz file. So that shouldn't make a difference, I think at least... I've ran this code with .las files on windows and ubuntu before with no problems

[akamoske]

Just ran through all of the download methods and nothing changed. I'll keep digging into it though

[akamoske]

Not sure what the problem is. I even wrote some code to source a .laz file from the NEON API. I took a file from SERC that I have processed using this code (it is from the 2017 flights). And I get the exact same errors!!!! UGHHHH. Not sure what the deal is, considering that I followed their tutorial. Not sure what the deal is.

Here is the code:

## Get missing R packages dependencies (if needed)
list.of.packages <- c("uuid","rlas","devtools", "httr", "jsonlite",
                      "dplyr", "downloader", "geoNEON", "neonUtilities")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[, "Package"])]
if (length(new.packages)) {
  print("installing : ")
  print(new.packages)
  install.packages(new.packages, repos = "http://cran.rstudio.com/", dependencies = TRUE)
}

## install LiDARforestR from github
if (!("LiDARforestR" %in% installed.packages()[, "Package"])) {
  devtools::install_github("akamoske/LiDARforestR")
}

library(httr)
library(jsonlite)
library(dplyr)
library(downloader)


## Setup scratch folder to contain datasets
scratch_folder <- file.path("~/scratch/neon_data/")
if (! file.exists(scratch_folder)) dir.create(scratch_folder,recursive=TRUE)
setwd(file.path(scratch_folder))
getwd()

req.aop <- GET("http://data.neonscience.org/api/v0/products/DP1.30003.001")
avail.aop <- fromJSON(content(req.aop, as="text"), simplifyDataFrame=T, flatten=T)
cam.urls <- unlist(avail.aop$data$siteCodes$availableDataUrls)

cam <- GET(cam.urls[intersect(grep("SERC", cam.urls), grep("2017", cam.urls))][1])
cam.files <- fromJSON(content(cam, as="text"))

head(cam.files$data$files$name)

download(cam.files$data$files$url[grep("365000_4309000", cam.files$data$files$name)],
         paste(getwd(), "/SERC.laz", sep=""))

[serbinsh]

@akamoske huh?

The downloaded source packages are in
	‘/tmp/RtmpiXq1Ta/downloaded_packages’
Warning message:
packages ‘geoNEON’, ‘neonUtilities’ are not available (for R version 3.5.0)

[serbinsh]

FYI - https://github.com/NEONScience/NEON-geolocation/tree/master/geoNEON

library(devtools)
install_github('NEONScience/NEON-geolocation/geoNEON', dependencies=TRUE)
library(geoNEON)

library(devtools)
install_github("NEONScience/NEON-utilities/neonUtilities", dependencies=TRUE)
library(neonUtilities)

[serbinsh]

OK so on my linux cluster it seems like I was able to get all of the code to run:

Files download OK

drwxrwsr-x. 2 sserbin sserbin 66 Aug 11 08:52 neon_data
[sserbin@modex scratch]$ cd neon_data/
[sserbin@modex neon_data]$ ll
total 19376
-rw-rw-r--. 1 sserbin sserbin 19839606 Aug 11 08:52 SERC.laz
[sserbin@modex neon_data]$

LiDARforestR works:

> ## Convert .laz or .las files into a list of voxelized lidar arrays
> laz.data <- laz.to.array(laz.files.path = file.path(scratch_folder),
+                          voxel.resolution = 10,
+                          z.resolution = 1)
[1] "Processing: /home/sserbin/scratch/neon_data//SERC.laz"
  |======================================================================| 100%
> ## Level each voxelized array in the list to mimic a canopy height model
> level.canopy <- canopy.height.levelr(lidar.array.list = laz.data)
>
> ## Estimate LAD for each voxel in leveled array in the list
> lad.estimates <- machorn.lad(leveld.lidar.array.list = level.canopy,
+                              voxel.height = 1,
+                              beer.lambert.constant = NULL)
>
> ## Convert the list of LAD arrays into a single raster stack
> lad.raster <- lad.array.to.raster.stack(lad.array.list = lad.estimates,
+                                         laz.array.list = laz.data,
+                                         epsg.code = 32611)
[1] "Raster number 1 is complete!"
> ## Convert the list of LAD arrays into a single raster stack
> lad.raster <- lad.array.to.raster.stack(lad.array.list = lad.estimates,
+                                         laz.array.list = laz.data,
+                                         epsg.code = 32618)
[1] "Raster number 1 is complete!"
> ## Create a single LAI raster from the LAD raster stack
> lai.raster <- raster::calc(lad.raster, fun = sum, na.rm = TRUE)

>
> ## Generate a quick raster plot of the resulting total canopy LAI values for each pixel
> plot(lai.raster)
> q()

@akamoske you could try the RCurl package and use their downloader?

[akamoske]

download.file(url = "https://ndownloader.figshare.com/files/7024955",
              destfile = file.path(scratch_folder,"SERC_fig.laz"),
              method = "auto",
              mode = "wb")

[akamoske]

I think that on windows you have to explicitly tell it that to write as a binary file with the mode = wb command

[akamoske]

Worked on a .laz and a .las file

[akamoske]

Does using the mode = "wb" work on mac and linux?

[akamoske]

I updated the README file to show these changes. I'm not going to push it to the main until I hear back from you that it works on your machines though. Here is the README text though, so that you can see what I was running on my Windows machine. Thanks for all your help through this! Last hurdle before submitting that LAD paper!

LiDARforestR

R package to estimate leaf area density (LAD) and leaf area index (LAI) from airborne LiDAR point clouds.

Information

The theory behind this R package is described in:

Kamoske A.G., Dahlin K.M., Stark S.C., and Serbin S.P. 2018. Leaf area density from airborne LiDAR: Comparing sensors and resolutions in a forest ecosystem. In preparation for submission to Forest Ecology and Management.

Corresponding Author

Code written by: Aaron G. Kamoske, PhD Student

• Michigan State University, Department of Geography
• ERSAM Lab
• [email protected]

Contributing Authors

Scott C. Stark

• Michigan State University, Department of Forestry
• Tropical Forestry Ecology Lab
• [email protected]

Shawn P. Serbin

• Brookhaven National Laboratory, Environmental and Climate Sciences Department
• [email protected]

Installation

The easiest way to install LiDARforestR is via install_github from the devtools package:

# If you haven't already installed this package and its dependencies
install.packages("devtools")

# If you alread have devtools installed or just installed it
library(devtools)

# Install LiDARforestR from GitHub
install_github("akamoske/LiDARforestR")

# Load the library
library(LiDARforestR)

Now all functions should be available.

Downloading example data

NEON provides a teaching LiDAR dataset that is easy to download via R. We can use this file as a test dataset here. Code to download this .las file follows:

# Install missing R package if needed
list.of.packages <- c("uuid","rlas","devtools")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[, "Package"])]
if (length(new.packages)) {
  print("installing : ")
  print(new.packages)
  install.packages(new.packages, repos = "http://cran.rstudio.com/", dependencies = TRUE)
}

# Create a scratch folder to contain example LiDAR dataset
scratch_folder <- file.path("~/scratch/neon_data/")
if (! file.exists(scratch_folder)) dir.create(scratch_folder,recursive=TRUE)
setwd(file.path(scratch_folder))
getwd()

# Download NEON example .las file
download.file(url = "https://ndownloader.figshare.com/files/7024955",
              destfile = file.path(scratch_folder,"neon_lidar_example.las"),
              method = "auto",
              mode = "wb")

Example of usage (after installation)

Once the package is loaded into your R session, this is the an example of how to use the functions in this package
to estimate LAD and LAI:

# Convert .laz or .las files into a list of voxelized lidar arrays
laz.data <- laz.to.array(laz.files.path = "./", 
                         voxel.resolution = 10, 
                         z.resolution = 1)

# Level each voxelized array in the list to mimic a canopy height model
level.canopy <- canopy.height.levelr(lidar.array.list = laz.data)

# Estimate LAD for each voxel in leveled array in the list 
lad.estimates <- machorn.lad(leveld.lidar.array.list = level.canopy, 
                             voxel.height = 1, 
                             beer.lambert.constant = NULL)

# Convert the list of LAD arrays into a single raster stack
lad.raster <- lad.array.to.raster.stack(lad.array.list = lad.estimates, 
                                        laz.array.list = laz.data, 
                                        epsg.code = 32618)

# Create a single LAI raster from the LAD raster stack
lai.raster <- raster::calc(lad.raster, fun = sum, na.rm = TRUE)

# Generate a quick raster plot of the resulting total canopy LAI values for each pixel
plot(lai.raster)

License

This project is licensed under the GNU GPUv2 License - see the LICENSE.md file for details

[serbinsh]

@akamoske Looks good and it seems that the "wb" option works fine on linux machines.

FYI - for this new example ,I think you need to change the ESPG code to 32611 (http://www.spatialreference.org/ref/epsg/wgs-84-utm-zone-11n/) if your example is using the NEON layer from California

Or did you change the example file to a SERC file in the example?

[akamoske]

Good eyes! I didn't notice that part. Just forgot to change a couple lines.

[akamoske]

LiDARforestR

R package to estimate leaf area density (LAD) and leaf area index (LAI) from airborne LiDAR point clouds.

Information

The theory behind this R package is described in:

Kamoske A.G., Dahlin K.M., Stark S.C., and Serbin S.P. 2018. Leaf area density from airborne LiDAR: Comparing sensors and resolutions in a forest ecosystem. In preparation for submission to Forest Ecology and Management.

Corresponding Author

Code written by: Aaron G. Kamoske, PhD Student

• Michigan State University, Department of Geography
• ERSAM Lab
• [email protected]

Contributing Authors

Scott C. Stark

• Michigan State University, Department of Forestry
• Tropical Forestry Ecology Lab
• [email protected]

Shawn P. Serbin

• Brookhaven National Laboratory, Environmental and Climate Sciences Department
• [email protected]

Installation

The easiest way to install LiDARforestR is via install_github from the devtools package:

# If you haven't already installed this package and its dependencies
install.packages("devtools")

# If you alread have devtools installed or just installed it
library(devtools)

# Install LiDARforestR from GitHub
install_github("akamoske/LiDARforestR")

# Load the library
library(LiDARforestR)

Now all functions should be available.

Downloading example data

NEON provides a teaching LiDAR dataset that is easy to download via R. We can use this file as a test dataset here. Code to download this .las file follows:

# Install missing R package if needed
list.of.packages <- c("uuid","rlas","devtools")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[, "Package"])]
if (length(new.packages)) {
  print("installing : ")
  print(new.packages)
  install.packages(new.packages, repos = "http://cran.rstudio.com/", dependencies = TRUE)
}

# Create a scratch folder to contain example LiDAR dataset
scratch_folder <- file.path("~/scratch/neon_data/")
if (! file.exists(scratch_folder)) dir.create(scratch_folder,recursive=TRUE)
setwd(file.path(scratch_folder))
getwd()

# Download NEON example .las file
download.file(url = "https://ndownloader.figshare.com/files/7024955",
              destfile = file.path(scratch_folder,"neon_lidar_example.las"),
              method = "auto",
              mode = "wb")

Example of usage (after installation)

Once the package is loaded into your R session, this is the an example of how to use the functions in this package
to estimate LAD and LAI:

# Convert .laz or .las files into a list of voxelized lidar arrays
laz.data <- laz.to.array(laz.files.path = "./", 
                         voxel.resolution = 10, 
                         z.resolution = 1)

# Level each voxelized array in the list to mimic a canopy height model
level.canopy <- canopy.height.levelr(lidar.array.list = laz.data)

# Estimate LAD for each voxel in leveled array in the list 
lad.estimates <- machorn.lad(leveld.lidar.array.list = level.canopy, 
                             voxel.height = 1, 
                             beer.lambert.constant = NULL)

# Convert the list of LAD arrays into a single raster stack
lad.raster <- lad.array.to.raster.stack(lad.array.list = lad.estimates, 
                                        laz.array.list = laz.data, 
                                        epsg.code = 32611)

# Create a single LAI raster from the LAD raster stack
lai.raster <- raster::calc(lad.raster, fun = sum, na.rm = TRUE)

# Generate a quick raster plot of the resulting total canopy LAI values for each pixel
plot(lai.raster)

License

This project is licensed under the GNU GPUv2 License - see the LICENSE.md file for details

[akamoske]

Also, I have a ton more functions to add to this overtime. I need to clean them up a little bit and Kyla recently did some updates to them so that they run super fast now. Things like extracting LAD information from the LAD stack based on point locations (output is a CSV which makes it super easy to filter by height), creating square plots based on a distance from a centroid and extracting averaged LAD profiles for each plot, making DTM, DSM and CHMs with interpolation for ground points, Beer-Lambert coefficient generation from the hemi photos, and then a bunch of stuff to make dealing with NEON's naming convention. Probably will start working on all that after going to HARV in a little bit and after comps.