Merge branch 'master' of github.com:EOxServer/eoxserver

eoxserver/render/browse/functions.py

@@ -28,12 +28,14 @@
 import logging
 from uuid import uuid4
 from functools import wraps
+from typing import List, SupportsFloat as Numeric
 
 import numpy as np
 
 from eoxserver.contrib import gdal
 from eoxserver.contrib import ogr
 from eoxserver.contrib import gdal_array
+from eoxserver.render.browse.util import convert_dtype
 
 
 logger = logging.getLogger(__name__)
@@ -212,10 +214,6 @@ def percentile(ds, perc, default=0):
     if histogram:
         min_, max_, _, buckets = histogram
         bucket_diff = (max_ - min_) / len(buckets)
-        nodata = band.GetNoDataValue()
-        if nodata is not None:
-            # Set bucket of nodata value to 0
-            buckets[round((nodata - min_) / bucket_diff)] = 0
         cumsum = np.cumsum(buckets)
         bucket_index = np.searchsorted(cumsum, cumsum[-1] * (perc / 100))
         return min_ + (bucket_index * bucket_diff)
@@ -256,14 +254,45 @@ def statistics_stddev(ds, default=0):
     return default
 
 
-def interpolate(ds, x1, x2, y1, y2):
-    """Perform linear interpolation for x between (x1,y1) and (x2,y2) """
+def interpolate(
+        ds:gdal.Dataset, x1:Numeric, x2:Numeric, y1:Numeric, y2:Numeric, clip:bool=False, nodata_range:List[Numeric]=None
+    ):
+    """Perform linear interpolation for x between (x1,y1) and (x2,y2) with
+    optional clamp and additional masking out multiple no data value ranges
+
+    Args:
+        ds (gdal.Dataset): input gdal dataset
+        x1 (Numeric): linear interpolate from min
+        x2 (Numeric): linear interpolate from max
+        y1 (Numeric): linear interpolate to min
+        y2 (Numeric): linear interpolate to max
+        clip (bool, optional): if set to True, performs clip (values below y1 set to y1, values above y2 set to y2). Defaults to False.
+        additional_no_data (List, optional): additionally masks out (sets to band no_data_value) a range of values. Defaults to []. Example [1,5]
+
+    Returns:
+        gdal.Dataset: Interpolated dataset
+    """
     band = ds.GetRasterBand(1)
-    x = band.ReadAsArray()
-    # NOTE: this formula uses large numbers which lead to overflows on uint16
-    x = x.astype("int64")
-    x = ((y2 - y1) * x + x2 * y1 - x1 * y2) / (x2 - x1)
-    return gdal_array.OpenNumPyArray(x, True)
+    nodata_value = band.GetNoDataValue()
+    orig_image = band.ReadAsArray()
+    # NOTE: the interpolate formula uses large numbers which lead to overflows on uint16
+    if orig_image.dtype != convert_dtype(orig_image.dtype):
+        orig_image = orig_image.astype(convert_dtype(orig_image.dtype))
+    interpolated_image = ((y2 - y1) * orig_image + x2 * y1 - x1 * y2) / (x2 - x1)
+    if clip:
+        # clamp values below min to min and above max to max
+        np.clip(interpolated_image, y1, y2, out=interpolated_image)
+
+    if nodata_value is not None:
+        # restore nodata pixels on interpolated array from original array
+        interpolated_image[orig_image == nodata_value] = nodata_value
+        if nodata_range:
+            # apply mask of additional nodata ranges from original array on interpolated array
+            interpolated_image[(orig_image >= nodata_range[0]) & (orig_image <= nodata_range[1])] = nodata_value
+
+    ds = gdal_array.OpenNumPyArray(interpolated_image, True)
+    ds.GetRasterBand(1).SetNoDataValue(nodata_value)
+    return ds
 
 
 def wrap_numpy_func(function):

eoxserver/render/browse/generate.py

@@ -116,6 +116,7 @@ class BandExpressionError(ValueError):
     _ast.Add,
     _ast.Sub,
     _ast.Num if hasattr(_ast, 'Num') else _ast.Constant,
+    _ast.List,
 
     _ast.BitAnd,
     _ast.BitOr,
