[WIP]: Bandmath

README.md

@@ -90,6 +90,37 @@ $ curl "http://localhost:8000/tiles/14/3851/6812@2x?url=https%3A%2F%2Fs3-us-west
 
 ![greyscale stretched](docs/greyscale_stretched.png)
 
+
+### `/ndtiles/{z}/{x}/{y}` - Tiles
+
+Calculated normalized diffference index between two bands. The most common example is NDVI
+
+tile = (Band2 - Band1)/(Band2 _ Band 1)
+
+#### Parameters
+
+* `url` - a URL to a valid COG. Required.
+* `nodata` - a custom NODATA value.
+* `band1` - Band1 for Normalized Difference
+* `band2` - Band2 for Normalized Difference
+
+`@2x` can be added to the filename (after the `{y}` coordinate) to request
+retina tiles. The map preview will detect support for retina displays and
+request tiles accordingly.
+
+PNGs or JPEGs will be rendered depending on the presence of NODATA values in the
+source image (surfaced as transparency in the output).
+
+#### Examples
+
+```bash
+$ curl "http://localhost:8000/tiles/14/3851/6812@2x?url=https%3A%2F%2Fs3-us-west-2.amazonaws.com%2Fplanet-disaster-data%2Fhurricane-harvey%2FSkySat_Freeport_s03_20170831T162740Z3.tif" | imgcat
+```
+
+![RGB](docs/nd_test.jpeg)
+
+
+
 ### `/tiles` - TileJSON
 
 #### Parameters

virtual/catalogs.py

@@ -13,12 +13,17 @@
 
 class VirtualCatalog(Catalog):
 
-    def __init__(self, uri, rgb=None, nodata=None, linear_stretch=None):
+    def __init__(self, uri, rgb=None, nodata=None, linear_stretch=None, band1=None, band2=None):
         self._uri = uri
         self._rgb = rgb
         self._nodata = nodata
         self._linear_stretch = linear_stretch
         self._meta = {}
+        self._band1 = int(band1)
+        self._band2 = int(band2)
+
+        LOG.info("band1:{}".format(band1))
+        LOG.info("band2:{}".format(band2))
 
         with get_source(self._uri) as src:
             self._bounds = warp.transform_bounds(src.crs, WGS84_CRS, *src.bounds)
@@ -74,6 +79,10 @@ def get_sources(self, bounds, resolution):
         if self._linear_stretch is not None:
             recipes["linear_stretch"] = "per_band"
 
