Merge pull request #522 from EOxServer/band-statistics

Band statistics and virtual assets

eoxserver/render/browse/functions.py

@@ -206,6 +206,64 @@ def pansharpen(pan_ds, *spectral_dss):
     return out_ds
 
 
+def percentile(ds, perc, default=0):
+    band = ds.GetRasterBand(1)
+    histogram = band.GetDefaultHistogram()
+    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)
+    return default
+
+
+def _has_stats(band):
+    return 'STATISTICS_MINIMUM' in band.GetMetadata()
+
+def statistics_min(ds, default=0):
+    band = ds.GetRasterBand(1)
+    if _has_stats(band):
+        min_, _, _, _ = band.GetStatistics(True, False)
+        return min_
+    return default
+
+
+def statistics_max(ds, default=0):
+    band = ds.GetRasterBand(1)
+    if _has_stats(band):
+        _, max_, _, _ = band.GetStatistics(True, False)
+        return max_
+    return default
+
+def statistics_mean(ds, default=0):
+    band = ds.GetRasterBand(1)
+    if _has_stats(band):
+        _, _, mean, _ = band.GetStatistics(True, False)
+        return mean
+    return default
+
+
+def statistics_stddev(ds, default=0):
+    band = ds.GetRasterBand(1)
+    if _has_stats(band):
+        _, _, _, stddev = band.GetStatistics(True, False)
+        return stddev
+    return default
+
+
+def interpolate(ds, x1, x2, y1, y2):
+    """Perform linear interpolation for x between (x1,y1) and (x2,y2) """
+    band = ds.GetRasterBand(1)
+    x = band.ReadAsArray()
+    x = ((y2 - y1) * x + x2 * y1 - x1 * y2) / (x2 - x1)
+    return gdal_array.OpenNumPyArray(x, True)
+
+
 def wrap_numpy_func(function):
     @wraps(function)
     def inner(ds, *args, **kwargs):
@@ -252,6 +310,12 @@ def inner(ds, *args, **kwargs):
     'roughness': roughness,
     'contours': contours,
     'pansharpen': pansharpen,
+    'percentile': percentile,
+    'statistics_min': statistics_min,
+    'statistics_max': statistics_max,
+    'statistics_mean': statistics_mean,
+    'statistics_stddev': statistics_stddev,
+    'interpolate': interpolate,
 }
 
 

eoxserver/render/browse/util.py

@@ -20,8 +20,8 @@ def warp_fields(coverages, field_name, bbox, crs, width, height):
         field.data_type
     )
     nil_value = None
+    band = out_ds.GetRasterBand(1)
     if field.nil_values:
-        band = out_ds.GetRasterBand(1)
         nil_value = float(field.nil_values[0][0])
         band.SetNoDataValue(nil_value)
         band.Fill(nil_value)
@@ -40,6 +40,23 @@ def warp_fields(coverages, field_name, bbox, crs, width, height):
 
     out_ds.SetProjection(sr.ExportToWkt())
 
+    # set initial statistics
+    stats = coverages[0].get_statistics_for_field(field_name)
+    if stats:
+        band.SetStatistics(
+            stats.minimum or 0,
+            stats.maximum or 0,
+            stats.mean or 0,
+            stats.stddev or 0,
+        )
+        histogram = stats.histogram
+        if histogram:
+            band.SetDefaultHistogram(
+                histogram.min or 0,
+                histogram.max or 0,
+                histogram.buckets,
+            )
+
     for coverage in coverages:
         location = coverage.get_location_for_field(field_name)
         band_index = coverage.get_band_index_for_field(field_name)

eoxserver/render/coverage/objects.py

@@ -451,11 +451,29 @@ def format(self):
         return self._format
 
 
+class Histogram(object):
+    def __init__(self, min, max, buckets):
+        self.min = min
+        self.max = max
+        self.buckets = buckets
+
+
+class Statistics(object):
+    def __init__(self, mean, minimum, maximum, stddev, valid_percent, histogram):
+        self.mean = mean
+        self.minimum = minimum
+        self.maximum = maximum
+        self.stddev = stddev
+        self.valid_percent = valid_percent
+        self.histogram = histogram
+
+
 class ArraydataLocation(Location):
