feat(heatmap): add heatmap sub-layer for collection layers

Heatmap layers render the footprints of the products within a collection in a heatmap: altering the color depending on how many footprints intersect at a given pixel. The "styles" paramater alters the colorscale applied, whereas the "dim_range" changes the data range for the colorscale.

Updated documentation for heatmap layers.
Added simple test for heatmap layer.

autotest/autotest_services/tests/wms/base.py

@@ -97,6 +97,7 @@ class WMS13GetMapTestCase(testbase.RasterTestCase):
     frmt = "image/jpeg"
     time = None
     dim_bands = None
+    dim_range = None
 
     swap_axes = True
 
@@ -128,6 +129,9 @@ def getRequest(self):
         if self.dim_bands:
             params += "&dim_bands=%s" % self.dim_bands
 
+        if self.dim_range:
+            params += "&dim_range=%s" % self.dim_range
+
         if self.httpHeaders is None:
             return (params, "kvp")
         else:

autotest/autotest_services/tests/wms/test_v13.py

@@ -357,6 +357,20 @@ class WMS13GetMapDatasetStyledTestCase(wmsbase.WMS13GetMapTestCase):
     styles = ("color",)
     frmt = "image/png"
 
+
+# ==============================================================================
+# Heatmap
+# ==============================================================================
+
+class WMS13GetMapCollectionHeatmapTestCase(wmsbase.WMS13GetMapTestCase):
+    layers = ["MER_FRS_1P_reduced__heatmap"]
+    height = 50
+    width = 100
+    bbox = [-3.75, 32.158895, 28.326165, 46.3]
+    dim_range = "0 5"
+    frmt = "image/png"
+
+
 #===============================================================================
 # Feature Info
 #===============================================================================

documentation/users/instance.rst

@@ -152,8 +152,19 @@ EOXS_MAPSERVER_LAYER_FACTORIES
           'eoxserver.render.mapserver.factories.MaskLayerFactory',
           'eoxserver.render.mapserver.factories.MaskedBrowseLayerFactory',
           'eoxserver.render.mapserver.factories.OutlinesLayerFactory',
+          'eoxserver.render.mapserver.factories.HeatmapLayerFactory',
       ]
 
+
+DEFAULT_EOXS_MAPSERVER_HEATMAP_RANGE_DEFAULT = (0, 10)
+  The default range for heatmap layers when none are provided via ``dim_range``.
+
+  Default:
+
+  .. code-block:: python
+
+      (0, 10)
+
 EOXS_COVERAGE_METADATA_FORMAT_READERS
   The list of coverage metadata readers that will be employed to read metadata
   when a new coverage is registered.

documentation/users/operations/management.rst

@@ -408,6 +408,7 @@ This results in a catalog of the following available layers:
    geometries.
  - ``Collection__outlined``: this is a combination of the previous two layers:
    each Product is rendered in ``TRUE_COLOR`` with its outlines highlighted.
+ - ``Collection__heatmap``: this renders the heatmap of the Products footprints.
  - ``Collection__TRUE_COLOR``, ``Collection__FALSE_COLOR``,
    ``Collection__NDVI``: these are the browse visualizations with the
    definintions from earlier.
@@ -435,6 +436,8 @@ The following list shows all of these rendering options with an example product
     +-----------------------------------+---------------------------------------------------+
     | ``Collection__outlined``          | .. figure:: images/product_outlined.png           |
     +-----------------------------------+---------------------------------------------------+
+    | ``Collection__heatmap``           | .. figure:: images/product_heatmap.png            |
+    +-----------------------------------+---------------------------------------------------+
     | ``Collection__validity``          | .. figure:: images/product_validity.png           |
     +-----------------------------------+---------------------------------------------------+
     | ``Collection__masked_validity``   | .. figure:: images/product_masked_validity.png    |

documentation/users/services/wms.rst

@@ -122,7 +122,11 @@ parameters that are available with GetMap requests.
     |                           |     mask of the provided ``mask-name``.                   |                                  |                                |
     |                           |   - ``<browse-type-name>``: renders the product(s)        |                                  |                                |
     |                           |     according to the browse types instructions (or uses   |                                  |                                |
-    |                           |     an already existing browse if available.              |                                  |                                |
+    |                           |     an already existing browse if available)              |                                  |                                |
+    |                           |                                                           |                                  |                                |
+    |                           | - ``Collection``                                          |                                  |                                |
+    |                           |   - ``heatmap``: renders the contained products in a      |                                  |                                |
+    |                           |     heatmap.                                              |                                  |                                |
     +---------------------------+-----------------------------------------------------------+----------------------------------+--------------------------------+
     | styles                    | The style for each of the rendered layers to be           |                                  | M                              |
     |                           | rendered with. This must be either empty or a             |                                  |                                |

eoxserver/render/map/objects.py

@@ -26,7 +26,9 @@
 # ------------------------------------------------------------------------------
 
 from weakref import proxy
