In bmcgaughey1/USGSlidar: Provides Functions to Access USGS LIDAR Data Collection

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Overall Goal

The USGSlidar package was developed to make it somewhat easy to query the collection of lidar data available from the U.S. Geological Survey's 3D Elevation Program (and other sources) to find out if data are available for specific locations of interest. The package accesses some of the project and tile index files available from 3DEP and uses them to search for data. The package provides links to lidar data files that can be downloaded for free directly from 3DEP servers.

The package also facilitates access to the USGS lidar collection maintained in Entwine format. This collection contains a subset of the data available from USGS but it is constantly being updated to match the data available in the national map. The index for the Entwine collection does not contain detailed project information so there is extra effort required to merge the Entwine index with information in the USGS 3DEP collection.

Packages

Load the packages needed for the examples.

library(ggplot2)
library(ggspatial)
library(USGSlidar)

Example 1 -- Query USGS Entwine lidar collection

Using the USGS Entwine data collection, query for coverage for a buffered point location. The Entwine data collection is in the web mercator projection (EPSG:3857). The point coordinate is also in the web mercator projection. The point is buffered using a circular buffer with a radius of 500m.

The Entwine data collection can be accessed using PDAL to retrieve data for specific areas of interest without needing to know the structure of the data (tile arrangement or tile naming). Unfortunately, the Entwine index file contains very few details for each acquisition so it is hard to find data for specific years. Other USGS data collections contain this information. I have developed R code that matches index files from the data hosted on the rockyweb server with the projects in the Entwine collection to produce a new index file that can be used with the USGSlidar package. This index is hosted here. The second example uses this enhanced index.

The example point is located within the 2014 Glacier Peak lidar acquisition in Washington state, USA.

For the Entwine data collection, all data are from lidar so the lidarOnly option in queryUSGSProjectIndex is not needed. In addition, this option relies on a specific field in the USGS WESM index so it should only be used with indexes that include this field.

This example will download the Entwine index produced by Howard Butler into the current working directory.

# retrieve the entwine boundary index and store in current directory
fetchUSGSProjectIndex(type = "ENTWINE")

# query to find the project polygon(s) that covers the buffered point
# the return value is the project area polygon(s) that intersect buffered point
project_polys <- queryUSGSProjectIndex(-13499097
                                       , 6139669
                                       , buffer = 500
                                       , shape = "circle"
                                       , crs = 3857
                                       , lidarOnly = FALSE
)

# now query to get the sample area attributed with the project polygon information
# this is the same call as above but with return="aoi"
point_poly <- queryUSGSProjectIndex(-13499097
                                    , 6139669
                                    , buffer = 500
                                    , shape = "circle"
                                    , return = "aoi"
                                    , crs = 3857
                                    , lidarOnly = FALSE
)

# finally, query using a point in UTM 10N within the same lidar acquisition to attribute
# the point with lidar project information. return point will be in web mercator projection.
point <- queryUSGSProjectIndex(640377.7
                                , 5330460
                                , buffer = 0
                                , shape = "circle"
                                , return = "aoi"
                                , crs = "EPSG:26910"
                                , returncrs = 3857
                                , lidarOnly = FALSE
)

# display the project area polygon and the buffered point
ggplot() +
  ggtitle("2014 Glacier Peak lidar acquisition") +
  annotation_map_tile(type = "osm") +
  layer_spatial(project_polys, fill = "cornsilk", alpha = 0.5) +
  layer_spatial(point, col = "green", size = 2) +
  annotation_spatial(point_poly, fill = "red")

cat("Attributes for the acquistion in the ENTWINE index:\n",
    paste(paste(colnames(project_polys), collapse= "\n "), "\n"))

Example 2 -- Query USGS Entwine lidar collection using enhanced ENTWINE index

This example is identical to Example 1 except it uses the enhanced index for ENTWINE data that includes information about the acquisition including the start and end dates for data collection. Using this index, you can search for data covering specific dates for a given location.

This example will download the Enhanced Entwine index produced by Robert J. McGaughey into the current working directory.

Be sure to look at the output below the map to see the additional columns of information available in the enhanced index.

# retrieve the enhanced entwine boundary index and store in current directory
fetchUSGSProjectIndex(type = "ENTWINEPLUS")

# query to find the project polygon(s) that covers the buffered point
# the return value is the project area polygon(s) that intersect buffered point
project_polys <- queryUSGSProjectIndex(-13499097
                                       , 6139669
                                       , buffer = 500
                                       , shape = "circle"
                                       , crs = 3857
)

# now query to get the sample area attributed with the project polygon information
# this is the same call as above but with return="aoi"
point_poly <- queryUSGSProjectIndex(-13499097
                                    , 6139669
                                    , buffer = 500
                                    , shape = "circle"
                                    , return = "aoi"
                                    , crs = 3857
)

# finally, query using a point in UTM 10N within the same lidar acquisition to attribute
# the point with lidar project information. return point will be in web mercator projection.
point <- queryUSGSProjectIndex(640377.7
                                , 5330460
                                , buffer = 0
                                , shape = "circle"
                                , return = "aoi"
                                , crs = 26910
                                , returncrs = 3857
)

# display the project area polygon and the buffered point
ggplot() +
  ggtitle("2014 Glacier Peak lidar acquisition") +
  annotation_map_tile(type = "osm") +
  layer_spatial(project_polys, fill = "cornsilk", alpha = 0.5) +
  layer_spatial(point, col = "green", size = 2) +
  annotation_spatial(point_poly, fill = "red")

cat("Attributes for the acquistion in the enhanced ENTWINE index:\n",
    paste(paste(colnames(project_polys), collapse= "\n "), "\n"))