Thanks for all your work on this package, it has been super helpful for pulling vector layers into R and I recently decided to give it a try for some raster layers that I'm interested in using for an analysis. However, I seem to be running into an issue with the 'get_map_layer' function.

I'm trying to pull burn severity rasters from MTBS for the Sierra Nevada region. Here's the code:

library(tidyverse)
library(sf)
library(arcpullr)

# create sf object; sierra nevada ecoregion
sierra_nevada_sf <-
  get_spatial_layer(
    url = 'https://edits.nationalmap.gov/arcgis/rest/services/landforms/landforms_epa_eco_regions/MapServer/0',
    where = "US_L3NAME = 'Sierra Nevada'"
  ) |>
  # just use a simple bounding box for now
  st_bbox() |> 
  st_as_sfc() |> 
  st_as_sf()

# query mtbs rasters
ca_fire <- get_map_layer('https://apps.fs.usda.gov/arcx/rest/services/RDW_Wildfire/MTBS_CONUS/MapServer', sierra_nevada_sf)

Which gives me the following error:

Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer",  : 
  Cannot create a RasterLayer object from this file.
In addition: Warning message:
`C:/Users/jmwin/AppData/Local/Temp/RtmpcZctRJ/file871c40e54c95' not recognized as a supported file format. (GDAL error 4) 

Initially, I was hoping to download each year of burn severity data separately. It seems this might be possible according to a previous issue, but I get the same error above:

# doesn't work either; 1984 data
ca_fire <- 
  get_map_layer(
  'https://apps.fs.usda.gov/arcx/rest/services/RDW_Wildfire/MTBS_CONUS/MapServer/', 
  sierra_nevada_sf, 
  layers = 'show:38'
  )

This led me back to the vignette, but I also get the error when trying to reproduce the Wiscland2 code:

# wi_landcover is part of the package; this doesn't work
wi_landcover2 <- get_map_layer(landcover_url, wis_poly)

Here's my session info if it is helpful:

> sessionInfo()
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22621)

Matrix products: default


locale:
[1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8    LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                           LC_TIME=English_United States.utf8    

time zone: America/Chicago
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] arcpullr_0.2.8  sf_1.0-14       lubridate_1.9.3 forcats_1.0.0   stringr_1.5.0   dplyr_1.1.3     purrr_1.0.2     readr_2.1.4    
 [9] tidyr_1.3.0     tibble_3.2.1    ggplot2_3.4.4   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] gtable_0.3.4       raster_3.6-23      htmlwidgets_1.6.2  lattice_0.21-8     RMariaDB_1.2.2     tzdb_0.4.0        
 [7] vctrs_0.6.3        tools_4.3.1        generics_0.1.3     curl_5.0.2         proxy_0.4-27       fansi_1.0.4       
[13] pkgconfig_2.0.3    KernSmooth_2.23-21 assertthat_0.2.1   lifecycle_1.0.3    compiler_4.3.1     tictoc_1.2        
[19] munsell_0.5.0      terra_1.7-46       keyring_1.3.1      codetools_0.2-19   htmltools_0.5.6    class_7.3-22      
[25] pillar_1.9.0       crayon_1.5.2       classInt_0.4-10    DT_0.29            tidyselect_1.2.0   digest_0.6.33     
[31] stringi_1.7.12     rprojroot_2.0.3    fastmap_1.1.1      grid_4.3.1         here_1.0.1         colorspace_2.1-0  
[37] cli_3.6.1          magrittr_2.0.3     utf8_1.2.3         e1071_1.7-13       withr_2.5.2        scales_1.2.1      
[43] sp_2.0-0           bit64_4.0.5        timechange_0.2.0   httr_1.4.7         bit_4.0.5          hms_1.1.3         
[49] rlang_1.1.1        Rcpp_1.0.11        glue_1.6.2         DBI_1.1.3          rstudioapi_0.15.0  vroom_1.6.3       
[55] jsonlite_1.8.7     R6_2.5.1           units_0.8-4 

Thanks!

Hi there, I'm so excited to have found this package! I'm sure it will save me from future headaches 😄

I ran into issues pulling data for my own example with PASDA landcover data, so I instead tried reproducing the example included in the documentation. But I ran into trouble there too. Here are the steps I followed:

# WDNR Server
image_server <- "https://dnrmaps.wi.gov/arcgis_image/rest/services/"

# WI Landcover Type URL
landcover_path <- "DW_Land_Cover/EN_Land_Cover2_Lev2/MapServer"
landcover_url <- paste0(image_server, landcover_path)

# defining the WI shape
wis_poly <- arcpullr::wis_poly

# pulling map data
wi_landcover <- get_map_layer(landcover_url, wis_poly)

And then I got the following error:

Warning: [rast] unknown extent
Warning: [rast] unknown extent

Any idea what might be going on here? Let me know if I can provide any addiitonal context.

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.

Hi there, I keep getting this error when using the package:

Error in .rasterObjectFromFile(x, objecttype = "RasterBrick", ...) :
Cannot create a RasterLayer object from this file.

I'm trying to pull a portion of this raster using a bounding box: https://server6.tplgis.org/arcgis6/rest/services/Heat_Severity_2021/ImageServer/query?outFields=*&where=1%3D1

Any help would be greatly appreciated!

library(arcpullr)
library(sf)

Define the URL for the ArcGIS REST API

heat_severity_url <- "https://server6.tplgis.org/arcgis6/rest/services/Heat_Severity_2021/ImageServer/"

Define the sf object for the bounding box

Define the sf object with the original coordinates

sf_object <- st_sf(st_as_sfc(st_bbox(c(xmin = -123.156738, xmax = -122.036133, ymin = 45.201393,
ymax = 45.744527), crs = st_crs(3857))))

Use the get_image_layer function to pull the raster layer from the ArcGIS REST API

heat_severity_layer <- get_image_layer(heat_severity_url, sf_object)

Hello, thank you for this helpful package package!

Im currently working on extracting elevation data (from https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer) based on the bounding box of a multiline sf object. When I input the multiline sf object, the tool does as expected and extracts the elevation data for each line feature. While this is useful, want to use a bounding box to get the elevation profile for the entire area surrounding our multiline feature to use as a basemap in tmap which seems possible given the option to use a bbox as input for this function.

However, in lines 374-378 it looks like the function doesn't accept a bbox as an input object. Whenever I run the code with a bbox input in place of an sf object, I get the following error "Error in body_config(body, match.arg(encode)) : object 'bbox_coords' not found". When looking at the code behind the function (lines 334:447 in arcpullr/R/retrieve_layers), I see that the function creates a bbox from the input sf object and then passes that into bbox_coords, but I dont see where the function passes the bbox input to bbox coords. Am I misunderstanding how to pass in the bbox as in input or is the get_image_layer unable to accept bbox as an input in place of an sf?

Thanks for the assistance!

Thanks for writing this very handy package. I think though I run into a problem with retrieving spatial data due to queries becoming too large if all records are requested in one query.

Looking at the code, I see that features are retrieved in chunks:

get_esri_features <- function(query_url, fields, where, token='', head, ...) {
  ids <- get_object_ids(query_url, where, token, ...)
  if(is.null(ids)){
    warning("No records match the search critera")
    return()
  }
  if (isTRUE(head)) {
    id_splits <- ids[1:5]
  } else if (head > 0 & head < 500) {
    id_splits <- ids[1:head]
  } else {
    id_splits <- split(ids, ceiling(seq_along(ids)/500))
  }
  results <- lapply(
    id_splits,
    get_esri_features_by_id,
    query_url,
    fields,
    token,
    ...
  )
  merged <- unlist(results, recursive=FALSE)
  return(merged)
}

I think that the default chunk size of 500 is too big for some of the things I'd like to retrieve. I tested it with smaller values, and it seemed to work. Note that you don't get any errors or warnings, just an empty result. I wouldn't be suprised if there were some server setting that limits the size of the returned query, but I'm not really familiar with the ArcGIS server.

Here is my test:

> subbioregions.sf <- get_spatial_layer("https://spatial-gis.information.qld.gov.au/arcgis/rest/services/Boundaries/AdminBoundariesFramework/MapServer/3",
+                                       idsplit=500)
> subbioregions.sf
Simple feature collection with 0 features and 0 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  WGS 84
[1] geoms
<0 rows> (or 0-length row.names)
> subbioregions.sf <- get_spatial_layer("https://spatial-gis.information.qld.gov.au/arcgis/rest/services/Boundaries/AdminBoundariesFramework/MapServer/3",
+                                       idsplit=5)
> subbioregions.sf
Simple feature collection with 133 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 137.9946 ymin: -29.1779 xmax: 153.5529 ymax: -9.067217
Geodetic CRS:  WGS 84
First 10 features:
   q_reg    q_reg_name q_sub            q_sub_name objectid                          geoms
1    BRB Brigalow Belt  11.1     Townsville Plains        1 MULTIPOLYGON (((146.615 -19...
2    BRB Brigalow Belt 11.10          Basalt Downs        2 MULTIPOLYGON (((147.6078 -2...
3    BRB Brigalow Belt 11.11   Isaac - Comet Downs        3 MULTIPOLYGON (((148.6707 -2...
4    BRB Brigalow Belt 11.12 Nebo - Connors Ranges        4 MULTIPOLYGON (((149.0068 -2...
5    BRB Brigalow Belt 11.13  South Drummond Basin        5 MULTIPOLYGON (((146.7731 -2...
6    BRB Brigalow Belt 11.14    Marlborough Plains        6 MULTIPOLYGON (((149.464 -22...
7    BRB Brigalow Belt 11.15    Claude River Downs        7 MULTIPOLYGON (((146.2638 -2...
8    BRB Brigalow Belt 11.16            Woorabinda        8 MULTIPOLYGON (((149.0985 -2...
9    BRB Brigalow Belt 11.17          Boomer Range        9 MULTIPOLYGON (((149.6715 -2...
10   BRB Brigalow Belt 11.18   Mount Morgan Ranges       10 MULTIPOLYGON (((149.8721 -2...

I added the argument idsplit in get_esri_features which then gets carried through by your ... mechanism.

Do you think this could be added to the code? Happy to make PR if you'd like.

Im having trouble this week getting a RasterBrick from the get_image_layer function. Last week this function was working but it is currently throwing an error. Heres a reproducible example using package data and the national map elevation image server.

`wis_counties<-arcpullr::wis_counties
adams<-wis_counties[wis_counties$county=="adams",]
adams_bbox<-sf::st_bbox(adams)
url<-"https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer/"