+        if (self._band1 is not None) and (self._band2 is not None):
+            recipes['nd_bands'] = (self._band1, self._band2)
+            LOG.info("nd_bands: {}".format(recipes['nd_bands']))
+
         yield Source(
             url=self._uri,
             name=self._name,

virtual/marblecutter_overload_nd.py

@@ -0,0 +1,157 @@
+# coding=utf-8
+from __future__ import absolute_import, division, print_function
+
+import logging
+
+import mercantile
+from rasterio.crs import CRS
+
+
+LOG = logging.getLogger(__name__)
+TILE_SHAPE = (256, 256)
+WEB_MERCATOR_CRS = CRS.from_epsg(3857)
+WGS84_CRS = CRS.from_epsg(4326)
+
+
+# coding=utf-8
+
+import unicodedata
+
+import numpy as np
+from rasterio.transform import Affine
+from virtual import mosaic
+from marblecutter.stats import Timer
+from marblecutter.utils import Bounds, PixelCollection
+from marblecutter import get_resolution_in_meters, NoDataAvailable
+
+
+
+def render_nd_tile(tile, catalog, transformation=None, format=None, scale=1, data_band_count=3):
+    """Render a tile into Web Mercator."""
+    bounds = Bounds(mercantile.xy_bounds(tile), WEB_MERCATOR_CRS)
+    shape = tuple(map(int, Affine.scale(scale) * TILE_SHAPE))
+
+    catalog.validate(tile)
+
+    return render_nd(
+        bounds,
+        shape,
+        WEB_MERCATOR_CRS,
+        catalog=catalog,
+        format=format,
+        data_band_count=data_band_count,
+        transformation=transformation,
+    )
+
+
+def performBandMath(pixels, catalog):
+    LOG.info("bandmath")
+    if (catalog.band1 is not None) and (catalog.band2 is not None):
+        LOG.info("Perform Bandmath")
+        band1 = int(catalog.band1) if isinstance(catalog.band1, str) else catalog.band1
+        band2 = int(catalog.band2) if isinstance(catalog.band2, str) else catalog.band2
+
+        LOG.info("pixels.data shape = {}".format(pixels.data.shape))
+        data = np.asarray(pixels.data)
+        LOG.info("data shape = {}".format(pixels.data.shape))
+        LOG.info("pixels = {}".format(pixels))
+
+
+        bandTile1 = data[:, :, band1]
+        bandTile2 = data[:, :, band2]
+
+        newBand = (bandTile2.astype(float) - bandTile1.astype(float)) / (bandTile2.astype(float) + bandTile1.astype(float))
+        newData = data
+        for idx in range(np.shape(data)[2]):
+            newData[:,:,idx] = newBand
+
+        np.seterr(divide='ignore', invalid='ignore')
+        new_pixel_collection =  PixelCollection(newData,
+                         pixels.bounds,
+                         pixels.band)
+
+    else:
+        ##TODO introduce new error
+        raise NoDataAvailable()
+
+    LOG.info("return Bandmath")
+    return new_pixel_collection
+
+def render_nd(
+    bounds,
+    shape,
+    target_crs,
+    format,
+    data_band_count,
+    catalog=None,
+    sources=None,
+    transformation=None
+):
+    """Render data intersecting bounds into shape using an optional
+    transformation. And perform normative diference using bands specified in catalog"""
+    resolution_m = get_resolution_in_meters(bounds, shape)
+    stats = []
+
+    if sources is None and catalog is None:
+        raise Exception("Either sources or a catalog must be provided.")
+
+    if transformation:
+        bounds, shape, offsets = transformation.expand(bounds, shape)
+
+    if sources is None and catalog is not None:
+        with Timer() as t:
+            sources = catalog.get_sources(bounds, resolution_m)
+        stats.append(("Get Sources", t.elapsed))
+
+    # TODO try to avoid materializing the iterator
+    sources = list(sources)
+
+    if sources is None or len(sources) == 0:
+        raise NoDataAvailable()
+
+    with Timer() as t:
+        sources_used, pixels = mosaic.composite(
+            sources, bounds, shape, target_crs, data_band_count
+        )
+    stats.append(("Composite", t.elapsed))
+
+    if pixels.data is None:
+        raise NoDataAvailable()
+
+    data_format = "raw"
+
+    if transformation:
+        with Timer() as t:
+            pixels, data_format = transformation.transform(pixels)
+        stats.append(("Transform", t.elapsed))
+
+        with Timer() as t:
+            pixels = transformation.postprocess(pixels, data_format, offsets)
+
+        stats.append(("Post-process", t.elapsed))
+
+    with Timer() as t:
+        (content_type, formatted) = format(pixels, data_format)
+    stats.append(("Format", t.elapsed))
+
+    headers = {
+        "Content-Type": content_type,
+        "Server-Timing": [
+            'op{};desc="{}";dur={:0.2f}'.format(i, name, time)
+            for (i, (name, time)) in enumerate(stats)
+        ]
+        + [
+            'src{};desc="{} - {}"'.format(
+                i,
+                unicodedata.normalize("NFKD", unicode(name)).encode(
+                    "ascii", "ignore"
+                ).replace(
+                    '"', '\\"'
+                ),
+                url,
+            )
+            for (i, (name, url)) in enumerate(sources_used)
+        ],
+    }
+
+    return (headers, formatted)

virtual/mosaic.py

@@ -0,0 +1,123 @@
+# noqa
+# coding=utf-8
+from __future__ import absolute_import, division, print_function
+
+import logging
+import multiprocessing
+from concurrent import futures
+
+import numpy as np
+
+from rasterio import warp
+
+from virtual import recipes
+from marblecutter.utils import Bounds, PixelCollection
+
+LOG = logging.getLogger(__name__)
+
+
+def composite(sources, bounds, shape, target_crs, band_count):
+    """Composite data from sources into a single raster covering bounds, but in
+    the target CRS."""
+    # avoid circular dependencies
+    from marblecutter import _nodata, get_resolution_in_meters, get_source, read_window
+
+    # TODO this belongs in render
+    if bounds.crs == target_crs:
+        canvas_bounds = bounds
+    else:
+        canvas_bounds = Bounds(
+            warp.transform_bounds(bounds.crs, target_crs, *bounds.bounds), target_crs
+        )
+
+    resolution = get_resolution_in_meters(bounds, shape)
+    sources = recipes.preprocess(sources, resolution=resolution)
+
+    def _read_window(source):
+        with get_source(source.url) as src:
+            LOG.info(
+                "Compositing %s (%s) as band %s", source.url, source.name, source.band
+            )
+
+            # read a window from the source data
+            # TODO ask for a buffer here, get back an updated bounding box
+            # reflecting it
+            try:
+                window_data = read_window(src, canvas_bounds, shape, source.recipes)
+            except Exception as e:
+                LOG.exception("Error reading %s: %s", source.url, e)
+                return
+
+            return source, PixelCollection(
+                window_data.data, window_data.bounds, source.band
+            )
+
+    # iterate over available sources, sorted by decreasing "quality"
+    with futures.ThreadPoolExecutor(
+        max_workers=multiprocessing.cpu_count() * 5
+    ) as executor:
+        ws = executor.map(_read_window, sources)
+
+    sources_used = []
+
+    ws = recipes.postprocess(ws)
+
+    canvas_data = np.ma.zeros(
+        (band_count,) + shape, dtype=np.float32, fill_value=_nodata(np.float32)
+    )
+    canvas_data.mask = True
+
+    canvas = PixelCollection(canvas_data, canvas_bounds)
+
+    for source, window_data in filter(None, ws):
+        window_data = recipes.apply(source.recipes, window_data, source=source)
+
+        if window_data.data is None:
+            continue
+
+        # paste the resulting data onto a canvas
+        canvas = paste(
+            PixelCollection(
+                window_data.data.astype(np.float32),
+                window_data.bounds,
+                window_data.band,
+            ),
+            canvas,
+        )
+        sources_used.append(source)
+
+        if not canvas.data.mask.any():
+            # stop if all pixels are valid
+            break
+
+    return map(lambda s: (s.name, s.url), sources_used), canvas
+
+
+def paste(window_pixels, canvas_pixels):
+    """ "Reproject" src data into the correct position within a larger image"""
+    window_data, (window_bounds, window_crs), band = window_pixels
+    canvas, (canvas_bounds, canvas_crs), _ = canvas_pixels
+    if window_crs != canvas_crs:
+        raise Exception("CRSes must match: {} != {}".format(window_crs, canvas_crs))
+
+    if window_bounds != canvas_bounds:
+        raise Exception(
+            "Bounds must match: {} != {}".format(window_bounds, canvas_bounds)
+        )
+
+    if band is None and window_data.shape != canvas.shape:
+        raise Exception(
+            "Data shapes must match: {} != {}".format(window_data.shape, canvas.shape)
+        )
+
+    if band is None:
+        merged = np.ma.where(canvas.mask & ~window_data.mask, window_data, canvas)
+        merged.fill_value = canvas.fill_value
+    else:
+        merged_band = np.ma.where(
+            canvas.mask[band] & ~window_data.mask, window_data, canvas[band]
+        )
+        canvas[band] = merged_band[0]
+        merged = canvas
+
+    return PixelCollection(merged, canvas_pixels.bounds)

virtual/nd_formats.py

@@ -0,0 +1,36 @@
+# coding=utf-8
+from __future__ import absolute_import
+
+import logging
+
+import numpy as np
+
+from marblecutter import _isimage
+from marblecutter.utils import PixelCollection
+from marblecutter.formats.jpeg import JPEG
+from marblecutter.formats.png import PNG
+
+LOG = logging.getLogger(__name__)
+JPEG_FORMAT = JPEG()
+PNG_FORMAT = PNG()
+
+
+def Optimalnd():
+
+    def _format(pixels, data_format):
+        if not _isimage(data_format):
+            raise Exception("Must be an image format")
+
+        alpha = pixels.data[:, :, 3]
+
+        if np.equal(alpha, 255).all():
+            # solid
+            rgb_pixels = PixelCollection(
+                pixels.data[:, :, 0:3], pixels.bounds, pixels.band
+            )
+            return JPEG_FORMAT(rgb_pixels, "RGB")
+        else:
+            # partially transparent
+            return PNG_FORMAT(pixels, data_format)
+
+    return _format