From 13be9f5a04cb711a22d11cf5de832ca75867f3df Mon Sep 17 00:00:00 2001 From: teslov Date: Thu, 7 Oct 2021 20:52:09 +0300 Subject: [PATCH 1/7] Added support of retrieve metadata for cutline specified as geojson feature or polygon in boundaries header --- app/api/tiler.py | 29 ++++++++++++++++++++++------- 1 file changed, 22 insertions(+), 7 deletions(-) diff --git a/app/api/tiler.py b/app/api/tiler.py index 78b2e85f..639423e8 100644 --- a/app/api/tiler.py +++ b/app/api/tiler.py @@ -1,13 +1,15 @@ import json +import base64 import numpy import rio_tiler.utils from rasterio.enums import ColorInterp +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,14 @@ 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.headers.get('boundaries') if formula == '': formula = None if bands == '': bands = None if defined_range == '': defined_range = None + if boundaries_feature == '': boundaries_feature = None + else: + boundaries_feature = json.load(boundaries_feature) try: expr, hrange = lookup_formula(formula, bands) if defined_range is not None: @@ -144,21 +150,26 @@ 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, geometry_crs="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: - data, mask = src.preview(expression=expr) + 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]) @@ -166,7 +177,11 @@ class Metadata(TaskNestedView): stats = {b: ImageStatistics(**s) for b, s in stats.items()} metadata = RioMetadata(statistics=stats, **src.info().dict()) else: - metadata = src.metadata(pmin=pmin, pmax=pmax, hist_options=histogram_options, nodata=nodata) + 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_bbox}) + else: + metadata = src.metadata(pmin=pmin, pmax=pmax, hist_options=histogram_options, nodata=nodata) info = json.loads(metadata.json()) except IndexError as e: # Caught when trying to get an invalid raster metadata From de35720e4003797ca5ca690a5d9344830ae1590d Mon Sep 17 00:00:00 2001 From: teslov Date: Thu, 7 Oct 2021 20:59:28 +0300 Subject: [PATCH 2/7] Added getting boundaries for tile as query param --- app/api/tiler.py | 33 ++++++++++++++++++++++++++++----- 1 file changed, 28 insertions(+), 5 deletions(-) diff --git a/app/api/tiler.py b/app/api/tiler.py index 639423e8..7bfa193e 100644 --- a/app/api/tiler.py +++ b/app/api/tiler.py @@ -297,6 +297,13 @@ 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 + else: + boundaries_feature = json.load(boundaries_feature) + if formula == '': formula = None if bands == '': bands = None @@ -339,6 +346,12 @@ 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, geometry_crs="epsg:4326") + boundaries_bbox = featureBounds(boundaries_feature) + else: + boundaries_cutline = None + boundaries_bbox = None # Handle N-bands datasets for orthophotos (not plant health) if tile_type == 'orthophoto' and expr is None: ci = src.dataset.colorinterp @@ -370,12 +383,22 @@ class Tiles(TaskNestedView): try: with COGReader(url) as src: if expr is not None: - tile = src.tile(x, y, z, expression=expr, tilesize=tilesize, nodata=nodata, - padding=padding, - resampling_method=resampling) + 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: - tile = src.tile(x, y, z, indexes=indexes, tilesize=tilesize, nodata=nodata, - padding=padding, resampling_method=resampling) + 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) except TileOutsideBounds: raise exceptions.NotFound("Outside of bounds") From f5903f74b51f816115b6088fe085822c037f8016 Mon Sep 17 00:00:00 2001 From: teslov Date: Thu, 7 Oct 2021 21:03:48 +0300 Subject: [PATCH 3/7] Changed header to query params at metadata --- app/api/tiler.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/app/api/tiler.py b/app/api/tiler.py index 7bfa193e..f09616d8 100644 --- a/app/api/tiler.py +++ b/app/api/tiler.py @@ -123,7 +123,7 @@ 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.headers.get('boundaries') + boundaries_feature = self.request.query_params.get('boundaries') if formula == '': formula = None if bands == '': bands = None From fad63797aff60ef8c67be438d697a5e007f2c6c4 Mon Sep 17 00:00:00 2001 From: teslov Date: Thu, 7 Oct 2021 21:04:49 +0300 Subject: [PATCH 4/7] Fixed json.load to loads function --- app/api/tiler.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/app/api/tiler.py b/app/api/tiler.py index f09616d8..904deb6f 100644 --- a/app/api/tiler.py +++ b/app/api/tiler.py @@ -130,7 +130,7 @@ class Metadata(TaskNestedView): if defined_range == '': defined_range = None if boundaries_feature == '': boundaries_feature = None else: - boundaries_feature = json.load(boundaries_feature) + boundaries_feature = json.loads(boundaries_feature) try: expr, hrange = lookup_formula(formula, bands) if defined_range is not None: @@ -302,7 +302,7 @@ class Tiles(TaskNestedView): if boundaries_feature == '': boundaries_feature = None else: - boundaries_feature = json.load(boundaries_feature) + boundaries_feature = json.loads(boundaries_feature) if formula == '': formula = None From 8deed3b0eca0ad3698f96db797c6ef85560c6f37 Mon Sep 17 00:00:00 2001 From: teslov Date: Thu, 7 Oct 2021 21:09:41 +0300 Subject: [PATCH 5/7] Fixed errors during creating cutline --- app/api/tiler.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/app/api/tiler.py b/app/api/tiler.py index 904deb6f..5bc4fa9b 100644 --- a/app/api/tiler.py +++ b/app/api/tiler.py @@ -124,12 +124,11 @@ class Metadata(TaskNestedView): 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 - else: + if boundaries_feature is not None: boundaries_feature = json.loads(boundaries_feature) try: expr, hrange = lookup_formula(formula, bands) @@ -151,7 +150,7 @@ class Metadata(TaskNestedView): 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, geometry_crs="epsg:4326") + boundaries_cutline = create_cutline(src.dataset, geometry=boundaries_feature, geometry_crs="epsg:4326") boundaries_bbox = featureBounds(boundaries_feature) else: boundaries_cutline = None @@ -298,10 +297,9 @@ class Tiles(TaskNestedView): 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 - else: + if boundaries_feature is not None: boundaries_feature = json.loads(boundaries_feature) @@ -347,11 +345,9 @@ class Tiles(TaskNestedView): 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, geometry_crs="epsg:4326") - boundaries_bbox = featureBounds(boundaries_feature) + boundaries_cutline = create_cutline(src.dataset, geometry=boundaries_feature, geometry_crs="epsg:4326") else: boundaries_cutline = None - boundaries_bbox = None # Handle N-bands datasets for orthophotos (not plant health) if tile_type == 'orthophoto' and expr is None: ci = src.dataset.colorinterp From 542077bdbfd7045c92780c39ef3ea93c2311ccce Mon Sep 17 00:00:00 2001 From: teslov Date: Thu, 7 Oct 2021 21:47:25 +0300 Subject: [PATCH 6/7] Fixed misstypes --- app/api/tiler.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/app/api/tiler.py b/app/api/tiler.py index 5bc4fa9b..06a12585 100644 --- a/app/api/tiler.py +++ b/app/api/tiler.py @@ -3,6 +3,7 @@ import base64 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 @@ -150,7 +151,7 @@ class Metadata(TaskNestedView): with COGReader(raster_path) as src: band_count = src.metadata()['count'] if boundaries_feature is not None: - boundaries_cutline = create_cutline(src.dataset, geometry=boundaries_feature, geometry_crs="epsg:4326") + boundaries_cutline = create_cutline(src.dataset, boundaries_feature, CRS.from_string('EPSG:4326')) boundaries_bbox = featureBounds(boundaries_feature) else: boundaries_cutline = None @@ -178,7 +179,7 @@ class Metadata(TaskNestedView): 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_bbox}) + , 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()) @@ -216,7 +217,6 @@ class Metadata(TaskNestedView): info['color_maps'] = [] info['algorithms'] = algorithms - if colormaps: for cmap in colormaps: try: @@ -345,7 +345,7 @@ class Tiles(TaskNestedView): 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, geometry=boundaries_feature, geometry_crs="epsg:4326") + 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) From d563052d07de50f7de233afae77d043c2ac32b2a Mon Sep 17 00:00:00 2001 From: teslov Date: Fri, 8 Oct 2021 09:53:29 +0300 Subject: [PATCH 7/7] Removed unused import --- app/api/tiler.py | 1 - 1 file changed, 1 deletion(-) diff --git a/app/api/tiler.py b/app/api/tiler.py index 06a12585..7e84bf0c 100644 --- a/app/api/tiler.py +++ b/app/api/tiler.py @@ -1,5 +1,4 @@ import json -import base64 import numpy import rio_tiler.utils from rasterio.enums import ColorInterp