2

I am using the ArcGIS REST API to get lidar data in R, and while I know the image server I'm pulling from has two layers, I can't figure out how to choose the one I want. I know there are two layers because I can load them both in QGIS as WCS layers. The one I want has the ID "0" and is a single band raster with the raw elevation data, and the second (ID "1") is a 3 band RGB image of the elevations with a predefined color palette (without the elevation data itself).

When I fetch the data with the R package arcpullr without specifying anything, I get the wrong layer (with ID of "1"):

library(arcpullr)

# Creating an SF area of interest layer that will be used as bounding box
aoi <- 
  st_sf(st_as_sfc(st_bbox(c(xmin = -8071116, xmax = -8066084, ymin = 5523883, 
                            ymax = 5528915), crs = st_crs(3857))))

# Extracting a 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)

I tried adding the service layer ID to the url:

endpoint <- "https://maps.vcgi.vermont.gov/arcgis/rest/services/EGC_services/IMG_VCGI_LIDARNDSM_WM_CACHE_v1/ImageServer/0/"

but I get the error, Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘extent’ for signature ‘"NULL"’.

And I have tried adding the layer as an argument to get_image_layer() like so:

ndsm <- get_image_layer(url = endpoint, sf_object = aoi, ID = "0")

but the argument is just ignored and I still get back the wrong layer.

I can't find any information about the separate layers in the service webpage. In QGIS, the one I want, with the ID of "0" is named "x57cab711_bbd5_4ef5_b452_19ce02cd64bfy0.afr", and the other is named "x57cab711_bbd5_4ef5_b452_19ce02cd64bfy0.afr_nDSMcache_PrecipMinMaxStrch_0p5_47_v1".

What am I missing?

nealmaker
  • 157
  • 8
  • Any progress on this question? I'm getting the same error even trying to replicate the example for `get_image_layer()` from the vignette. I initially thought it might just be a crs mismatch, but I've tried using both the supplied example polygon 'wis_poly' as well as querying for the equivalent poly straight from the widnr rest endpoint using `get_spatial_layer()`. Neither approach is working for me to retrieve an image layer. – W. Kessler Feb 01 '23 at 17:53
  • No progress. The `get_image_layer()` vignette example is working for me. I think that error comes up when `arcpuller` can't connect to the endpoint you gave it. – nealmaker Feb 03 '23 at 00:16

1 Answers1

0

I doubt this is a satisfactory answer to your question, but I've been corresponding with Paul Frater, the creator of Arcpullr about basically this issue. I asked him how to pull the raw elevation data from an image service that was dishing out Lidar derived DEMs. I also noticed that the rasterbrick returned by get_image_layer() returned a 3band RGB raster and nothing else. When I asked how to get the raw raster values he said:

"That’s a good question….my initial gut reaction is because the object is being pulled from an ImageServer (hence only color representation at each pixel). I would think to get actual values you would need to be pulling a raster layer from a MapServer layer type."

So, seems like at least with arcpullr, the data needs to be a part of a published map service, not simply an image service. Doesn't seem like this should be the case though....

EDIT: I did find this resource, which outlines some of the parameters in the export image functionality of the REST API. If you know it's multiband, bandIds parameter or maybe sliceID

I also poked around in the arcpullr github, and I think the issue is, though, that the arcpullr function get_image_layer is, on the back end simply grabbing the href object from the returned JSON and passing it to the raster package function raster() as a URL. The result is just a 3band RGB png image. No underlying values are to be found. Hope this helps!

EDIT 2: Okay, making progress. I found (finally) that changing the default format in the get_image_layer(format = 'tiff') call from png to tiff, really a geotiff did the trick. granted the layer have been testing on has been a single band raster, but I imagine that this is going to do the trick for both of us.

W. Kessler
  • 355
  • 3
  • 9