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.
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...
As @Spacedman wrote in the comment, building footprint data seems to be available from the ESRI service.
The documentation about the REST API is at https://developers.arcgis.com/rest/services-reference/enterprise/query-feature-service-.htm. Basic command for getting data from a BBOX that is expressed as EPSG:4326 coordinates, and output also in EPSG:4326, is like this:
You can use the URL with ogrinfo and ogr2ogr. In this test with your coordinates the URL is url-encoded because otherwise I GDAL was throwing an error.
ogrinfo "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" -al |more
INFO: Open of `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'
using driver `ESRIJSON' successful.
Layer name: ESRIJSON
Geometry: Polygon
Feature Count: 11978
Extent: (-76.120200, 38.049758) - (-75.678439, 38.420649)
Layer SRS WKT:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
BlockgroupID: String (12.0)
OGRFeature(ESRIJSON):0
BlockgroupID (String) = 240399302002
POLYGON ((-75.9405880296284 38.1716720350858,-75.9408220407599 38.1716180090955,-75.940838030772 38.1716600293136,-75.9406040196404 38.1717149733605,-75.9405880296284 38.1716720350858))
OGRFeature(ESRIJSON):1
BlockgroupID (String) = 240399302001
POLYGON ((-75.9427420100167 38.1518889956417,-75.9426950281273 38.1519089871405,-75.9426629582717 38.1518639886001,-75.9427100299926 38.151843997089,-75.9427420100167 38.1518889956417))
...
Realize this question is a bit old, but doing this currently using the arcgislayers library is so slick I simply have to share, as I've been struggling with different methods of using Microsoft Building footprints in R for hundreds of small areas:
library(arcgis)
library(arcgislayers)
library(sf)
library(mapview)
url <- "https://services.arcgis.com/P3ePLMYs2RVChkJx/ArcGIS/rest/services/MSBFP2/FeatureServer/0"
buildings <- arcgislayers::arc_open(url)
list_fields(buildings)
bb <- st_bbox(c(xmin = -76.12, ymin = 38.05, xmax = -75.68, ymax = 38.42 ), crs = st_crs(4326))
buildings_of_interest <- arcgislayers::arc_select(buildings, filter_geom = sf::st_bbox(bb))
mapview::mapview(buildings_of_interest,map.types = "OpenStreetMap")
I would love to get this working using rstac as well rather than using the ESRI hosted REST service.
st_readbut the exact query string to access within a bounding box needs creating...