Convert NHDPlus elevation derivatives into standalone mreged tif files. See comment for README.

convert.R

# Usage:
# docker run --mount type=bind,source="$(pwd)",target=/data -w /data dblodgett/hydrogeoenv-custom:latest Rscript ./convert.R

temp_source <- tempfile(fileext = ".R")

download.file("https://raw.githubusercontent.com/dblodgett-usgs/hyRefactor/master/R/download_fdr_fac.R",
              temp_source)

library("rvest")
library("xml2")
library("httr")

source(temp_source)

dir.create("nhd_data", showWarnings = FALSE)

product <- c("DEM", "hydroDEM", "FDRFAC")

out_dir <- "/data/nhd_data/"

download_elev(product, out_dir, regions = "01")

sevz <- list.files(out_dir, pattern = "*.7z", full.names = TRUE)

for(out_fi in sevz) {
  system(paste0("7z -y -o", path.expand(out_dir), " x ",
                path.expand(out_fi), " -aos"), intern = TRUE)
}

convert_to_tif <- function(aux_files, search_string, out_dir, nodata = -32768) {
  
  dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
  
  in_files <- aux_files[grepl(search_string, aux_files)]
  
  in_paths <- unique(gsub(".aux", "", in_files))
  
  out_paths <- paste0("/data/", 
                      out_dir, "/",
                      gsub("/", "_", dirname(in_paths)), 
                      '.tif')
  
  in_paths <- paste0("/data/", in_paths)
  
  run_list <- lapply(seq_len(length(in_paths)), function(x) list(in_paths[x], out_paths[x]))
  
  pbapply::pblapply(run_list, function(x) {
    
    if(!file.exists(x[[2]])) {
      command <- paste("gdalwarp -dstnodata", nodata, "-co TILED=YES -co BLOCKXSIZE=256",
                       "-co BLOCKYSIZE=256 -co COMPRESS=DEFLATE -co ZLEVEL=4", x[[1]], 
                       x[[2]])
      
      message(command)
      
      system(command)
    }
    
  })
  
  out_file <- paste0("/data/", out_dir, "_all.tif")
  
  message(paste("Merging", length(out_paths), "files.", paste(out_paths, collapse = "\n")))
  
  files <- paste(out_paths, collapse = " ")
  
  if(!file.exists(out_file)) {
  
  command <-paste('gdal_merge.py -o', out_file, 
                  '-a_nodata', nodata, '-n', nodata, '-co BIGTIFF=YES -co TILED=YES', 
                  '-co BLOCKXSIZE=256 -co BLOCKYSIZE=256 -co COMPRESS=DEFLATE', 
                  '-co ZLEVEL=4 -init', nodata, files)

  } else {
    command <- "echo 'file exists'"
  }
  
  command
  
}

all_files <- list.files("nhd_data", full.names = TRUE, 
                        recursive = TRUE, include.dirs = TRUE)

all_files <- all_files[!grepl(".*NHDPlusCI.*|.*NHDPlusHI.*|.*NHDPlusPI.*", all_files)]

# "DEM" files ending in "elev_cm.aux"
dem <- convert_to_tif(all_files,  "elev_cm.aux$|elev_cm$", "DEM")

message(dem)

system(dem)

# hydroDEM files ending in hydrodem.aux
hydrodem <- convert_to_tif(all_files, "hydrodem.aux$", "hydroDEM")

message(hydrodem)

system(hydrodem)

# fdr ending in fdr.aux
fdr <- convert_to_tif(all_files, "fdr.aux$|fdr$", "fdr", nodata = 0)

message(fdr)

system(fdr)

# fac ending in fac.aux
fac <- convert_to_tif(all_files, "fac.aux$|fac$", "fac")

message(fac)

system(fac)

Merge NHDPlusV2 Elevation Derivatives

The R script in this repository is designed to be run in the docker container: dblodgett-usgs/hydrogeoenv:custom or code.chs.usgs.gov:5001/wma/hydrogeoenv:custom-latest.

Create a working directory, place convert.R in it and run one of the following commands.

On windows, use:

docker run --mount type=bind,source="%CD%",target=/data -w /data code.chs.usgs.gov:5001/wma/hydrogeoenv:custom-latest Rscript ./convert.R

On Linux, use:

docker run --mount type=bind,source="$(pwd)",target=/data -w /data code.chs.usgs.gov:5001/wma/hydrogeoenv:custom-latest Rscript ./convert.R