R – Downloading Areas Designated by Bounding Box of Microsoft’s USBuildingFootprints

buildingr

I would like to plot some houses in a relatively small area of MD, the bounding box is: xmin: -76.12, ymin: 38.05, xmax: -75.68, ymax: 38.42

Microsoft's USBuildingFootprints dataset can be downloaded by state here https://github.com/Microsoft/USBuildingFootprints, but is there another way to download just a requested area, in order to get around the large files and slow download speeds?

I would prefer a solution using R or an internet browser.

Best Answer

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...
Related Question