https://github.com/cran/FedData
Tip revision: 177f9ef89b4de323fc857719df80ce4409273fd5 authored by R. Kyle Bocinsky on 22 April 2019, 12:40:03 UTC
version 2.5.7
version 2.5.7
Tip revision: 177f9ef
NHD_FUNCTIONS.R
#' Download and crop the National Hydrography Dataset.
#'
#' \code{get_nhd} returns a list of Spatial* objects extracted
#' from the National Hydrography Dataset.
#'
#' @param template A Raster* or Spatial* object to serve
#' as a template for cropping.
#' @param label A character string naming the study area.
#' @param raw.dir A character string indicating where raw downloaded files should be put.
#' The directory will be created if missing. Defaults to './RAW/NHD/'.
#' @param extraction.dir A character string indicating where the extracted and cropped NHD shapefiles should be put.
#' The directory will be created if missing. Defaults to './EXTRACTIONS/NHD/'.
#' @param force.redo If an extraction for this template and label already exists, should a new one be created?
#' @return A list of Spatial* objects extracted from the National Hydrography Dataset.
#' @export
#' @examples
#' \dontrun{
#' # Extract data for the Village Ecodynamics Project 'VEPIIN' study area:
#' # http://village.anth.wsu.edu
#' vepPolygon <- polygon_from_extent(raster::extent(672800,740000,4102000,4170000),
#' proj4string='+proj=utm +datum=NAD83 +zone=12')
#'
#' # Get the NHD (USA ONLY)
#' NHD <- get_nhd(template=vepPolygon, label='VEPIIN')
#'
#' # Plot the VEP polygon
#' plot(vepPolygon)
#'
#' # Plot the NHD data
#' plot(NHD$NHDFlowline, add=T)
#' plot(NHD$NHDLine, add=T)
#' plot(NHD$NHDArea, col='black', add=T)
#' plot(NHD$NHDWaterbody, col='black', add=T)
#' }
get_nhd <- function(template,
label,
raw.dir = "./RAW/NHD",
extraction.dir = paste0("./EXTRACTIONS/", label, "/NHD"),
force.redo = FALSE) {
raw.dir <- normalizePath(paste0(raw.dir,"/."), mustWork = FALSE)
extraction.dir <- normalizePath(paste0(extraction.dir,"/."), mustWork = FALSE)
dir.create(raw.dir, showWarnings = FALSE, recursive = TRUE)
dir.create(extraction.dir, showWarnings = FALSE, recursive = TRUE)
if (!force.redo & length(list.files(extraction.dir)) > 0) {
files <- list.files(extraction.dir)
files <- files[grepl("shp", files)]
files <- files[!grepl("template", files)]
files <- files[grepl(label, files)]
files <- gsub(".shp", "", files)
files <- files[order(files)]
shapes <- lapply(files, function(file) {
rgdal::readOGR(extraction.dir, file, verbose = F)
})
names(shapes) <- gsub(paste0(label, "_NHD"), "", files)
return(shapes)
}
if (class(template) %in% c("RasterLayer", "RasterStack", "RasterBrick")) {
template <- spdf_from_polygon(sp::spTransform(polygon_from_extent(template), sp::CRS("+proj=longlat +ellps=GRS80")))
}
message("(Down)Loading the NHD HUC4 dataset.")
HUC4 <- get_huc4(template = template, raw.dir = raw.dir)
area.list <- HUC4$HUC4 %>% as.character()
# Get the spatial data for each area
message("(Down)Loading the NHD subregion data.")
subregionShapes <- lapply(area.list, function(area) {
return(get_nhd_subregion(template = template, area = area, raw.dir = raw.dir))
})
# Get all layer names
layers <- unique(unlist(lapply(subregionShapes, names)))
# Merge like datasets
message("Merging all NHD data in study area.")
allShapes <- lapply(layers, function(layer) {
shapes <- sapply(subregionShapes, "[[", layer)
null.shapes <- sapply(shapes, is.null)
shapes <- do.call("rbind", shapes[!null.shapes])
if (is.null(shapes))
return(shapes)
shapes <- raster::crop(shapes, sp::spTransform(template, sp::CRS(raster::projection(shapes))))
if (is.null(shapes))
return(shapes)
# shapes <- spTransform(shapes,CRS(projection(template)))
layer <- gsub("NHD", "", layer)
suppressWarnings(rgdal::writeOGR(shapes,
dsn = extraction.dir,
layer = paste0(label, "_NHD_", layer),
driver = "ESRI Shapefile",
overwrite_layer = TRUE))
return(shapes)
})
names(allShapes) <- gsub("NHD", "", layers)
# Remove null layers
allShapes <- allShapes[!sapply(allShapes, is.null)]
return(allShapes)
}
#' Download a zipped directory containing a shapefile of the HUC4 subregions of the NHD.
#'
#' @param raw.dir A character string indicating where raw downloaded files should be put.
#' @return A character string representing the full local path of the HUC4 zipped directory.
#' @export
#' @keywords internal
download_huc4 <- function(raw.dir) {
# url <- 'ftp://rockyftp.cr.usgs.gov/vdelivery/Datasets/Staged/WBD/FileGDB101/WBD_National.zip'
# url <- "ftp://ftp.igsb.uiowa.edu/gis_library/USA/huc_04.zip"
url <- "https://github.com/ropensci/FedData/raw/master/data-raw/nhd_huc4.gpkg.zip"
tryCatch(download_data(url = url, destdir = raw.dir), error = function(e){message("HUC4 download not available. Using local version.")})
return(normalizePath(paste(raw.dir, "/nhd_huc4.gpkg.zip", sep = "")))
# return(normalizePath(paste(raw.dir,'WBD_National.zip',sep='')))
}
#' Download and crop a shapefile of the HUC4
#' regions of the National Hydrography Dataset.
#'
#' \code{get_huc4} returns a \code{SpatialPolygonsDataFrame} of the HUC4 regions within
#' the specified \code{template}. If template is not provided, returns the entire HUC4 dataset.
#'
#' @param template A Raster* or Spatial* object to serve
#' as a template for cropping.
#' @param raw.dir A character string indicating where raw downloaded files should be put.
#' The directory will be created if missing.
#' @return A \code{SpatialPolygonsDataFrame} of the HUC4 regions within
#' the specified \code{template}.
#' @export
#' @keywords internal
get_huc4 <- function(template = NULL, raw.dir) {
tmpdir <- tempfile()
if (!dir.create(tmpdir))
stop("failed to create my temporary directory")
huc4File <- download_huc4(raw.dir)
utils::unzip(huc4File, exdir = tmpdir)
HUC4 <- sf::st_read(stringr::str_c(tmpdir,"/data-raw/nhd_huc4.gpkg"),
quiet = TRUE) %>%
as("Spatial")
# HUC4 <- rgdal::readOGR(tmpdir, layer = "huc_04", verbose = FALSE)
#
# HUC4@proj4string <- sp::CRS("+proj=utm +zone=15 +datum=NAD83 +ellps=WGS84")
# Get a list of NHD subregions within the project study area
if (!is.null(template)) {
if (class(template) %in% c("RasterLayer", "RasterStack", "RasterBrick")) {
template <- spdf_from_polygon(sp::spTransform(polygon_from_extent(template), sp::CRS("+proj=longlat +ellps=GRS80")))
}
HUC4 <- raster::crop(HUC4, sp::spTransform(template, sp::CRS(raster::projection(HUC4))))
}
unlink(tmpdir, recursive = TRUE)
return(HUC4)
}
#' Download a zipped NHD HUC4 subregion.
#'
#' HUC4 subregion are specified by a unique character string, best obtained using the \code{\link{get_huc4}} function.
#' \code{download_nhd_subregion} returns the path to the downloaded zip file.
#'
#' @param area A 4-character string indicating the HUC4 NHD subregion to download.
#' @param raw.dir A character string indicating where raw downloaded files should be put.
#' The directory will be created if missing.
#' @return A character string representing the full local path of the downloaded zip file.
#' @export
#' @keywords internal
download_nhd_subregion <- function(area, raw.dir) {
url <- paste0("https://prd-tnm.s3.amazonaws.com/StagedProducts/Hydrography/NHD/HU4/HighResolution/GDB/NHD_H_",area,"_HU4_GDB.zip")
destdir <- raw.dir
download_data(url = url, destdir = destdir)
return(normalizePath(paste0(destdir, "/", basename(url))))
}
#' Download and crop data from a zipped HUC4 subregion
#' of the National Hydrography Dataset.
#'
#' \code{get_nhd_subregion} returns a list of \code{SpatialPolygonsDataFrame}s of the layers of the HUC4 subregion,
#' within the specified \code{template}. If template is not provided, returns the entire HUC4 subregion.
#'
#' @param template A Raster* or Spatial* object to serve
#' as a template for cropping.
#' @param area A 4-character string indicating the HUC4 NHD subregion to download and crop.
#' @param raw.dir A character string indicating where raw downloaded files should be put.
#' The directory will be created if missing.
#' @return A \code{SpatialPolygonsDataFrame} of the HUC4 regions within
#' the specified \code{template}.
#' @export
#' @keywords internal
get_nhd_subregion <- function(template = NULL, area, raw.dir) {
tmpdir <- tempfile()
if (!dir.create(tmpdir))
stop("failed to create my temporary directory")
file <- download_nhd_subregion(area = area, raw.dir = raw.dir)
utils::unzip(file, exdir = tmpdir)
# Get the path to the geodatabase
dsn <- list.files(tmpdir, full.names = T, pattern = "*.gdb")
# dsn <- dsn[grepl(".gdb", dsn)]
dsn <- normalizePath(dsn)
# List all layers in the geodatabase
layers <- rgdal::ogrListLayers(dsn)
layers <- layers[grepl("NHD", layers)]
# Get each layer in the geodatabase
shapes <- lapply(layers, function(layer) {
tryCatch(suppressWarnings(rgdal::readOGR(dsn = dsn, layer = layer, verbose = F)), error = function(e) NULL)
})
names(shapes) <- layers
# Rename the features to prepare for merging
shapes <- lapply(shapes, function(shape) {
tryCatch(sp::spChFIDs(shape, paste(area, "_", shape$Permanent_Identifier, sep = "")), error = function(e) NULL)
})
# Crop each feature to the template area
if (!is.null(template)) {
if (class(template) %in% c("RasterLayer", "RasterStack", "RasterBrick")) {
template <- spdf_from_polygon(sp::spTransform(polygon_from_extent(template), sp::CRS("+proj=longlat +ellps=GRS80")))
}
shapes <- lapply(shapes, function(shape) {
tryCatch(raster::crop(shape, sp::spTransform(template, sp::CRS(raster::projection(shape)))), error = function(e) NULL)
})
}
shapes <- shapes[!sapply(shapes, is.null)]
unlink(tmpdir, recursive = TRUE)
return(shapes)
}