2

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.

3
  • 1
    Via the ArcGIS REST Server connection in at services.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services I can see (in QGIS) lots of layers including MSBFP. I can select "only request features overlapping the view extent" and I see buildings, so that must be subsetting by box. If you can fabricate a query for that REST server with a bounding box you should be able to get a subset....
    – Spacedman
    Dec 12, 2021 at 21:52
  • 1
    I think its accessible via the ESRIJSON OGR driver gdal.org/drivers/vector/esrijson.html which means you might be able to get it with st_read but the exact query string to access within a bounding box needs creating...
    – Spacedman
    Dec 12, 2021 at 22:06
  • It seems like this approach would work with any of the layers on the ArcGIS Rest server. Maybe it would be helpful to users to point that out. I could add it in the question or just ask another generalized question about downloading any layer using R?
    – Cannon
    Dec 14, 2021 at 8:54

2 Answers 2

3

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...
2

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:

https://services.arcgis.com/P3ePLMYs2RVChkJx/ArcGIS/rest/services/MSBFP2/FeatureServer/0/query?f=json&returnGeometry=true&spatialRel=esriSpatialRelIntersects&geometry={"xmin":-76.12,"ymin":38.05,"xmax":-75.68,"ymax":38.42,"spatialReference":{"wkid":4326}}&geometryType=esriGeometryEnvelope&outSR=4326

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))
...
2
  • 1
    I've coded that string into an R solution!
    – Spacedman
    Dec 13, 2021 at 13:57
  • Thanks for the help. I would also accept your answer if picking more than one was allowed.
    – Cannon
    Dec 14, 2021 at 8:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.