To do this in R, construct a URL string like this with your bounding box encoded in it. Use paste0 or some other string concatenation function to build it. This is based on @user30184 answer:
url="https://services.arcgis.com/P3ePLMYs2RVChkJx/ArcGIS/rest/services/MSBFP2/FeatureServer/0/query?f=json&returnGeometry=true&spatialRel=esriSpatialRelIntersects&geometry={%22xmin%22:-76.12,%22ymin%22:38.05,%22xmax%22:-75.68,%22ymax%22:38.42,%22spatialReference%22:{%22wkid%22:4326}}&geometryType=esriGeometryEnvelope&outSR=4326"
Then you can use st_read from the sf package:
bdgs = st_read(url)
Getting about 11,000 buildings:
dim(bdgs)
## [1] 11978 2
Your sf needs to be linked with a GDAL/OGR that supports these servers:
> st_drivers()["ESRIJSON",]
name long_name write copy is_raster is_vector vsi
ESRIJSON ESRIJSON ESRIJSON FALSE FALSE FALSE TRUE TRUE
If that returns NA values then your GDAL/OGR needs upgrading/rebuilding.
Here's an R function that downloads the buildings given some bounding coordinates in WGS84::
get_buildings <- function(xmin, ymin, xmax, ymax){
URL = paste0("https://services.arcgis.com/P3ePLMYs2RVChkJx/",
"ArcGIS/rest/services/MSBFP2/FeatureServer/0/",
"query?f=json&returnGeometry=true&",
"spatialRel=esriSpatialRelIntersects&",
"geometry={",
"%22xmin%22:",xmin,
",%22ymin%22:",ymin,
",%22xmax%22:",xmax,
",%22ymax%22:",ymax,
",%22spatialReference%22:{%22wkid%22:4326}}&",
"geometryType=esriGeometryEnvelope&outSR=4326")
st_read(URL, quiet=TRUE)
}
giving:
> bb = get_buildings(-76.12, 38.05, -75.68, 38.42)
> bb
Simple feature collection with 11978 features and 1 field
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -76.1202 ymin: 38.04976 xmax: -75.67844 ymax: 38.42065
Geodetic CRS: WGS 84
First 10 features:
BlockgroupID geometry
1 240399303001 POLYGON ((-75.74032 38.0858...
2 240399303002 POLYGON ((-75.7885 38.1083,...
3 240199709002 POLYGON ((-76.11691 38.3538...
4 240199709002 POLYGON ((-76.11554 38.3529...
5 240450108004 POLYGON ((-75.91377 38.2617...
6 240450108004 POLYGON ((-75.87707 38.3167...
Best Answer
Use "x <- sf::st_read" it's much faster than rgdal and will happily read much larger data sets.
When you plot it, be sure to use "plot (st_geometry(x))" first, since it will otherwise do a lot of repeated drawing for the default multi panel layout. It defaults to faceting on every column, up to some max.
Indexing with "[" works the same as with sp, so use that to plot just a few rows before trying the entire thing if it's really massive.