Why is arcpullr pulling an empty RasterBrick from ArcGIS REST API?-CodePudding

I am having trouble getting raster data in R from the ArcGIS REST API. It is located at the endpoint https://maps.vcgi.vermont.gov/arcgis/rest/services/EGC_service/IMG_VCGI_LIDARNDSM_WM_CACHE_v1/ImageServer and I want to fetch just the data for my area of interest. I can successfully load the data in QGIS as a web coverage service using the url "https://maps.vcgi.vermont.gov/arcgis/services/EGC_services/IMG_VCGI_LIDARNDSM_WM_CACHE_v1/ImageServer/WCSServer", but want something I can turn into a raster layer or raster brick in R.

I would prefer a high-level method and have tried the arcpullr package, using an sf polygon object of my area of interest to define the bounding box (aoi in the code below, which is a simple feature collection with 1 feature and 1 field; Bounding box: xmin: 499682.2 ymin: 208467.7 xmax: 503271.3 ymax: 212056.7; Projected CRS: NAD83 / Vermont, 32145). This returns an empty raster brick (all values are NA) in the correct location:

endpoint <- "https://maps.vcgi.vermont.gov/arcgis/rest/services/EGC_services/IMG_VCGI_LIDARNDSM_WM_CACHE_v1/ImageServer/"
ndsm <- get_image_layer(url = endpoint, sf_object = aoi)

I have also tried to construct the request myself based on the API reference, using various queries, and have only been able to write empty image files (though the server returns status codes of '200'). Here is an example:

url <- "https://maps.vcgi.vermont.gov/arcgis/rest/services/EGC_services/IMG_VCGI_LIDARNDSM_WM_CACHE_v1/ImageServer/exportImage?f=html&bbox=499682.2,208467.7,503271.3,212056.7&imageSR=32145&bboxSR=32145&format=png"
file = tempfile(fileext = ".png")
httr::GET(url = url, httr::write_disk(file))

I would appreciate any ideas for getting arcpullr to work, suggestions for other relatively easy approaches, or (as a last resort) resources for learning more about APIs.

CodePudding user response：

The problem you are experiencing is that you are not using the right crs. You simply need to convert your aoi from EPSG: 32145 to EPSG: 3857 (the spatial projection information is provided on the Service webpage).

So, please find below a little reprex.

Reprex

Code

library(arcpullr)

# Creating the aoi in `EPSG: 32145` and converting it into the right crs (i.e. `EPSG: 3857`)
aoi <- st_sf(st_as_sfc(st_bbox(c(xmin = 499682.2, xmax = 503271.3, ymin = 208467.7, ymax = 212056.7), crs = st_crs(32145)))) %>% 
  st_transform(., st_crs(3857))

# Extracting the raster corresponding to the `aoi`
endpoint <- "https://maps.vcgi.vermont.gov/arcgis/rest/services/EGC_services/IMG_VCGI_LIDARNDSM_WM_CACHE_v1/ImageServer/"
ndsm <- get_image_layer(url = endpoint, sf_object = aoi)

Output (the rasterbrick contains the pixel values)

ndsm
#> class      : RasterBrick 
#> dimensions : 400, 400, 160000, 3  (nrow, ncol, ncell, nlayers)
#> resolution : 12.58247, 12.58247  (x, y)
#> extent     : -8071116, -8066083, 5523882, 5528915  (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      : X_ags_6919e084_8cf5_46f6_b95c_1090872a66d0.1, X_ags_6919e084_8cf5_46f6_b95c_1090872a66d0.2, X_ags_6919e084_8cf5_46f6_b95c_1090872a66d0.3 
#> min values :                                            1,                                           49,                                            0 
#> max values :                                          254,                                          254,                                          149

Visualization

raster::plotRGB(ndsm)

^{Created on 2022-02-21 by the reprex package (v2.0.1)}