Merge pull request #1068 from AIDI-solar/cutline

Added support of getting tiles and metadata with geoJSON polygon specified at request param

app/api/tiler.py

@@ -2,12 +2,14 @@ import json
 import numpy
 import rio_tiler.utils
 from rasterio.enums import ColorInterp
+from rasterio.crs import CRS
+from rasterio.features import bounds as featureBounds
 import urllib
 import os
 from django.http import HttpResponse
 from rio_tiler.errors import TileOutsideBounds
 from rio_tiler.utils import has_alpha_band, \
-    non_alpha_indexes, render, mapzen_elevation_rgb
+    non_alpha_indexes, render, create_cutline
 from rio_tiler.utils import _stats as raster_stats
 from rio_tiler.models import ImageStatistics, ImageData
 from rio_tiler.models import Metadata as RioMetadata
@@ -121,10 +123,13 @@ class Metadata(TaskNestedView):
         formula = self.request.query_params.get('formula')
         bands = self.request.query_params.get('bands')
         defined_range = self.request.query_params.get('range')
-
+        boundaries_feature = self.request.query_params.get('boundaries')
         if formula == '': formula = None
         if bands == '': bands = None
         if defined_range == '': defined_range = None
+        if boundaries_feature == '': boundaries_feature = None
+        if boundaries_feature is not None:
+            boundaries_feature = json.loads(boundaries_feature)
         try:
             expr, hrange = lookup_formula(formula, bands)
             if defined_range is not None:
@@ -144,27 +149,36 @@ class Metadata(TaskNestedView):
         try:
             with COGReader(raster_path) as src:
                 band_count = src.metadata()['count']
+                if boundaries_feature is not None:
+                    boundaries_cutline = create_cutline(src.dataset, boundaries_feature, CRS.from_string('EPSG:4326'))
+                    boundaries_bbox = featureBounds(boundaries_feature)
+                else:
+                    boundaries_cutline = None
+                    boundaries_bbox = None
                 if has_alpha_band(src.dataset):
                     band_count -= 1
                 nodata = None
                 # Workaround for https://github.com/OpenDroneMap/WebODM/issues/894
                 if tile_type == 'orthophoto':
                     nodata = 0
-                # info = src.metadata(pmin=pmin, pmax=pmax, histogram_bins=255, histogram_range=hrange, expr=expr,
-                #                      nodata=nodata)
                 histogram_options = {"bins": 255, "range": hrange}
-
                 if expr is not None:
+                    if boundaries_cutline is not None:
+                        data, mask = src.preview(expression=expr, vrt_options={'cutline': boundaries_cutline})
+                    else:
                         data, mask = src.preview(expression=expr)
                     data = numpy.ma.array(data)
                     data.mask = mask == 0
-                    expression_bloc = expr.lower().split(",")
                     stats = {
                         str(b + 1): raster_stats(data[b], percentiles=(pmin, pmax), bins=255, range=hrange)
                         for b in range(data.shape[0])
                     }
                     stats = {b: ImageStatistics(**s) for b, s in stats.items()}
                     metadata = RioMetadata(statistics=stats, **src.info().dict())
+                else:
+                    if (boundaries_cutline is not None) and (boundaries_bbox is not None):
+                        metadata = src.metadata(pmin=pmin, pmax=pmax, hist_options=histogram_options, nodata=nodata
+                                                , bounds=boundaries_bbox, vrt_options={'cutline': boundaries_cutline})
                     else:
                         metadata = src.metadata(pmin=pmin, pmax=pmax, hist_options=histogram_options, nodata=nodata)
                 info = json.loads(metadata.json())
@@ -202,7 +216,6 @@ class Metadata(TaskNestedView):
 
         info['color_maps'] = []
         info['algorithms'] = algorithms
-
         if colormaps:
             for cmap in colormaps:
                 try:
@@ -282,6 +295,12 @@ class Tiles(TaskNestedView):
         rescale = self.request.query_params.get('rescale')
         color_map = self.request.query_params.get('color_map')
         hillshade = self.request.query_params.get('hillshade')
+        boundaries_feature = self.request.query_params.get('boundaries')
+        if boundaries_feature == '':
+            boundaries_feature = None
+        if boundaries_feature is not None:
+            boundaries_feature = json.loads(boundaries_feature)
+
 
         if formula == '': formula = None
         if bands == '': bands = None
@@ -324,6 +343,10 @@ class Tiles(TaskNestedView):
             has_alpha = has_alpha_band(src.dataset)
             if z < minzoom - ZOOM_EXTRA_LEVELS or z > maxzoom + ZOOM_EXTRA_LEVELS:
                 raise exceptions.NotFound()
+            if boundaries_feature is not None:
+                boundaries_cutline = create_cutline(src.dataset, boundaries_feature, CRS.from_string('EPSG:4326'))
+            else:
+                boundaries_cutline = None
             # Handle N-bands datasets for orthophotos (not plant health)
             if tile_type == 'orthophoto' and expr is None:
                 ci = src.dataset.colorinterp
@@ -355,9 +378,19 @@ class Tiles(TaskNestedView):
         try:
             with COGReader(url) as src:
                 if expr is not None:
+                    if boundaries_cutline is not None:
+                        tile = src.tile(x, y, z, expression=expr, tilesize=tilesize, nodata=nodata,
+                                        padding=padding,
+                                        resampling_method=resampling, vrt_options={'cutline': boundaries_cutline})
+                    else:
                         tile = src.tile(x, y, z, expression=expr, tilesize=tilesize, nodata=nodata,
                                         padding=padding,
                                         resampling_method=resampling)
+                else:
+                    if boundaries_cutline is not None:
+                        tile = src.tile(x, y, z, tilesize=tilesize, nodata=nodata,
+                                        padding=padding,
+                                        resampling_method=resampling, vrt_options={'cutline': boundaries_cutline})
                     else:
                         tile = src.tile(x, y, z, indexes=indexes, tilesize=tilesize, nodata=nodata,
                                         padding=padding, resampling_method=resampling)