Ms-buildings within a bounding-box

[Daniele Mainardi (dmainardi96)]

We want to search for building footprints in a certain area using pystac-client with the "ms-buildings" collection. We need to create a .geojson tile with the footprints giving a restricted tile, like a city street.
This is the python code we use:

from pystac_client import Client
import planetary_computer as pc

# Search against the Planetary Computer STAC API
catalog = Client.open(
  "https://planetarycomputer.microsoft.com/api/stac/v1"
)

# Define your area of interest
aoi = {
  "type": "Polygon",
  "coordinates": [
    bounding_box
  ]
}

# Define your search with CQL2 syntax
search = catalog.search(filter_lang="cql2-json", filter={
  "op": "and",
  "args": [
    {"op": "s_intersects", "args": [{"property": "geometry"}, aoi]},
    {"op": "=", "args": [{"property": "collection"}, "ms-buildings"]}
  ]
})

# Grab the first item from the search results and sign the assets
first_item = next(search.get_items())
pc.sign_item(first_item).assets

asset = first_item.assets["data"]
df = geopandas.read_parquet(asset.href, storage_options=asset.extra_fields["table:storage_options"])

The problem is that the asset contains the url for the whole country and the geopandas read_parquet method takes a long time to scan every building. Is there a way to query only the buildings inside the tile?

Any help, hints, or guidance would be greatly appreciated!

[Tom Augspurger (TomAugspurger)]

the geopandas read_parquet method takes a long time to scan every building. Is there a way to query only the buildings inside the tile?

If you use an operation like clip, then the spatial partitioning should ensure that only data from the tiles intersecting with your AOI are read (see "Working with large files" at https://planetarycomputer.microsoft.com/dataset/ms-buildings#Example-Notebook).

Do you have a sense for how many partitions of the parquet dataset intersect with your area of interest?

[Daniele Mainardi (dmainardi96)]

This is an example of an area of interest I need.

I wanted to make a script that given a point, returns the building footprints in a surrounding area of approximately 150m x 150m.
The problem is that i need this operation done within ~5 seconds. I also tried to use dask_geopandas to read the parquet, but it's too slow for my purpose, since my area of interest intersects with the polygon of the entire nation, but i need only a very small portion of it.

I wonder if there's any tool i could use that can directly give me the building polygons in a small AoI in the fastest way possible.
Thank you for the answer!

[Tom Augspurger (TomAugspurger)]

It might be challenging to do that with how the files are currently laid out. At a minimum, you'll need to

1. Read the bounding boxes of each partition in the parquet dataset (one HTTP request per partition)
2. Read all the data of partitions overlapping with your AOI

In theory, you could skip the first stage by finding the quadkeys your AOI intersects, and use that to know which partitions to load up. But overall, it sounds like the layout on Blob Storage isn't really aligned with your access pattern.

[Tom Augspurger (TomAugspurger)]

Following up here, the most recent update to the ms-buildings collection should make this easier. The data are now partitioned by (RegionName, quadkey). You can either use the Delta Table metadata or the STAC API to query by quadkey / geometry. https://planetarycomputer.microsoft.com/dataset/ms-buildings#Example-Notebook has examples of both.