-from typing import List
+from typing import List, Optional, Tuple
+
+from django.contrib.gis.geos import GEOSGeometry
 
 from eoxserver.render.coverage.objects import (
     GRID_TYPE_TEMPORAL, GRID_TYPE_ELEVATION, Coverage, Mosaic,
@@ -293,6 +295,24 @@ def fill(self):
         return self._fill
 
 
+class HeatmapLayer(Layer):
+    """ Representation of a heatmap layer.
+    """
+    def __init__(self, name: str, style: str, footprints: List[GEOSGeometry],
+                 range: Optional[Tuple[float, float]] = None):
+        super(HeatmapLayer, self).__init__(name, style)
+        self._footprints = footprints
+        self._range = range
+
+    @property
+    def footprints(self) -> List[GEOSGeometry]:
+        return self._footprints
+
+    @property
+    def range(self) -> Optional[Tuple[float, float]]:
+        return self._range
+
+
 class Map(object):
     """ Abstract interpretation of a map to be drawn.
     """

eoxserver/render/mapserver/config.py

@@ -34,4 +34,11 @@
     'eoxserver.render.mapserver.factories.MaskLayerFactory',
     'eoxserver.render.mapserver.factories.MaskedBrowseLayerFactory',
     'eoxserver.render.mapserver.factories.OutlinesLayerFactory',
+    'eoxserver.render.mapserver.factories.HeatmapLayerFactory',
 ]
+
+
+# default for EOXS_MAPSERVER_HEATMAP_RANGE_DEFAULT: the default range for Heatmap
+# render requests
+
+DEFAULT_EOXS_MAPSERVER_HEATMAP_RANGE_DEFAULT = (0, 10)

eoxserver/render/mapserver/factories.py

@@ -34,11 +34,12 @@
     from itertools import zip_longest as izip_longest
 
 from django.conf import settings
+from django.contrib.gis.geos import GEOSGeometry
 from django.utils.module_loading import import_string
 
 from eoxserver.core.util.iteratortools import pairwise_iterative
 from eoxserver.contrib import mapserver as ms
-from eoxserver.contrib import vsi, vrt, gdal, osr
+from eoxserver.contrib import vsi, vrt, gdal, osr, ogr
 from eoxserver.render.browse.objects import (
     Browse, GeneratedBrowse, BROWSE_MODE_GRAYSCALE, BROWSE_MODE_RGBA
 )
@@ -47,13 +48,14 @@
 )
 from eoxserver.render.browse.defaultstyles import DEFAULT_RASTER_STYLES
 from eoxserver.render.map.objects import (
-    CoverageLayer, CoveragesLayer, MosaicLayer, OutlinedCoveragesLayer,
+    CoverageLayer, CoveragesLayer, HeatmapLayer, MosaicLayer, OutlinedCoveragesLayer,
     BrowseLayer, OutlinedBrowseLayer,
     MaskLayer, MaskedBrowseLayer, OutlinesLayer,
     Layer, Map,
 )
 from eoxserver.render.coverage.objects import Coverage, Field
 from eoxserver.render.mapserver.config import (
+    DEFAULT_EOXS_MAPSERVER_HEATMAP_RANGE_DEFAULT,
     DEFAULT_EOXS_MAPSERVER_LAYER_FACTORIES,
 )
 from eoxserver.render.colors import BASE_COLORS, COLOR_SCALES, OFFSITE_COLORS
@@ -658,6 +660,100 @@ def create(self, map_obj, layer):
         layer_obj.insertClass(class_obj)
 
 
+class HeatmapLayerFactory(BaseMapServerLayerFactory):
+    handled_layer_types = [HeatmapLayer]
+
+    def _create_vector_ds(self, footprints: List[GEOSGeometry]) -> gdal.Dataset:
+        """ Stores the given footprints in a GDAL vector dataset. It uses the in-memory
+            driver to reduce disk IO.
+        """
+        driver = gdal.GetDriverByName("Memory")
+        ds = driver.Create("", 0, 0, 0, gdal.GDT_Unknown)
+        layer = ds.CreateLayer("data")
+        from osgeo import osr
+        sr = osr.SpatialReference()
+        sr.ImportFromEPSG(4326)
+        for footprint in footprints:
+            feature = ogr.Feature(ogr.FeatureDefn())
+            geom = ogr.CreateGeometryFromWkt(footprint.wkt, sr)
+            feature.SetGeometryDirectly(geom)
+            layer.CreateFeature(feature)
+
+        return ds
+
+    def _rasterize_footprints(self, map_obj: ms.mapObj, vector_ds: gdal.Dataset,
+                              filename: str):
+        """ Rasterizes the footprints from a vector datasets into a raster dataset
+            of type Uint16 with the same dimension and spatial bounds as the provided map.
+            It uses additive mode in order to calculate how many geometries overlap with
+            a specific pixel.
+            The dataset is stored under a given filename (vsimem possible).
+        """
+        sr = osr.SpatialReference()
+        sr.ImportFromProj4(map_obj.getProjection())
+        extent = map_obj.extent
+
+        gdal.Rasterize(
+            filename,
+            vector_ds,
+            format="GTiff",
+            width=map_obj.width,
+            height=map_obj.height,
+            outputType=gdal.GDT_UInt16,
+            outputBounds=(extent.minx, extent.miny, extent.maxx, extent.maxy),
+            outputSRS=sr.ExportToWkt(),
+            initValues=[0],
+            burnValues=[1],
+            add=True
+        )
+
+    def create(self, map_obj: ms.mapObj, layer: Layer):
+        """_summary_
+
+        Args:
+            map_obj (ms.mapObj): _description_
+            layer (Layer): _description_
+
+        Returns:
+            _type_: _description_
+        """
+        assert isinstance(layer, HeatmapLayer)
+
+        vector_ds = self._create_vector_ds(layer.footprints)
+        filename_generator = FilenameGenerator('/vsimem/{uuid}.{extension}')
+        filename = filename_generator.generate("tif")
+        self._rasterize_footprints(map_obj, vector_ds, filename)
+
+        layer_obj = ms.layerObj(map_obj)
+        layer_obj.type = ms.MS_LAYER_RASTER
+        layer_obj.status = ms.MS_ON
+        layer_obj.data = filename
+
+        layer_obj.setProjection(map_obj.getProjection())
+
+        # TODO: make configuration possible for the default range here
+
+        default_range = getattr(
+            settings, 'EOXS_MAPSERVER_HEATMAP_RANGE_DEFAULT',
+            DEFAULT_EOXS_MAPSERVER_HEATMAP_RANGE_DEFAULT
+        )
+
+        range_ = layer.range or default_range
+        _create_raster_style(
+            DEFAULT_RASTER_STYLES[layer.style or "plasma"],
+            layer_obj,
+            range_[0],
+            range_[1],
+            [0],
+        )
+
+        return filename_generator
+
+    def destroy(self, map_obj, layer, filename_generator):
+        # cleanup temporary files
+        for filename in filename_generator.filenames:
+            vsi.unlink(filename)
+
 # ------------------------------------------------------------------------------
 # utils
 # ------------------------------------------------------------------------------