@@ -490,6 +491,12 @@ def _evaluate_expression(expr, fields_and_datasets, variables, cache):
     elif hasattr(_ast, 'Constant') and isinstance(expr, _ast.Constant):
         result = expr.value
 
+    elif hasattr(_ast, 'List') and isinstance(expr, _ast.List):
+        result = [
+            _evaluate_expression(
+                item, fields_and_datasets, variables, cache,
+            ) for item in expr.elts
+        ]
     else:
         raise BandExpressionError('Invalid expression node %s' % expr)
 

eoxserver/render/browse/util.py

@@ -1,8 +1,12 @@
 from uuid import uuid4
+import numpy as np
+import logging
 
 from eoxserver.contrib import gdal, osr
 from eoxserver.resources.coverages import crss
 
+logger = logging.getLogger(__name__)
+
 
 def create_mem_ds(width, height, data_type):
     driver = gdal.GetDriverByName('MEM')
@@ -78,3 +82,27 @@ def warp_fields(coverages, field_name, bbox, crs, width, height):
             gdal.Unlink(vrt_filename)
 
     return out_ds
+
+
+def convert_dtype(dtype:np.dtype):
+    """Maps numpy dtype to a larger itemsize
+    to avoid value overflow during mathematical operations
+
+    Args:
+        dtype (np.dtype): input dtype
+
+    Returns:
+        dtype (np.dtype): either one size larger dtype or original dtype
+    """
+    mapping = {
+        np.dtype(np.int8): np.dtype(np.int16),
+        np.dtype(np.int16): np.dtype(np.int32),
+        np.dtype(np.int32): np.dtype(np.int64),
+        np.dtype(np.uint8): np.dtype(np.int16),
+        np.dtype(np.uint16): np.dtype(np.int32),
+        np.dtype(np.uint32): np.dtype(np.int64),
+        np.dtype(np.float16): np.dtype(np.float32),
+        np.dtype(np.float32): np.dtype(np.float64),
+    }
+    output_dtype = mapping.get(dtype, dtype)
+    return output_dtype

eoxserver/render/mapserver/factories.py

@@ -414,20 +414,30 @@ def make_browse_layer_generator(self, map_obj, browses, map_,
                             range_ = _get_range(field)
 
                         for layer_obj in layer_objs:
-                            if browse.show_out_of_bounds_data and nodata_value == 0 and range_[0] >= 1:
-                                # NOTE: this trick only works with 0 as nodata value
-                                #       and if there actually is data below the min range
-
-                                # TODO: currently we only have data with exactly this case (integer values,
-                                #       0 as no data value). We can generalize this to float and for coverages
-                                #       as necessary.
-
-                                # only need 1:1 mapping if that's not the min
-                                lut_start= ("1:1," if range_[0] > 1 else "")
-                                lut = f"{lut_start}%d:1,%d:255" % tuple(range_)
+                            # NOTE: Only works if browsetype nodata is lower than browse_type_min by at least 1
+                            if browse.show_out_of_bounds_data:
+                                # final LUT for min,max 200,700 and nodata=0 should look like:
+                                # 0:0,1:1,200:1,700:256
+                                lut_inputs = {
+                                    range_[0]: 1,
+                                    range_[1]: 256,
+                                }
+                                if nodata_value is not None:
+                                    # no_data_value_plus +1 to ensure that only no_data_value is
+                                    # rendered as black (transparent)
+                                    nodata_value_plus = nodata_value + 1
+                                    lut_inputs[nodata_value] = 0
+                                    lut_inputs[nodata_value_plus] = 1
+
+                                # LUT inputs needs to be ascending
+                                sorted_inputs = {
+                                    k: v for k, v in sorted(list(lut_inputs.items()))
+                                }
+                                lut = ",".join("%d:%d" % (k,v) for k,v in sorted_inputs.items())
 
                                 layer_obj.setProcessingKey("LUT_%d" % i, lut)
                             else:
+                                # due to offsite 0,0,0 will make all pixels below or equal to min transparent
                                 layer_obj.setProcessingKey(
                                     "SCALE_%d" % i,
                                     "%s,%s" % tuple(range_)