How can I delineate multiple watersheds in R using the streamstats package?-CodePudding

There is an R package in development that I would like to use called streamstats. What it does is delineate a watershed (within the USA) for a latitude & longitude point along a body of water and provides watershed characteristics such as drainage area and proportions of various land covers. What I would like to do is extract some watershed characteristics of interest from a data frame of several lat & long positions.

I can get the package to do what I want for one point

devtools::install_github("markwh/streamstats")
library(streamstats)

setTimeout(120)

dat1 <- data.frame(matrix(ncol = 3, nrow = 3))
x <- c("state","lat","long")
colnames(dat1) <- x
dat1$state <- c("NJ","NY","VA")
dat1$lat <- c(40.99194,42.02458,38.04235)
dat1$long <- c(-74.28000,-75.11928,-79.88144)

test_dat <- dat1[1,]

ws1 <- delineateWatershed(xlocation = test_dat$long, ylocation = test_dat$lat, crs = 4326, 
                          includeparameters = "true", includeflowtypes = "true")

chars1 <- computeChars(workspaceID = ws1$workspaceID, rcode = "MA")
chars1$parameters

However what I would like is to be able to give the delineateWatershed function several watersheds at once (i.e., all 3 locations found in dat1) and combine the chars1$parameters output variables DRNAREA,FOREST,LC11DEV, and LC11IMP into a data frame. Maybe this could be achieved with a for loop?

The ideal output would look like this

  state      lat      long DRNAREA FOREST LC11DEV LC11IMP
1    NJ 40.99194 -74.28000     160   66.2   26.20    5.50
2    NY 42.02458 -75.11928     457   89.3    2.52    0.18
3    VA 38.04235 -79.88144     158     NA    4.63    0.20

CodePudding user response：

I would put what you have in a function then use purrr::pmap_df() to loop through each row in dat1 then bind all the results together. See also this answer

library(dplyr)
library(purrr)
library(tidyr)
library(streamstats)

setTimeout(120)

dat1 <- data.frame(matrix(ncol = 3, nrow = 2))
colnames(dat1) <- c("state", "lat", "long")
dat1$state <- c("NJ", "NY")
dat1$lat <- c(40.99194, 42.02458)
dat1$long <- c(-74.28000, -75.11928)
dat1
#>   state      lat      long
#> 1    NJ 40.99194 -74.28000
#> 2    NY 42.02458 -75.11928

Define a function for catchment delineation

catchment_delineation <- function(rcode_in, lat_y, long_x) {
  
  print(paste0("Processing for lat = ", lat_y, " and long = ", long_x))

  ws <- delineateWatershed(xlocation = long_x, ylocation = lat_y, crs = 4326, 
                           includeparameters = "true", includeflowtypes = "true")
  ws_properties <- computeChars(workspaceID = ws$workspaceID, rcode = rcode_in)

  # keep only what we need
  ws_properties_df <- ws_properties$parameters %>% 
    filter(code %in% c("DRNAREA", "FOREST", "LC11DEV", "LC11IMP")) %>% 
    mutate(ID = ws$workspaceID, 
           state = rcode_in,
           long = long_x,
           lat = lat_y)

  return(ws_properties_df)

}

Apply the function to each row in dat1 data frame

catchment_df <- pmap_df(dat1, ~ catchment_delineation(..1, ..2, ..3))
#> https://streamstats.usgs.gov/streamstatsservices/watershed.geojson?rcode=NJ&xlocation=-74.28&ylocation=40.99194&includeparameters=true&includeflowtypes=true&includefeatures=true&crs=4326https://streamstats.usgs.gov/streamstatsservices/parameters.json?rcode=NJ&workspaceID=NJ20210923064141811000&includeparameters=truehttps://streamstats.usgs.gov/streamstatsservices/watershed.geojson?rcode=NY&xlocation=-75.11928&ylocation=42.02458&includeparameters=true&includeflowtypes=true&includefeatures=true&crs=4326https://streamstats.usgs.gov/streamstatsservices/parameters.json?rcode=NY&workspaceID=NY20210923064248530000&includeparameters=true

catchment_df
#>                       ID                            name
#> 1 NJ20210923064141811000                   Drainage Area
#> 2 NJ20210923064141811000                  Percent Forest
#> 3 NJ20210923064141811000 Percent Developed from NLCD2011
#> 4 NJ20210923064141811000     Percent_Impervious_NLCD2011
#> 5 NY20210923064248530000                   Drainage Area
#> 6 NY20210923064248530000                  Percent Forest
#> 7 NY20210923064248530000 Percent Developed from NLCD2011
#> 8 NY20210923064248530000     Percent_Impervious_NLCD2011
#>                                                                          description
#> 1                                            Area that drains to a point on a stream
#> 2                                               Percentage of area covered by forest
#> 3                  Percentage of developed (urban) land from NLCD 2011 classes 21-24
#> 4 Average percentage of impervious area determined from NLCD 2011 impervious dataset
#> 5                                            Area that drains to a point on a stream
#> 6                                               Percentage of area covered by forest
#> 7                  Percentage of developed (urban) land from NLCD 2011 classes 21-24
#> 8 Average percentage of impervious area determined from NLCD 2011 impervious dataset
#>      code         unit  value state      long      lat
#> 1 DRNAREA square miles 160.00    NJ -74.28000 40.99194
#> 2  FOREST      percent  66.20    NJ -74.28000 40.99194
#> 3 LC11DEV      percent  26.20    NJ -74.28000 40.99194
#> 4 LC11IMP      percent   5.50    NJ -74.28000 40.99194
#> 5 DRNAREA square miles 457.00    NY -75.11928 42.02458
#> 6  FOREST      percent  89.30    NY -75.11928 42.02458
#> 7 LC11DEV      percent   2.52    NY -75.11928 42.02458
#> 8 LC11IMP      percent   0.18    NY -75.11928 42.02458

Reshape the result to desired format

catchment_reshape <- catchment_df %>% 
  select(state, long, lat, code, value) %>% 
  pivot_wider(names_from = code,
              values_from = value)
catchment_reshape
#> # A tibble: 2 x 7
#>   state  long   lat DRNAREA FOREST LC11DEV LC11IMP
#>   <chr> <dbl> <dbl>   <dbl>  <dbl>   <dbl>   <dbl>
#> 1 NJ    -74.3  41.0     160   66.2   26.2     5.5 
#> 2 NY    -75.1  42.0     457   89.3    2.52    0.18

^{Created on 2021-09-22 by the reprex package (v2.0.1)}