Merge pull request #562 from EOxServer/CLM-cloud-coverage-process

CLM cloud coverage process

.github/workflows/ci.yml

@@ -9,6 +9,7 @@ on:
       - 'MANIFEST.in'
       - 'pyproject.toml'
       - 'eoxserver/**'
+      - 'autotest/**'
 jobs:
   build-docker:
     runs-on: ubuntu-latest

autotest/autotest/data/fixtures/clm_cloud_coverages.json

@@ -0,0 +1,75 @@
+[
+{
+    "model": "coverages.coveragetype",
+    "pk": 5,
+    "fields": {
+        "name": "CLM"
+    }
+},
+{
+    "model": "coverages.eoobject",
+    "pk": 22,
+    "fields": {
+        "identifier": "CLM_CLOUD_COVERAGE_PRODUCT",
+        "begin_time": "2020-02-15T07:59:36.409Z",
+        "end_time": "2020-02-20T08:00:36.342Z",
+        "footprint": "SRID=4326;POLYGON ((16.6276837970076 47.466520638378,16.7499627020568 47.4647191364739,16.7516507329145 47.5154960582897,16.6292540083727 47.5173007469159,16.6276837970076 47.466520638378))",
+        "inserted": "2019-01-01T00:00:00.000Z",
+        "updated": "2019-01-01T00:00:00.000Z"
+    }
+},
+{
+    "model": "coverages.product",
+    "pk": 22,
+    "fields": {
+        "product_type": null,
+        "package": null
+    }
+},
+{
+    "model": "coverages.eoobject",
+    "pk": 23,
+    "fields": {
+        "identifier": "CLM_CLOUD_COVERAGE",
+        "begin_time": null,
+        "end_time": null,
+        "footprint": "SRID=4326;POLYGON ((16.6276837970076 47.466520638378,16.7499627020568 47.4647191364739,16.7516507329145 47.5154960582897,16.6292540083727 47.5173007469159,16.6276837970076 47.466520638378))",
+        "inserted": "2019-01-01T00:00:00.000Z",
+        "updated": "2019-01-01T00:00:00.000Z"
+    }
+},
+{
+    "model": "coverages.coverage",
+    "pk": 23,
+    "fields": {
+        "grid": 1,
+        "axis_1_origin": null,
+        "axis_2_origin": null,
+        "axis_3_origin": null,
+        "axis_4_origin": null,
+        "axis_1_size": 569,
+        "axis_2_size": 486,
+        "axis_3_size": null,
+        "axis_4_size": null,
+        "coverage_type": 5,
+        "parent_product": 22,
+        "collections": [],
+        "mosaics": []
+    }
+},
+{
+    "model": "coverages.arraydataitem",
+    "pk": 21,
+    "fields": {
+        "storage": null,
+        "location": "autotest/data/CLM_clipped.tif",
+        "format": "image/tiff",
+        "coverage": 23,
+        "field_index": 0,
+        "band_count": 1,
+        "subdataset_type": null,
+        "subdataset_locator": null,
+        "bands_interpretation": 0
+    }
+}
+]

autotest/autotest/expected/WPS20ExecuteCloudCoverageOnCLM.json

@@ -0,0 +1 @@
+{"result": {"2020-02-15T07:59:36.409000+00:00": 0.06714258582744015}}

autotest/autotest/expected/WPS20ExecuteCloudCoverageOnCLMFewClouds.json

@@ -0,0 +1 @@
+{"result": {"2020-02-15T07:59:36.409000+00:00": 0.25697363393198314}}

autotest/autotest/expected/WPS20ExecuteCloudCoverageOnCLMMoreCloudy.json

@@ -0,0 +1 @@
+{"result": {"2020-02-15T07:59:36.409000+00:00": 0.0}}

autotest/autotest/expected/WPS20ExecuteCloudCoverageTinyGeometry.json

@@ -1 +1 @@
-{"result": {"2020-02-15T07:59:36.409000+00:00": 0}}
+{"result": {"2020-02-15T07:59:36.409000+00:00": 0.0}}

autotest/autotest_services/tests/wps/test_v20_cloud_coverage.py

@@ -32,10 +32,8 @@
 )
 
 
