In bmcgaughey1/LidarIndexR: Provides functions to build index files for lidar projects

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)

LidarIndexR

The goal of LidarIndexR is to provide functions to build index files for lidar projects. There are functions specific to local file systems and functions for remote servers using FTP and HTTPS protocols. USGS changed the format of the XML generated for directory listing on the rockyweb server (or changed the server software) sometime between 2022 and 2024. Code was modified in June 2024 to correctly parse the new output. However, I can't guarantee that the directory functions work with all HTTP(S) servers due to differences in the HTML generated to show directory contents. The changes in 2024 should make the functions more reliable given different formats but I haven't tested functions using other servers.

Installation

LidarIndexR is NOT available on CRAN.

You can install the development version from GitHub. If you don't already have the devtools package installed, uncomment the first line.

# install.packages("devtools")
devtools::install_github("bmcgaughey1/LidarIndexR")

Example #1: Create a tile index from files on a remote ftp server

This example creates a tile index for a single lidar project stored in the USGS rockyweb server. In this case, you have to know the details for the project including the folder name for the project and the projection information for the data.

library(LidarIndexR)
library(sf)
library(ggplot2)
library(viridis)

## Basic example to build an index
URL <- "https://rockyweb.usgs.gov/vdelivery/Datasets/Staged/Elevation/LPC/Projects/WA_ElwhaRiver_2015"
pointFolder <- "WA_Elwha_TB_2015/LAZ"
outputFile <- "~/WA_ElwhaRiver_2015.gpkg"

# projection info...for most projects, this can be pulled from the WESM project index
pointCRS <- 26910

# create index
BuildIndexFromPoints(URL, pointFolder, outputFile, projString = pointCRS,
                     appendInfo = data.frame("Project" = "WA_ElwhaRiver_2015"))

# read the index and display
index <- st_read(outputFile)

ggplot(index) +
  ggtitle("WA_ElwhaRiver_2015 lidar tiles") +
  theme(plot.title = element_text(hjust = 0.5)) +
  geom_sf(aes(fill = PointCount), show.legend = TRUE) +
  scale_fill_viridis()

Example #2: Create a tile index using a project from the USGS WESM project index

This example creates an index for the same lidar project as Example #1. However, this example uses the USGS project index to populate the new tile index with information regarding the acquisition and to provide projection information for the project.

You will need to change the line that reads the index to reflect the location of the WESM index on your computer.

The example requires a local copy of the USGS project index from here. The commented code in the example downloads the WESM.gpkg index.

Note: The call to \code{BuildIndexFromUSGSProjectIndexItem} can take quite a bit of time to execute. It touches every point file on the server to read bounding box information for the tiles. For projects with thousands of tiles, it can take hours to build the index.

library(LidarIndexR)
library(sf)
library(utils)
library(ggplot2)
library(viridis)

# download the USGS WESM index
# utils::download.file("https://rockyweb.usgs.gov/vdelivery/Datasets/Staged/Elevation/metadata/WESM.gpkg",
#                      "WESM.gpkg", 
#                      method = "libcurl", 
#                      mode = "wb")

# read the FESM index from a local file
index <- st_read("G:/R_Stuff/EntwineIndex/WESM.gpkg")
pointFolder <- "laz"
outputFile <- "~/WA_Elwha_TB_2015.gpkg"

item <- index[which(index$workunit == "WA_Elwha_TB_2015"), ]

# create index
BuildIndexFromUSGSProjectIndexItem(item, pointFolder, outputFile)

# read the index and display
tindex <- st_read(outputFile)

ggplot(tindex) +
  ggtitle("WA_Elwha_TB_2015 lidar tiles") +
  theme(plot.title = element_text(hjust = 0.5)) +
  geom_sf(aes(fill = PointCount), show.legend = TRUE) +
  scale_fill_viridis()

Example #3: Create an overall polygon for a project from the tile index

This example builds on Example #2 to create a polygon showing the area covered by the point tiles. The default size for the raster used to merge polygons is 512 by 512. By increasing this to 2048 by 2048, the project polygon captures the small tiles that are otherwise omitted when using the default size.

# build on the last example to generate a project polygon

tproject <- BuildProjectPolygonFromIndex(
  outputFile,
  item$workunit,
  appendInfo = item[, c("workunit_id", "collect_start", "collect_end", "p_method", "horiz_crs")],
  quiet = FALSE,
  nx = 2048, ny = 2048
)

# display the project polygon
ggplot(tproject) +
  ggtitle(paste0(item$workunit, " project area")) +
  theme(plot.title = element_text(hjust = 0.5)) +
  geom_sf(aes(fill = PID), show.legend = FALSE)

# You will need to change the line that reads the index to reflect the location of the WESM index
# on your computer.

library(LidarIndexR)
library(sf)

# test code...not to be run as part of README

index <- st_read("G:/R_Stuff/PlotClipping/WESM.gpkg")
outputFile <- "~/AK_BrooksCamp_2012.gpkg"
appendCols <- c("opr_category", "lpc_pub_date")

item <- index[which(index$workunit == "AK_BROOKSCAMP_2012"), ]

# read the index
tindex <- st_read(outputFile)

info <- data.frame(item[, appendCols])

# drop geometry
info$geometry <- NULL

# duplicate rows
info <- do.call("rbind", replicate(nrow(tindex), info, simplify = FALSE))

# cbind to tile data frame
tindex <- cbind(tindex, info)

bmcgaughey1/LidarIndexR documentation built on April 14, 2025, 6:07 p.m.