More consistent hillshading across zoom levels

app/api/hillshade.py

@@ -71,12 +71,6 @@ def hillshade(self, elevation, vert_exag=1, dx=1, dy=1, fraction=1.):
             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)
 
@@ -115,19 +109,19 @@ def shade_normals(self, normals, fraction=1.):
         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

app/api/tiler.py

@@ -463,9 +463,9 @@ def get(self, request, pk=None, project_pk=None, tile_type="", z="", x="", y="",
                 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
+                dy = src.dataset.meta["transform"][4] * delta_scale
                 ls = LightSource(azdeg=315, altdeg=45)
 
                 # Remove elevation data from edge buffer tiles

app/raster_utils.py

@@ -233,9 +233,9 @@ def write_band(arr, dst, band_num):
 
                 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
+                    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