From 38149e527659b906c595ae106794c60e05b1fb5c Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Fri, 15 Nov 2019 16:52:23 -0500 Subject: [PATCH] Hillshade POC, tile artifacts present --- app/api/tiler.py | 67 +++++++++++++++++++--------- app/static/app/js/components/Map.jsx | 4 +- requirements.txt | 3 +- 3 files changed, 51 insertions(+), 23 deletions(-) diff --git a/app/api/tiler.py b/app/api/tiler.py index b8392f18..0064ecc4 100644 --- a/app/api/tiler.py +++ b/app/api/tiler.py @@ -5,13 +5,12 @@ from rasterio import MemoryFile from rio_tiler.errors import TileOutsideBounds from rio_tiler.mercator import get_zooms from rio_tiler import main -from rio_tiler.utils import array_to_image, get_colormap, expression, linear_rescale, _chunks -from rio_color.operations import parse_operations -from rio_color.utils import scale_dtype, to_math_type +from rio_tiler.utils import array_to_image, get_colormap, expression, linear_rescale, _chunks, _apply_discrete_colormap from rio_tiler.profiles import img_profiles import numpy import mercantile +from matplotlib.colors import LightSource from .tasks import TaskNestedView from rest_framework import exceptions @@ -22,7 +21,7 @@ def get_tile_url(task, tile_type, query_params): url = '/api/projects/{}/tasks/{}/{}/tiles/{{z}}/{{x}}/{{y}}.png'.format(task.project.id, task.id, tile_type) params = {} - for k in ['expr', 'rescale', 'color_map']: + for k in ['expr', 'rescale', 'color_map', 'hillshade']: if query_params.get(k): params[k] = query_params.get(k) @@ -53,7 +52,7 @@ def get_raster_path(task, tile_type): return task.get_asset_download_path(tile_type + ".tif") -def postprocess(tile, mask, rescale = None, color_formula = None): +def rescale_tile(tile, mask, rescale = None): if rescale: rescale_arr = list(map(float, rescale.split(","))) rescale_arr = list(_chunks(rescale_arr, 2)) @@ -69,16 +68,17 @@ def postprocess(tile, mask, rescale = None, color_formula = None): ) tile = tile.astype(numpy.uint8) - if color_formula: - # make sure one last time we don't have - # negative value before applying color formula - tile[tile < 0] = 0 - for ops in parse_operations(color_formula): - tile = scale_dtype(ops(to_math_type(tile)), numpy.uint8) - return tile, mask +def apply_colormap(tile, color_map = None): + if color_map is not None and isinstance(color_map, dict): + tile = _apply_discrete_colormap(tile, color_map) + elif color_map is not None: + tile = numpy.transpose(color_map[tile][0], [2, 0, 1]).astype(numpy.uint8) + + return tile + class TileJson(TaskNestedView): def get(self, request, pk=None, project_pk=None, tile_type=""): """ @@ -138,7 +138,7 @@ class Metadata(TaskNestedView): raster_path, coords.x, coords.y, coords.z, expr=expr, tilesize=256, nodata=None ) - rtile, rmask = postprocess(tile, mask, rescale=rescale) + rtile, rmask = rescale_tile(tile, mask, rescale) del tile del mask @@ -190,12 +190,12 @@ class Tiles(TaskNestedView): #nodata = self.request.query_params.get('nodata') indexes = None - color_formula = None nodata = None expr = self.request.query_params.get('expr') rescale = self.request.query_params.get('rescale') color_map = self.request.query_params.get('color_map') + hillshade = self.request.query_params.get('hillshade') is not None # TODO: disable color_map # TODO: server-side expressions @@ -225,20 +225,45 @@ class Tiles(TaskNestedView): if tile.shape[0] == 4: mask = None - rtile, rmask = postprocess( - tile, mask, rescale=rescale, color_formula=color_formula - ) - del tile - del mask - if color_map: try: color_map = get_colormap(color_map, format="gdal") except FileNotFoundError: raise exceptions.ValidationError("Not a valid color_map value") + rtile, rmask = rescale_tile(tile, mask, rescale=rescale) + + if hillshade: + if tile.shape[0] != 1: + raise exceptions.ValidationError("Cannot compute hillshade of non-elevation raster (multiple bands found)") + + with rasterio.open(url) as src: + dx = src.meta["transform"][0] + dy = -src.meta["transform"][4] + + ls = LightSource() + elevation = tile[0] # First band + rgb = apply_colormap(rtile, color_map) + + # BxMxN (0..255) --> MxNxB (0..1) + rgb = rtile.reshape(rgb.shape[1], rgb.shape[2], rgb.shape[0]) / 255.0 + + rtile = ls.shade_rgb(rgb, + elevation, + vert_exag=1, + dx=dx, + dy=dy, + blend_mode="overlay" + ) + del rgb + + # MxNxB (0..1) --> BxMxN (0..255) + rtile = (255.0 * rtile).astype(numpy.uint8).reshape(rtile.shape[2], rtile.shape[0], rtile.shape[1]) + else: + rtile = apply_colormap(rtile, color_map) + options = img_profiles.get(driver, {}) return HttpResponse( - array_to_image(rtile, rmask, img_format=driver, color_map=color_map, **options), + array_to_image(rtile, rmask, img_format=driver, **options), content_type="image/{}".format(ext) ) \ No newline at end of file diff --git a/app/static/app/js/components/Map.jsx b/app/static/app/js/components/Map.jsx index 752ef2fe..dc206e8d 100644 --- a/app/static/app/js/components/Map.jsx +++ b/app/static/app/js/components/Map.jsx @@ -93,8 +93,10 @@ class Map extends React.Component { const { url, meta, type } = tile; let metaUrl = url + "metadata"; - if (type == "plant") metaUrl += "?expr=" + encodeURIComponent("(b2-b1)/(b2+b1-b3)") + "&rescale=0.02,0.1&color_map=rdylgn" + if (type == "plant") metaUrl += "?expr=" + encodeURIComponent("(b2-b1)/(b2+b1-b3)") + "&rescale=0.02,0.1&color_map=rdylgn"; + if (type == "dsm") metaUrl += "?rescale=140%2C170&hillshade=1"; + console.log(type, metaUrl); this.tileJsonRequests.push($.getJSON(metaUrl) .done(mres => { const { scheme, name, maxzoom } = mres; diff --git a/requirements.txt b/requirements.txt index 7522556a..7fd19620 100644 --- a/requirements.txt +++ b/requirements.txt @@ -55,4 +55,5 @@ vine==1.1.4 webcolors==1.5 rasterio==1.1.0 rio-tiler==1.3.0 -rio-color==1.0.0 \ No newline at end of file +rio-color==1.0.0 +mercantile==1.1.2 \ No newline at end of file