-class WPS20ExecuteCloudCoverageNonCloudy(
-    ContentTypeCheckMixIn, testbase.JSONTestCase
-):
-    fixtures = testbase.JSONTestCase.fixtures + ["cloud_coverages.json"]
+class WPS20ExecuteCloudCoverageNonCloudy(ContentTypeCheckMixIn, testbase.JSONTestCase):
+    fixtures = testbase.JSONTestCase.fixtures + ["scl_cloud_coverages.json"]
 
     expectedContentType = "application/json; charset=utf-8"
 
@@ -65,7 +63,7 @@ def getRequest(self):
 class WPS20ExecuteCloudCoverageCloudyGeometry(
     ContentTypeCheckMixIn, testbase.JSONTestCase
 ):
-    fixtures = testbase.JSONTestCase.fixtures + ["cloud_coverages.json"]
+    fixtures = testbase.JSONTestCase.fixtures + ["scl_cloud_coverages.json"]
 
     expectedContentType = "application/json; charset=utf-8"
 
@@ -95,7 +93,7 @@ def getRequest(self):
 class WPS20ExecuteCloudCoverageTinyGeometry(
     ContentTypeCheckMixIn, testbase.JSONTestCase
 ):
-    fixtures = testbase.JSONTestCase.fixtures + ["cloud_coverages.json"]
+    fixtures = testbase.JSONTestCase.fixtures + ["scl_cloud_coverages.json"]
 
     expectedContentType = "application/json; charset=utf-8"
 