-    def __init__(self, path, env, format, start_field, end_field):
+    def __init__(self, path, env, format, start_field, end_field, band_statistics):
         super(ArraydataLocation, self).__init__(path, env, format)
         self._start_field = start_field
         self._end_field = end_field
+        self._band_statistics = band_statistics
 
     @property
     def start_field(self):
@@ -472,6 +490,10 @@ def field_count(self):
     def field_index_to_band_index(self, field_index):
         return field_index - self.start_field
 
+    def field_statistics(self, field_index):
+        band_index = self.field_index_to_band_index(field_index)
+        return self._band_statistics[band_index]
+
 
 class Coverage(object):
     """ Representation of a coverage for internal processing.
@@ -574,46 +596,43 @@ def extent(self):
         elif self.footprint:
             return self.footprint.extent
 
-    def get_location_for_field(self, field_or_identifier):
+    def lookup_field(self, field_or_identifier):
         if isinstance(field_or_identifier, Field):
             field = field_or_identifier
             if field not in self.range_type:
                 return None
+            return field
         else:
             try:
-                field = next(
+                return next(
                     field
                     for field in self.range_type
                     if field.identifier == field_or_identifier
                 )
             except StopIteration:
                 return None
 
+    def get_location_for_field(self, field_or_identifier):
+        field = self.lookup_field(field_or_identifier)
+
         index = field.index
         for location in self.arraydata_locations:
             if index >= location.start_field and index <= location.end_field:
                 return location
 
     def get_band_index_for_field(self, field_or_identifier):
-        if isinstance(field_or_identifier, Field):
-            field = field_or_identifier
-            if field not in self.range_type:
-                return None
-        else:
-            try:
-                field = next(
-                    field
-                    for field in self.range_type
-                    if field.identifier == field_or_identifier
-                )
-            except StopIteration:
-                return None
+        field = self.lookup_field(field_or_identifier)
 
         index = field.index
         for location in self.arraydata_locations:
             if index >= location.start_field and index <= location.end_field:
                 return index - location.start_field + 1
 
+    def get_statistics_for_field(self, field_or_identifier):
+        field = self.lookup_field(field_or_identifier)
+        location = self.get_location_for_field(field)
+        return location.field_statistics(field.index)
+
     @classmethod
     def from_model(cls, model):
         # use coverages EO metadata by default and fill up with
@@ -629,13 +648,30 @@ def from_model(cls, model):
             footprint = model.parent_product.footprint
         eo_metadata = EOMetadata(begin_time, end_time, footprint)
 
-        arraydata_locations = [
-            ArraydataLocation(
-                get_vsi_path(item), get_vsi_env(item.storage), item.format,
-                item.field_index, item.field_index + (item.band_count - 1)
+        arraydata_locations = []
+        for item in model.arraydata_items.all():
+            statistics = [None] * item.band_count
+            for band_statistics in item.array_statistics.all():
+                statistics[band_statistics.band_index - 1] = Statistics(
+                    band_statistics.mean,
+                    band_statistics.minimum,
+                    band_statistics.maximum,
+                    band_statistics.stddev,
+                    band_statistics.valid_percent,
+                    Histogram(
+                        band_statistics.histogram.get('min', 0),
+                        band_statistics.histogram.get('max', 0),
+                        band_statistics.histogram.get('buckets', []),
+                    ),
+                )
+
+            arraydata_locations.append(
+                ArraydataLocation(
+                    get_vsi_path(item), get_vsi_env(item.storage), item.format,
+                    item.field_index, item.field_index + (item.band_count - 1),
+                    statistics
+                )
             )
-            for item in model.arraydata_items.all()
-        ]
 
         metadata_locations = [
             Location(

eoxserver/resources/coverages/management/commands/stac.py

@@ -68,6 +68,14 @@ def add_arguments(self, parser):
                 'Optional.'
             )
         )
+        register_parser.add_argument(
+            '--create-type', '-c', dest='create_type',
+            action="store_true", default=False,
+            help=(
+                'Whether to automatically create a product type from the STAC '
+                'Item. Optional.'
+            )
+        )
         register_parser.add_argument(
             "--replace", "-r",
             dest="replace", action="store_true", default=False,
@@ -113,7 +121,7 @@ def handle(self, subcommand, *args, **kwargs):
         elif subcommand == "types":
             self.handle_types(*args, **kwargs)
 
-    def handle_register(self, location, stdin, type_name, replace,
+    def handle_register(self, location, stdin, type_name, create_type, replace,
                         *args, **kw):
         if stdin:
             location = '.'
@@ -123,8 +131,15 @@ def handle_register(self, location, stdin, type_name, replace,
             with open(location) as f:
                 values = json.load(f)
 
+        if create_type:
+            product_type, is_new = \
+                create_product_type_from_stac_item(values, type_name)
+            self.print_msg("Created new product type %s" % product_type.name)
+            type_name = product_type.name
+
         product, replaced = register_stac_product(
-            location, values, type_name, replace=replace
+            values, type_name, replace=replace,
+            file_href=location if not stdin else None,
         )
         self.print_msg(
             "Successfully %s product %s" % (

eoxserver/resources/coverages/migrations/0011_bandstatistics.py

@@ -0,0 +1,32 @@
+# Generated by Django 2.2.17 on 2022-03-25 10:28
+
+from django.db import migrations, models
+import django.db.models.deletion
+import jsonfield.fields
+
+
+class Migration(migrations.Migration):
+
+    dependencies = [
+        ('coverages', '0010_polarisation_channels_expand'),
+    ]
+
+    operations = [
+        migrations.CreateModel(
+            name='BandStatistics',
+            fields=[
+                ('id', models.AutoField(auto_created=True, primary_key=True, serialize=False, verbose_name='ID')),
+                ('band_index', models.PositiveSmallIntegerField(default=1)),
+                ('mean', models.FloatField(blank=True, null=True)),
+                ('minimum', models.FloatField(blank=True, null=True)),
+                ('maximum', models.FloatField(blank=True, null=True)),
+                ('stddev', models.FloatField(blank=True, null=True)),
+                ('valid_percent', models.FloatField(blank=True, null=True)),
+                ('histogram', jsonfield.fields.JSONField(blank=True, null=True)),
+                ('arraydata_item', models.ForeignKey(on_delete=django.db.models.deletion.CASCADE, related_name='array_statistics', to='coverages.ArrayDataItem')),
+            ],
+            options={
+                'unique_together': {('arraydata_item', 'band_index')},
+            },
+        ),
+    ]

eoxserver/resources/coverages/models.py

@@ -44,6 +44,7 @@
 from django.utils.timezone import now
 from django.utils.encoding import python_2_unicode_compatible
 from model_utils.managers import InheritanceManager
+from jsonfield import JSONField
 
 from eoxserver.backends import models as backends
 from eoxserver.core.util.timetools import isoformat
@@ -435,7 +436,6 @@ class Meta:
         unique_together = [('eo_object', 'semantic')]
 
 
-
 class Browse(backends.DataItem):
     product = models.ForeignKey(Product, on_delete=models.CASCADE, related_name='browses', **mandatory)
     browse_type = models.ForeignKey(BrowseType, on_delete=models.CASCADE, **optional)
@@ -484,6 +484,21 @@ class Meta:
         unique_together = [('coverage', 'field_index')]
 
 
+class BandStatistics(models.Model):
+    arraydata_item = models.ForeignKey(ArrayDataItem, on_delete=models.CASCADE, related_name='array_statistics', **mandatory)
+    band_index = models.PositiveSmallIntegerField(default=1, **mandatory)
+
+    mean = models.FloatField(**optional)
+    minimum = models.FloatField(**optional)
+    maximum = models.FloatField(**optional)
+    stddev = models.FloatField(**optional)
+    valid_percent = models.FloatField(**optional)
+    histogram = JSONField(**optional)
+
+    class Meta:
+        unique_together = [('arraydata_item', 'band_index')]
+
+
 # ==============================================================================
 # Additional Metadata Models for Collections, Products and Coverages
 # ==============================================================================