get_map_layer() throwing error

[joshuajdthompson]

First off, I love this package!

I'm having trouble replicating your workflow for get_map_layer().

md_state <- get_spatial_layer("https://geodata.md.gov/imap/rest/services/Boundaries/MD_PoliticalBoundaries/FeatureServer/0")
# this works fine 

get_service_type("https://geodata.md.gov/imap/rest/services/Elevation/MD_Bathymetry/MapServer/1")
# correctly identifies object type 

md_bathyDepth  <- get_map_layer("https://geodata.md.gov/imap/rest/services/Elevation/MD_Bathymetry/MapServer/1",md_state)
# throws an error

I am getting the following error:

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘raster’ for signature ‘"NULL"’

I've been furiously Googling but have not come up with a solution to debug but no luck yet.

[pfrater]

Hi @joshuajdthompson,

Apologies for the delayed response. I'm glad you find the package helpful. Thanks for being a user.

I think the problem you're having when calling get_map_layer is that you are including the /1 on the end (i.e. end the url at MapServer). The function calls the "Export Map" supported operation on the REST API. So that's normally how this function works.

I was able to pull your raster layer once I did this.

[ozjimbob]

There is a deeper issue here - what if the layer we want to retrieve is not layer 1 - but layer 10, 20 etc.? I've been trying to figure out how to have get_map_layer retrieve anything but the first layer in a set of layers. For example this map server endpoint:

https://services.thelist.tas.gov.au/arcgis/rest/services/Public/EmergencyManagementPublic/MapServer

Includes 68 different layers - they aren't just different bands of an image, they are completely separate data sets. Querying that URL using get_map_layer() retrieves just layer 0 (eg. "Nearby Safer Places"). If I want to retrieve, instead, layer 10 ("BOM Soil Dryness Index"), there doesn't appear to be a way to do it. Adding /10 to the URL (which seems to be the obvious way to do it if you're browsing the Services Directory) gives the error mentioned in the first post. Adding a parameter id=10 to the get_map_layer function doesn't do it and layer 0 is still retrieved. Adding a parameter layer="BOM Soil Dryness Index" doesn't work either - both those parameters are just guesses based on what you see in WMS URLs, so possibly the solution is we need to add a different parameter than those?

[pfrater]

It looks like these sub-layers are a collection of FeatureLayers, so you should just be able to use the get_spatial_layer() or get_layer_by_*() functions (where * is poly, point, etc.) for each "sub-layer". i.e.

url <- "https://services.thelist.tas.gov.au/arcgis/rest/services/Public/EmergencyManagementPublic/MapServer/1"
test <- get_spatial_layer(url, head = TRUE)

[ozjimbob]

Some are, many are actually rasters, eg. example code for layer 10, the "BOM Soil Dryness Index", based on your suggestion:

library(arcpullr)
library(sf)

server <- "https://services.thelist.tas.gov.au/arcgis/rest/services/"
path <- "Public/EmergencyManagementPublic/MapServer/10"
url <- paste0(server,path)
wi_SDI <- get_spatial_layer(url, head=TRUE)

Error returned:

Error in get_spatial_layer(url, head = TRUE) : 
  return_geometry is NULL and layer geometry type 
(e.g. 'esriGeometryPolygon' or 'esriGeometryPoint' or 'esriGeometryPolyline' )

Presumably because layer 10 is a raster, not a polygon, point or polyline.

[zac-driscoll]

Hey folks,

Sorry I'm late to this party. I dug into this a bit and found that the export operation for the map services has an argument called "layers". This argument determines which layer is shown on the exported map. I was able to use this to retrieve specific layers from your Chesapeake Bay example. Could you try this and see if it works for you?

library(arcpullr)

md_state <- get_spatial_layer("https://geodata.md.gov/imap/rest/services/Boundaries/MD_PoliticalBoundaries/FeatureServer/0")

#gets layer ID 1
md_bathyDepth  <- get_map_layer("https://geodata.md.gov/imap/rest/services/Elevation/MD_Bathymetry/MapServer",md_state,
                                layers = "show:1") 

plot_layer(md_bathyDepth , plot_pkg = "base")

#gets layers with  ID 1 and2
md_bathyDepth2  <- get_map_layer("https://geodata.md.gov/imap/rest/services/Elevation/MD_Bathymetry/MapServer",md_state,
                                layers = "show:1,2") 

plot_layer(md_bathyDepth2 , plot_pkg = "base")