From e89c838aa2a00b11e29ea9550990d795a2f674aa Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Sat, 11 Jan 2025 15:57:57 -0500 Subject: [PATCH 1/3] More consistent hillshading across zoom levels --- app/api/tiler.py | 2 +- app/raster_utils.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/app/api/tiler.py b/app/api/tiler.py index a5700772..f26eedda 100644 --- a/app/api/tiler.py +++ b/app/api/tiler.py @@ -463,7 +463,7 @@ class Tiles(TaskNestedView): if tile.data.shape[0] != 1: raise exceptions.ValidationError( _("Cannot compute hillshade of non-elevation raster (multiple bands found)")) - delta_scale = (maxzoom + ZOOM_EXTRA_LEVELS + 1 - z) * 4 + delta_scale = (maxzoom + ZOOM_EXTRA_LEVELS + 1 - z) ** 2 dx = src.dataset.meta["transform"][0] * delta_scale dy = -src.dataset.meta["transform"][4] * delta_scale ls = LightSource(azdeg=315, altdeg=45) diff --git a/app/raster_utils.py b/app/raster_utils.py index d4fee27b..ed8c6808 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -233,7 +233,7 @@ def export_raster(input, output, **opts): intensity = None if hillshade is not None and hillshade > 0: - delta_scale = (ZOOM_EXTRA_LEVELS + 1) * 4 + delta_scale = ZOOM_EXTRA_LEVELS ** 2 dx = src.meta["transform"][0] * delta_scale dy = -src.meta["transform"][4] * delta_scale ls = LightSource(azdeg=315, altdeg=45) From 35c191efe42694ec672ac541d7558213455313ad Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Sat, 11 Jan 2025 17:19:39 -0500 Subject: [PATCH 2/3] Avoid sign flip computation --- app/api/hillshade.py | 6 ------ app/api/tiler.py | 2 +- app/raster_utils.py | 2 +- 3 files changed, 2 insertions(+), 8 deletions(-) diff --git a/app/api/hillshade.py b/app/api/hillshade.py index 50006228..08840955 100644 --- a/app/api/hillshade.py +++ b/app/api/hillshade.py @@ -71,12 +71,6 @@ class LightSource: A 2d array of illumination values between 0-1, where 0 is completely in shadow and 1 is completely illuminated. """ - - # Because most image and raster GIS data has the first row in the array - # as the "top" of the image, dy is implicitly negative. This is - # consistent to what `imshow` assumes, as well. - dy = -dy - # compute the normal vectors from the partial derivatives e_dy, e_dx = np.gradient(vert_exag * elevation, dy, dx) diff --git a/app/api/tiler.py b/app/api/tiler.py index f26eedda..d04c7d1a 100644 --- a/app/api/tiler.py +++ b/app/api/tiler.py @@ -465,7 +465,7 @@ class Tiles(TaskNestedView): _("Cannot compute hillshade of non-elevation raster (multiple bands found)")) delta_scale = (maxzoom + ZOOM_EXTRA_LEVELS + 1 - z) ** 2 dx = src.dataset.meta["transform"][0] * delta_scale - dy = -src.dataset.meta["transform"][4] * delta_scale + dy = src.dataset.meta["transform"][4] * delta_scale ls = LightSource(azdeg=315, altdeg=45) # Remove elevation data from edge buffer tiles diff --git a/app/raster_utils.py b/app/raster_utils.py index ed8c6808..310766ee 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -235,7 +235,7 @@ def export_raster(input, output, **opts): if hillshade is not None and hillshade > 0: delta_scale = ZOOM_EXTRA_LEVELS ** 2 dx = src.meta["transform"][0] * delta_scale - dy = -src.meta["transform"][4] * delta_scale + dy = src.meta["transform"][4] * delta_scale ls = LightSource(azdeg=315, altdeg=45) intensity = ls.hillshade(arr[0], dx=dx, dy=dy, vert_exag=hillshade) intensity = intensity * 255.0 From c6225fff117fb949dd5105d5ffdadef2687de0af Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Sat, 11 Jan 2025 17:39:54 -0500 Subject: [PATCH 3/3] Disable hillshade contrast stretching --- app/api/hillshade.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/app/api/hillshade.py b/app/api/hillshade.py index 08840955..6bfab983 100644 --- a/app/api/hillshade.py +++ b/app/api/hillshade.py @@ -109,19 +109,19 @@ class LightSource: intensity = normals.dot(self.direction.astype(np.float32)) # Apply contrast stretch - imin, imax = intensity.min(), intensity.max() + # imin, imax = np.nanmin(intensity), np.nanmax(intensity) intensity *= fraction # Rescale to 0-1, keeping range before contrast stretch # If constant slope, keep relative scaling (i.e. flat should be 0.5, # fully occluded 0, etc.) - if (imax - imin) > 1e-6: + # if (imax - imin) > 1e-6: # Strictly speaking, this is incorrect. Negative values should be # clipped to 0 because they're fully occluded. However, rescaling # in this manner is consistent with the previous implementation and # visually appears better than a "hard" clip. - intensity -= imin - intensity /= (imax - imin) + # intensity -= imin + # intensity /= (imax - imin) intensity = np.clip(intensity, 0, 1) return intensity \ No newline at end of file