Skip to content

Add Support for SpatioTemporal Asset Catalogs #21

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
DanielJDufour wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Add Support for SpatioTemporal Asset Catalogs #21

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 14 additions & 0 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,6 +136,20 @@ See tile parameters.

`http://localhost:8000/preview?url=https%3A%2F%2Fs3-us-west-2.amazonaws.com%2Fplanet-disaster-data%2Fhurricane-harvey%2FSkySat_Freeport_s03_20170831T162740Z3.tif`

### `/stac` - SpatioTemporal Asset Catalog
Generate a tile from the automatic mosaicing of STAC assets

#### Parameters
* `url` - a URL to a STAC API Feature Collection of STAC Features
* `datetime` - an optional datetime filter to pass along to the STAC API
* `expr` - an optional band math expression where a band is referenced by `name[band_number]`. For example `Analytic[0]` would refer to the first band of the assset called `Analytic`.
* `jq` - an optional [jq](https://stedolan.github.io/jq/) expression for filtering and ordering the STAC catalog search results. For example `sort_by(.properties.datetime) | reverse | group_by(.properties["eo:row"] + .properties["eo:column"]) | map(.[0])` is used to get the newest capture for each Landsat 8 scene.

#### Known Limitations of STAC Support
- Mask is not carried over when doing band expressions across multiple assets
- Band math does not work across assets with different resolutions


## Deploying to AWS

marblecutter-virtual is deployed using the [AWS Serverless Application Model
Expand Down
6 changes: 5 additions & 1 deletion requirements.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,9 +1,13 @@
affine
cachetools ~= 2.0.0
flask-cors
jq
# marblecutter[web] ~= 0.3.1
https://github.com/mojodna/marblecutter/archive/5b9040b.tar.gz#egg=marblecutter[web]
rasterio[s3] ~= 1.0
mercantile
numpy
rasterio[s3] ~= 1.0
shapely
serverless-wsgi
# temporary workaround for https://github.com/mapbox/rasterio/issues/1651
mock
16 changes: 16 additions & 0 deletions templates/test.html
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
<html>

<body style="background: red">
<!-- Red, Green, Blue -->
<img src="/stac/14/4096/6048?url=https://sat-api.developmentseed.org/collections/landsat-8-l1/items&expr=B4,B3,B2&jq=sort_by(.properties.datetime)%20%7C%20reverse%20%7C%20group_by(.properties%5B%22eo%3Arow%22%5D%20%2B%20.properties%5B%22eo%3Acolumn%22%5D)%20%7C%20map(.%5B0%5D)">

<!-- NIR, Green, Blue -->
<img src="/stac/14/4096/6048?url=https://sat-api.developmentseed.org/collections/landsat-8-l1/items&expr=B5,B3,B2&jq=sort_by(.properties.datetime)%20%7C%20reverse%20%7C%20group_by(.properties%5B%22eo%3Arow%22%5D%20%2B%20.properties%5B%22eo%3Acolumn%22%5D)%20%7C%20map(.%5B0%5D)">

<!-- NDVI -->
<img src="/stac/14/4096/6048?url=https://sat-api.developmentseed.org/collections/landsat-8-l1/items&expr=(B5[0]-B4)%2F(B5[0]%2BB4[0])&jq=sort_by(.properties.datetime)%20%7C%20reverse%20%7C%20group_by(.properties%5B%22eo%3Arow%22%5D%20%2B%20.properties%5B%22eo%3Acolumn%22%5D)%20%7C%20map(.%5B0%5D)">

<!-- NDVI with First-Band for each Asset Implied -->
<img src="/stac/14/4096/6048?url=https://sat-api.developmentseed.org/collections/landsat-8-l1/items&expr=(B5-B4)%2F(B5%2BB4)&jq=sort_by(.properties.datetime)%20%7C%20reverse%20%7C%20group_by(.properties%5B%22eo%3Arow%22%5D%20%2B%20.properties%5B%22eo%3Acolumn%22%5D)%20%7C%20map(.%5B0%5D)">
</body>
</html>
197 changes: 195 additions & 2 deletions virtual/web.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,28 @@

import logging

from affine import Affine
from cachetools.func import lru_cache
from collections import defaultdict
from concurrent import futures
from datetime import datetime
from flask import Flask, Markup, jsonify, redirect, render_template, request
from flask_cors import CORS
from itertools import groupby
from marblecutter import WEB_MERCATOR_CRS, get_source, read_window
from marblecutter.tiling import TILE_SHAPE
from marblecutter.mosaic import MAX_WORKERS
from marblecutter import NoCatalogAvailable, tiling
from marblecutter.formats.optimal import Optimal
from marblecutter.transformations import Image
from marblecutter.web import bp, url_for
from mercantile import Tile
from marblecutter.utils import Bounds, Source, PixelCollection
import mercantile
import numpy as np
import jq
import re
import requests
from shapely.geometry import box

try:
from urllib.parse import urlparse, urlencode
Expand All @@ -33,6 +47,14 @@
app.url_map.strict_slashes = False
CORS(app, send_wildcard=True)

class Timer:
def __init__(self, description):
self.description = description
def __enter__(self, *args):
self.start = datetime.now()
def __exit__(self, *args):
duration = (datetime.now() - self.start).total_seconds()
LOG.info(f'{self.description} took {duration} seconds')

@lru_cache()
def make_catalog(args):
Expand Down Expand Up @@ -85,6 +107,19 @@ def bounds():

return jsonify({"url": catalog.uri, "bounds": catalog.bounds})

@app.route("/test")
def test():
return (
render_template(
"test.html",
tilejson_url=Markup(
url_for("meta", _external=True, _scheme="", **request.args)
),
),
200,
{"Content-Type": "text/html"},
)


@app.route("/preview")
def preview():
Expand All @@ -106,12 +141,170 @@ def preview():
{"Content-Type": "text/html"},
)

@app.route("/stac/<int:z>/<int:x>/<int:y>")
@app.route("/stac/<int:z>/<int:x>/<int:y>@<int:scale>x")
def render_png_from_stac_catalog(z, x, y, scale=1):
with Timer("rendering png from stac catalog"):
stac_url = request.args.get("url", None)
jq_filter = request.args.get("jq", None)
stac_expr = request.args.get("expr", None)
stac_datetime = request.args.get("datetime", None)

# size of the tile, usually (256, 256)
shape = tuple(map(int, Affine.scale(scale) * TILE_SHAPE))

if stac_expr:
# captures asset-band combos
# like B5[0] in (B5[0] - B4[0])/(B5[0] + B4[0])
# or like NIR in (NIR - RED) / (NIR + RED)
asset_band_regex = "(?P<asset>[A-Za-z][A-Za-z0-9]+)(?:\[(?P<band>\d+)\])?"
matches = list(set(re.findall(asset_band_regex, stac_expr)))

# sorted list of assets and bands
# like [('B4', 0), ('B5', 0)]
# or like [('NIR', 0), ('RED', 0)]
asset_bands = sorted(list(set([(asset, int(band) if band else 0) for asset, band in matches])))

# sorted list of asset names
# like ['B4', 'B5']
# or like ['NIR', 'RED']
asset_names = sorted(list(set([asset for asset, band in asset_bands])))

# convert expr from a format for running band math across multiple assets
# into a format for running band math for a single file for the combined assets
# from: (NIR - RED) / (NIR + RED)
# to: (b1 - b2) / (b1 + b2)
def repl(m):
asset, band = m.groups()
band = int(band) if band else 0
# add one to index number because single-file band math expression
# requires that band indexes starts at 1, i.e. b1, b2, b3...
return 'b' + str(asset_bands.index((asset, band)) + 1)
expr = re.sub(asset_band_regex, repl, stac_expr)
else:
asset_names = None
expr = None

tile = mercantile.Tile(x, y, z)

tile_bounds = mercantile.bounds(tile)
tile_bbox = [tile_bounds.west, tile_bounds.south, tile_bounds.east, tile_bounds.north]

# we use the parent tile for searching because sometimes
# a search engine might not return results
# when the tile is really small
parent_tile = mercantile.parent(tile)
search_bounds = mercantile.bounds(parent_tile)

search_bbox = [
search_bounds.west,
search_bounds.south,
search_bounds.east,
search_bounds.north
]

tile_polygon = box(*tile_bbox)

params = {
'bbox': str(search_bbox).replace(' ', ''),
'limit': 500,
}
if stac_datetime: params['datetime'] = stac_datetime

with Timer("querying stac"):
response = requests.get(stac_url, params=params)

features = response.json()['features']
LOG.info(f'number of features: {len(features)}')

# filter features to those that overlap tile
features = [feature for feature in features if box(*feature['bbox']).intersects(tile_polygon)]

feature_count = len(features)
LOG.info(f'number of features after filtering by feature extent: {feature_count}')

if jq_filter:
features = jq.compile(jq_filter).input(features).first()
LOG.info(f'number of features after filtering by jq expression: {len(features)}')

canvas_bounds = Bounds(bounds=mercantile.xy_bounds(tile), crs=WEB_MERCATOR_CRS)
LOG.info(f'canvas bounds: {canvas_bounds}')

assets = []
for fid, feature in enumerate(features):
images = {}
if asset_names:
for asset_name in asset_names:
images[asset_name] = feature['assets'][asset_name]['href']
elif 'visual' in feature['assets']:
images['visual'] = feature['assets']['visual']['href']
else:
raise "Not sure what assets to use to create the image"

for asset_name, href in images.items():
assets.append({
"fid": fid,
"name": asset_name,
"url": href
})

if expr and len(asset_names) > 0:
def add_pixels_to_asset(asset):
try:
url = asset['url']
with Timer(f'reading pixels for {url}'):
source = Source(url=url, name=url, resolution=None)
with get_source(url) as src:
with Timer(f'reading window for {url}'):
asset['pixels'] = read_window(src, canvas_bounds, shape, source)
except Exception as e:
LOG.error(e)
raise e

with Timer(f'reading all the pixels'):
with futures.ThreadPoolExecutor(max_workers=MAX_WORKERS) as executor:
executor.map(add_pixels_to_asset, assets)

sources = []
for fid, assets in groupby(assets, lambda asset: asset['fid']):
assets = list(assets)
name_to_pixels = dict([(asset['name'], asset['pixels']) for asset in assets])
windows = tuple([name_to_pixels[name].data[band] for name, band in asset_bands])
stacked = np.ma.stack(windows)
names = [asset['name'] for asset in assets]
source = Source(
url=None,
name=str(fid) + '{' + ','.join(names) + '}',
resolution=None,
expr=expr,
pixels=PixelCollection(stacked, canvas_bounds),
recipes={ "expr": expr, "imagery": True } if expr else {}
)
sources.append(source)
else:
sources = [Source(
url=asset['url'],
name=asset['name'],
resolution=None,
recipes={ "expr": expr, "imagery": True } if expr else {}
) for asset in assets]

headers, data = tiling.render_tile_from_sources(
tile,
sources,
format=IMAGE_FORMAT,
transformation=IMAGE_TRANSFORMATION,
scale=scale,
)

return data, 200, headers


@app.route("/tiles/<int:z>/<int:x>/<int:y>")
@app.route("/tiles/<int:z>/<int:x>/<int:y>@<int:scale>x")
def render_png(z, x, y, scale=1):
catalog = make_catalog(request.args)
tile = Tile(x, y, z)
tile = mercantile.Tile(x, y, z)

headers, data = tiling.render_tile(
tile,
Expand Down