如何从 R 中的 ESRI REST API MapServer 下载栅格图层？

Question

我想在这个端点下载栅格图层：https://coast.noaa.gov/arcgis/rest/services/FloodExposureMapper/CFEM_CoastalFloodHazardComposite/MapServer/0，但我一生都无法弄清楚如何指定有效的查询。具体来说，我想要栅格和关联的像元值，而不仅仅是具有 RGB 值或类似值的图像。

我尝试使用

arcpullr

包的

get_map_layer()

功能：

al = tigris::states(cb = TRUE, year = 2022) %>% filter(str_detect(NAME, "Alabama")) %>% st_transform(3857)

arcpullr::get_map_layer(url, sf_object = al, bbox_crs = 3857)

以及更多使用

httr

手动构建 API 调用。一切都无济于事；我收到各种“无效 URL”或“无效或缺少输入参数”错误，这些错误并不是特别有用。

有什么建议吗？

Answer 1

似乎你必须将 url 指向 MapServer，而不是图层，例如：

url <- "https://coast.noaa.gov/arcgis/rest/services/FloodExposureMapper/CFEM_CoastalFloodHazardComposite/MapServer/"

al <- geodata::gadm("USA", path = "data") 
al <- terra::subset(al, al$NAME_1 == "Alabama") |>
  terra::project(y = "EPSG:3857")

a <- arcpullr::get_map_layer(url, sf_object = sf::st_as_sf(al))
a

class      : RasterLayer 
dimensions : 400, 400, 160000  (nrow, ncol, ncell)
resolution : 1583.956, 1583.956  (x, y)
extent     : -9966272, -9332689, 3531506, 4165088  (xmin, xmax, ymin, ymax)
crs        : +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs 
source     : memory
names      : file283b0482e977c 
values     : 1, 10  (min, max)

让我们绘制它：

terra::plot(a)

由于某些原因，栅格被翻转了，所以让我们稍微处理一下：

a <- terra::rast(a) |>
  terra::flip()
a <- terra::mask(a, mask = al) |>
  terra::crop(al)
terra::plot(a)
terra::plot(al, add = TRUE)

请注意，{httr} 设置的默认超时可能不足以获取数据（因此不是我这边的正确表示）。此外，在检索图层期间，您可能会收到如下警告：

Warning messages:
1: In get_sf_crs(sf_object) : NAs introduced by coercion
2: [rast] unknown extent 
3: [rast] unknown extent 

或强制安装 {rvest} 和/或 {magick} 软件包。