This code will download ground-level LIDAR digital elevation model, the best available, for any part of state when you supply it with a shapefle. It queries the rest server with the dem indexes, and then downloads the relevant files, and joins them into one image.
library(tidyverse)
library(sf)
library(raster)
library(arcpullr)
shape_to_download <- read_sf('/path/to/shapefile_of_area_to_download.shp')
lidar.url <- tibble()
for (i in seq(2,14)) {
if (nrow(lidar.url) == 0)
lidar.url <- get_layer_by_poly(str_c('https://elevation.its.ny.gov/arcgis/rest/services/Dem_Indexes/MapServer/',i),
shape_to_download, sp_rel='intersects'
)
}
# download and merge all lidar titles that intersect
map(lidar.url$DIRECT_DL, \(x) {
lidar.file <- tempfile(fileext = '.tif')
download.file(x, destfile = lidar.file)
lidar <- rast(lidar.file)
}) %>% sprc %>% merge %>% writeRaster('output.tif')