@@ -902,7 +998,7 @@ def _create_raster_style_ramp(raster_style, layer, minvalue=0, maxvalue=255,
         next_perc, next_color = next_item
 
         cls = ms.classObj()
-        cls.setExpression("([pixel] > %s AND [pixel] < %s)" % (
+        cls.setExpression("([pixel] > %s AND [pixel] <= %s)" % (
             (minvalue + prev_perc * interval),
             (minvalue + next_perc * interval)
         ))
@@ -923,7 +1019,7 @@ def _create_raster_style_ramp(raster_style, layer, minvalue=0, maxvalue=255,
         high_nil = min(high_nil_values)
         cls = ms.classObj()
         cls.setExpression(
-            "([pixel] > %s AND [pixel] < %s)" % (maxvalue, high_nil)
+            "([pixel] > %s AND [pixel] <= %s)" % (maxvalue, high_nil)
         )
         cls.group = name
         style = ms.styleObj()

eoxserver/services/ows/wms/layermapper.py

@@ -25,14 +25,14 @@
 # THE SOFTWARE.
 # ------------------------------------------------------------------------------
 
-from django.db.models import Case, Value, When, IntegerField, BooleanField
+from django.db.models import Case, Value, When, IntegerField
 from django.conf import settings
 
 from eoxserver.core.config import get_eoxserver_config
 from eoxserver.core.decoders import config, enum
 from eoxserver.core.util.timetools import isoformat
 from eoxserver.render.map.objects import (
-    CoverageLayer, CoveragesLayer, OutlinedCoveragesLayer, MosaicLayer,
+    CoverageLayer, CoveragesLayer, HeatmapLayer, OutlinedCoveragesLayer, MosaicLayer,
     OutlinesLayer, BrowseLayer, OutlinedBrowseLayer,
     MaskLayer, MaskedBrowseLayer,
     LayerDescription,
@@ -176,6 +176,17 @@ def get_layer_description(self, eo_object, default_raster_styles, default_geomet
                     )
                 )
 
+            if isinstance(eo_object, models.Collection):
+                sub_layers.append(
+                    LayerDescription(
+                        "%s%sheatmap" % (
+                            eo_object.identifier, self.suffix_separator,
+                        ),
+                        styles=default_raster_styles,
+                        dimensions=dimensions,
+                    )
+                )
+
             return LayerDescription(
                 name=eo_object.identifier,
                 bbox=eo_object.footprint.extent if eo_object.footprint else None,
@@ -424,6 +435,25 @@ def lookup_layer(self, layer_name, suffix, style, filters_expressions,
                     masked_browses=masked_browses
                 )
 
+            elif suffix == 'heatmap':
+                return HeatmapLayer(
+                    name=full_name,
+                    style=style,
+                    footprints=[
+                        # TODO: once testdata is available use products instead
+                        # of coverages
+                        product.footprint for product in self.iter_coverages(
+                            eo_object, filters_expressions, sort_by,
+                            #limit=limit_products
+                        )
+                        # product.footprint for product in self.iter_products(
+                        #     eo_object, filters_expressions, sort_by,
+                        #     limit=limit_products
+                        # )
+                    ],
+                    range=ranges[0] if ranges else None
+                )
+
             else:
                 # either browse type or mask type
                 browse_type = self.get_browse_type(eo_object, suffix)