adams_elevation<-get_image_layer(url,adams)
adams_bbox_elevation<-get_image_layer(url,bbox=adams_bbox, bbox_crs="4326")`

Email of April 1:

arcpullr has strong dependence on at least one of rgdal, rgeos or
maptools, but does not seem to use any functionality in code. The usage
may have been in raster, which now uses terra instead, or may be in
examples or vignettes. Please move all use of rgdal, rgeos or maptools
to Suggests: and protect any use against these packages not being present.
The packages will be archived in October 2023. See
https://r-spatial.org/r/2022/04/12/evolution.html,
https://r-spatial.org/r/2022/12/14/evolution2.html and
https://rsbivand.github.io/csds_jan23/bivand_csds_ssg_230117.pdf and perhaps
view
https://www.youtube.com/watch?v=TlpjIqTPMCA&list=PLzREt6r1NenmWEidssmLm-VO_Y
mAh4pq9&index=1.

Best wishes,

On Mon, 19 Dec 2022, Roger Bivand wrote:

Dear Maintainer,

arcpullr depends on (depends, imports or suggests) raster and one or
more of the retiring packages rgdal, rgeos or maptools
(https://r-spatial.org/r/2022/04/12/evolution.html,
https://r-spatial.org/r/2022/12/14/evolution2.html). Since raster
3.6.3, all use of external FOSS library functionality has been
transferred to terra, making the retiring packages very likely
redundant. It would help greatly if you could remove dependencies on the
retiring packages as soon as possible.

Best wishes,

See also https://r-spatial.org/r/2023/04/10/evolution3.html.

Hello, and thank you for an incredibly useful package! I have run into an issue where a specific dataset that I am working with has one feature which somehow does not have geometry associated with it, leading to a crash in esri2sfPoint() here; presumably, if it were a different type of geometry (i.e., not a point, but perhaps a polygon) the same problem could occur. Obviously, this is not an ideal situation to be in from a data standpoint, and until now I had never actually seen a feature with no geometry- but given that it is (apparently) possible to build such a dataset in ArcGIS, it seems like something that arcpullr ought to support.

The feature in question is Object ID 6926 in this dataset: https://maps.ihs.gov/server/rest/services/OPHS/enterprise/MapServer/7/

I am not sure what the best way to handle this sort of situation would be from the standpoint of arcpullr's API and behavior. I don't know whether sf objects are able to include features that have no associated geometry, but if they are, I suppose one option would be to include the feature and its attributes? Another option would be to just skip these features altogether. In either event, I would imagine a warning to the user would be in order. Thoughts?

Ask: please create reproducible examples in help documentation and vignettes

Many of your examples in your package are not reproducible. They reference objects which do not exist or are not exported by your package.

For example the example in get_map_layer() is

wi_landcover<- get_map_layer(wi_landcover_url, wis_poly)
plot_layer(wi_landcover, outline_poly = wis_poly)

wi_landcover_url is an internal object that can be accessed like arcpullr:::wi_landcover_url. wis_poly, however, is not even an internal object.