@@ -120,3 +118,91 @@ def getRequest(self):
         </wps:Execute>
         """
         return (params, "xml")
+
+
+class WPS20ExecuteCloudCoverageOnCLM(ContentTypeCheckMixIn, testbase.JSONTestCase):
+    fixtures = testbase.JSONTestCase.fixtures + ["clm_cloud_coverages.json"]
+
+    expectedContentType = "application/json; charset=utf-8"
+
+    def getRequest(self):
+        params = """<wps:Execute
+        version="2.0.0"
+        service="WPS"
+        response="raw"
+        mode="sync"
+        xmlns:wps="http://www.opengis.net/wps/2.0"
+        xmlns:ows="http://www.opengis.net/ows/2.0" >
+          <ows:Identifier>CloudCoverage</ows:Identifier>
+          <wps:Input id="begin_time"><wps:Data>2020-01-01</wps:Data></wps:Input>
+          <wps:Input id="end_time"><wps:Data>2020-05-31</wps:Data></wps:Input>
+          <wps:Input id="geometry">
+            <wps:Data>
+              <wps:ComplexData mimeType="text/plain">POLYGON((16.69737831605056 47.47325091482521,16.711481970707467 47.50227740341751,16.665977994941493 47.4984370585647,16.65892616761306 47.47649970533886,16.69737831605056 47.47325091482521))</wps:ComplexData>
+            </wps:Data>
+          </wps:Input>
+          <wps:Output id="result" >
+          </wps:Output>
+        </wps:Execute>
+        """
+        return (params, "xml")
+
+
+class WPS20ExecuteCloudCoverageOnCLMMoreCloudy(
+    ContentTypeCheckMixIn, testbase.JSONTestCase
+):
+    fixtures = testbase.JSONTestCase.fixtures + ["clm_cloud_coverages.json"]
+
+    expectedContentType = "application/json; charset=utf-8"
+
+    def getRequest(self):
+        params = """<wps:Execute
+        version="2.0.0"
+        service="WPS"
+        response="raw"
+        mode="sync"
+        xmlns:wps="http://www.opengis.net/wps/2.0"
+        xmlns:ows="http://www.opengis.net/ows/2.0" >
+          <ows:Identifier>CloudCoverage</ows:Identifier>
+          <wps:Input id="begin_time"><wps:Data>2020-01-01</wps:Data></wps:Input>
+          <wps:Input id="end_time"><wps:Data>2020-05-31</wps:Data></wps:Input>
+          <wps:Input id="geometry">
+            <wps:Data>
+              <wps:ComplexData mimeType="text/plain">POLYGON((16.689481892710717 47.49088479184637,16.685847177764813 47.49102703651446,16.685862090779757 47.49375416408526,16.689267315989525 47.49363955046273,16.689481892710717 47.49088479184637))</wps:ComplexData>
+            </wps:Data>
+          </wps:Input>
+          <wps:Output id="result" >
+          </wps:Output>
+        </wps:Execute>
+        """
+        return (params, "xml")
+
+
+class WPS20ExecuteCloudCoverageOnCLMFewClouds(
+    ContentTypeCheckMixIn, testbase.JSONTestCase
+):
+    fixtures = testbase.JSONTestCase.fixtures + ["clm_cloud_coverages.json"]
+
+    expectedContentType = "application/json; charset=utf-8"
+
+    def getRequest(self):
+        params = """<wps:Execute
+        version="2.0.0"
+        service="WPS"
+        response="raw"
+        mode="sync"
+        xmlns:wps="http://www.opengis.net/wps/2.0"
+        xmlns:ows="http://www.opengis.net/ows/2.0" >
+          <ows:Identifier>CloudCoverage</ows:Identifier>
+          <wps:Input id="begin_time"><wps:Data>2020-01-01</wps:Data></wps:Input>
+          <wps:Input id="end_time"><wps:Data>2020-05-31</wps:Data></wps:Input>
+          <wps:Input id="geometry">
+            <wps:Data>
+              <wps:ComplexData mimeType="text/plain">POLYGON((16.689481892710717 47.49088479184637,16.685862090779757 47.49375416408526,16.717934765940697 47.50045333252222,16.689481892710717 47.49088479184637))</wps:ComplexData>
+            </wps:Data>
+          </wps:Input>
+          <wps:Output id="result" >
+          </wps:Output>
+        </wps:Execute>
+        """
+        return (params, "xml")

eoxserver/services/ows/wps/processes/get_cloud_coverage.py

@@ -30,13 +30,15 @@
 import functools
 from datetime import datetime
 from uuid import uuid4
+from typing import List, Callable, Optional
 
 from osgeo import ogr, osr
 
 from eoxserver.core import Component
 from eoxserver.contrib import gdal
 from eoxserver.resources.coverages import models
 from eoxserver.backends.access import gdal_open
+from eoxserver.services.ows.wps.exceptions import InvalidInputValueError
 from eoxserver.services.ows.wps.parameters import (
     LiteralData,
     ComplexData,
@@ -46,11 +48,11 @@
 )
 import logging
 
+
 logger = logging.getLogger(__name__)
 
 
 class CloudCoverageProcess(Component):
-
     identifier = "CloudCoverage"
     title = "Cloud coverage information about images of an AOI/TOI"
     description = ""
@@ -90,6 +92,18 @@ class CloudCoverageProcess(Component):
     SCL_LAYER_THIN_CIRRUS = 10
     SCL_LAYER_SATURATED_OR_DEFECTIVE = 1
 
+    # https://labo.obs-mip.fr/multitemp/sentinel-2/majas-native-sentinel-2-format/#English
+    # anything nonzero should be cloud, however that includes also cloud shadows which
+    # have a lot of false positives (or shadows that are not visible to the naked eye)
+    # for now only use the upper 4 bits which are:
+    # bit 4 (16) : clouds detected via mono-temporal thresholds
+    # bit 5 (32) : clouds detected via multi-temporal thresholds
+    # bit 6 (64) : thinnest clouds
+    # bit 7 (128) : high clouds detected by 1.38 µm
+    # sometimes bit 4 counts also seems to count things as cloud which don't appear to
+    # be clouds
+    CLM_MASK_ONLY_CLOUD = 0b11110000
+
     @staticmethod
     def execute(
         begin_time,
@@ -100,20 +114,36 @@ def execute(
         wkt_geometry = geometry[0].text
 
         # TODO Use queue object for more complex query if parent_product__footprint is not enough
-        coverages = models.Coverage.objects.filter(
+        relevant_coverages = models.Coverage.objects.filter(
             parent_product__begin_time__lte=end_time,
             parent_product__end_time__gte=begin_time,
             parent_product__footprint__intersects=wkt_geometry,
-            coverage_type__name="SCL",
         ).order_by("parent_product__begin_time")
 
-        logger.info("Matched %s coverages for cloud coverage", coverages.count())
+        if coverages_clm := relevant_coverages.filter(coverage_type__name="CLM"):
+            logger.info("Matched %s CLM covs for cloud coverage", coverages_clm.count())
+            calculation_fun = cloud_coverage_ratio_for_CLM
+            coverages = coverages_clm
+            # CLM is a bitmask, this value would mean that all types of cloud were found
+            # hopefully this never occurs naturally, so we can use it as no_data
+            no_data_value = 0b11111111
+
+        elif coverages_scl := relevant_coverages.filter(coverage_type__name="SCL"):
+            logger.info("Matched %s SCL covs for cloud coverage", coverages_scl.count())
+            calculation_fun = cloud_coverage_ratio_for_SCL
+            coverages = coverages_scl
+            no_data_value = None
+
+        else:
+            raise InvalidInputValueError("No coverage data found")
 
         with concurrent.futures.ThreadPoolExecutor(max_workers=5) as e:
             cloud_coverage_ratios = e.map(
                 functools.partial(
                     cloud_coverage_ratio_in_geometry,
+                    calculation_fun=calculation_fun,
                     wkt_geometry=wkt_geometry,
+                    no_data_value=no_data_value,
                 ),
                 [coverage.arraydata_items.get() for coverage in coverages],
             )
@@ -134,7 +164,51 @@ def execute(
 def cloud_coverage_ratio_in_geometry(
     data_item: models.ArrayDataItem,
     wkt_geometry: str,
+    calculation_fun: Callable[[List[int]], float],
+    no_data_value: Optional[int],
 ) -> float:
+    histogram = _histogram_in_geometry(
+        data_item=data_item,
+        wkt_geometry=wkt_geometry,
+        no_data_value=no_data_value,
+    )
+    return calculation_fun(histogram)
+
+
+def cloud_coverage_ratio_for_CLM(histogram: List[int]) -> float:
+    num_is_cloud = sum(
+        value
+        for index, value in enumerate(histogram)
+        if index & CloudCoverageProcess.CLM_MASK_ONLY_CLOUD > 0
+    )
+
+    num_pixels = sum(histogram)
+    return ((num_is_cloud / num_pixels)) if num_pixels != 0 else 0.0
+
+
+def cloud_coverage_ratio_for_SCL(histogram: List[int]) -> float:
+    num_cloud = sum(
+        histogram[scl_value]
+        for scl_value in [
+            CloudCoverageProcess.SCL_LAYER_CLOUD_MEDIUM_PROBABILITY,
+            CloudCoverageProcess.SCL_LAYER_CLOUD_HIGH_PROBABILITY,
+            CloudCoverageProcess.SCL_LAYER_THIN_CIRRUS,
+            CloudCoverageProcess.SCL_LAYER_SATURATED_OR_DEFECTIVE,
+        ]
+    )
+
+    num_no_data = histogram[CloudCoverageProcess.SCL_LAYER_NO_DATA]
+
+    num_pixels = sum(histogram) - num_no_data
+
+    return num_cloud / num_pixels if num_pixels != 0 else 0.0
+
+
+def _histogram_in_geometry(
+    data_item: models.ArrayDataItem,
+    wkt_geometry: str,
+    no_data_value: Optional[int],
+) -> List[int]:
     # NOTE: this is executed in threads, but all gdal operations are contained
     #       in here, so each thread has separate gdal data
 
@@ -151,6 +225,7 @@ def cloud_coverage_ratio_in_geometry(
                 cutlineDSName=geometry_mem_path,
                 cropToCutline=True,
                 warpOptions=["CUTLINE_ALL_TOUCHED=TRUE"],
+                dstNodata=no_data_value,
             ),
         )
 
@@ -161,25 +236,9 @@ def cloud_coverage_ratio_in_geometry(
         include_out_of_range=True,
     )
 
-    num_cloud = sum(
-        histogram[scl_value]
-        for scl_value in [
-            CloudCoverageProcess.SCL_LAYER_CLOUD_MEDIUM_PROBABILITY,
-            CloudCoverageProcess.SCL_LAYER_CLOUD_HIGH_PROBABILITY,
-            CloudCoverageProcess.SCL_LAYER_THIN_CIRRUS,
-            CloudCoverageProcess.SCL_LAYER_SATURATED_OR_DEFECTIVE,
-        ]
-    )
-
-    num_no_data = histogram[CloudCoverageProcess.SCL_LAYER_NO_DATA]
-
-    num_pixels = sum(histogram) - num_no_data
-
-    cloud_coverage_ratio = num_cloud / num_pixels if num_pixels != 0 else 0
-
     gdal.Unlink(tmp_ds)
 
-    return cloud_coverage_ratio
+    return histogram
 
 
 @contextlib.contextmanager