Merging dev-20240917 to main (#98)

* fixed the wrong collection Name for LST_reference_data * fix for #95 Solved with sorting the unique_filename_day list, fix for #95 * small changes to util_moveasset.py * added illumination angle "Fix" bug #97 - updated the 'terrainShadowMask' band to now contain the values from 0-90° representing the illumination angle based on the solar position and the SRTM 30m DEM. Values above 90° have been clamped to 90. - the functions addTerrainShadow and addTerrainShadow_predefined have been adapted to add the mask (calculated from SwissALTI3D or SwissSURFACE3D) as the value 100 to the terrainShadowMask band which now contains the illumination angle. - adapted the function addMaskedPixelCount to newly set the threshold for the terrainShadowMask to greater than 99. * include bad illumination for the terrainShadowMask * Replaced SRTM DEM with the Copernicus 30m DEM updated paths for slope and aspect to the new assets in projects/satromo-prod/assets/res * back to SRTM for illumination angle * Update prod_config.py --------- Co-authored-by: David Oesch <[email protected]>
swisstopo · Oct 14, 2024 · 358ef4e · 358ef4e
1 parent f85f7d0
commit 358ef4e
Show file tree

Hide file tree

Showing 5 changed files with 65 additions and 20 deletions.
diff --git a/configuration/prod_config.py b/configuration/prod_config.py
@@ -30,6 +30,7 @@
 # GDRIVE_MOUNT_DEV = r'M:\\satromo_export' # for GCS
 # under Windows, add \\ to escape the backslash like r'X:\\'
 S3_DESTINATION_DEV = r'X:\\'
+
 #  GITHUB
 GDRIVE_SOURCE_INT = "geedrivePROD:"
 GDRIVE_MOUNT_INT = "localgdrive"

diff --git a/main_functions/util_moveassets.py b/main_functions/util_moveassets.py
@@ -67,19 +67,16 @@ def initialize_gee_and_drive():
 # Remove whitespace and newline characters from the asset names
 asset_names = [name.strip() for name in asset_names]
 
-# Define the Earth Engine username for the destination account
-destination_username = "satromo-prod"
-
 # Copy each asset from the source to the destination
 for asset_name in asset_names:
     # Construct the source and destination asset paths
     # source_asset = "users/{}/SATROMO/{}".format(source_username, asset_name)
-    source_asset = "projects/satromo-int/assets/1991-2020_NDVI_STATS15/{}".format(
+    source_asset = "projects/satromo-int/assets/res/{}".format(
         asset_name)
     # destination_asset = 'projects/{}/assets/res/{}'.format(
     #     destination_username, asset_name)
-    destination_asset = 'projects/{}/assets/col/1991-2020_NDVI_SWISS/{}'.format(
-        destination_username, asset_name)
+    destination_asset = 'projects/satromo-prod/assets/res/{}'.format(
+        asset_name)
 
     # Copy the asset from the source to the destination
     ee.data.copyAsset(source_asset, destination_asset)

diff --git a/satromo_publish.py b/satromo_publish.py
@@ -717,6 +717,7 @@ def check_substrings_presence(file_merged, substring_to_check, additional_substr
     else:
         return False
 
+
 def check_asset_size(filename):
     """
     Checks the asset size of the product defined in the configuration 
@@ -973,10 +974,10 @@ def check_asset_size(filename):
                         # Authenticate with GDRIVE
                         initialize_drive()
 
-                        # os.remove(file_merged)
+                        # os.remove(file_merged
                         clean_up_gdrive(filename)
 
-                        # Remove each filename from the original list
+                        # Remove each filename from the original group and list
                         unique_filenames.remove(filename)
 
                 else:

diff --git a/step0_processors/step0_processor_s2_sr.py b/step0_processors/step0_processor_s2_sr.py
@@ -1,4 +1,5 @@
 import ee
+import numpy as np
 from main_functions import main_utils
 from .step0_utils import write_asset_as_empty
 
@@ -91,9 +92,15 @@ def generate_s2_sr_mosaic_for_single_date(day_to_process: str, collection: str,
     DEM_sa3d = ee.Image(
         "projects/satromo-prod/assets/res/SS3DR_SA3DRegio_10m_20kmBuffer_epsg32632")
 
+    # SRTM 30 - digital elevation model (slope and aspect) used for the atmospheric correction in sen2cor in a 30 m resolution
+    # source: https://developers.google.com/earth-engine/datasets/catalog/USGS_SRTMGL1_003
+    # processing: ee.Terrain.slope(DEM) and ee.Terrain.aspect(DEM) converted to radians
+    slope = ee.Image('projects/satromo-prod/assets/res/SRTM30m_slope_radians_epsg32632')
+    aspect = ee.Image('projects/satromo-prod/assets/res/SRTM30m_aspect_radians_epsg32632')
+
     # Terrain - very precise digital surface  model in a 10 m resolution
     # source: https://code.earthengine.google.com/ccfa64fe9827c93e2986e693983332e2
-    # processing:The shadow masks are  combined into a single image with multiple bands as asset per DOY.
+    # processing: The shadow masks are  combined into a single image with multiple bands as asset per DOY.
     terrain_shadow_collection = "projects/satromo-prod/assets/col/TERRAINSHADOW_SWISS/"
 
     # DX DY - Precalculated DX DY shifts
@@ -373,7 +380,39 @@ def maskCloudsAndShadowsSTwoCloudless(image):
             'cloud_mask_threshold': CLOUD_THRESHOLD         # threshold for cloud mask
         })
 
-    # This function detects and adds terrain shadows
+    # This function calculates and adds the illumination angle
+    def addIlluminationAngel(image):
+        # get the solar position
+        meanAzimuth = image.get('MEAN_SOLAR_AZIMUTH_ANGLE')
+        meanZenith = image.get('MEAN_SOLAR_ZENITH_ANGLE')
+
+        # Create an empty image to apply the expression
+        empty_image = ee.Image().float()
+
+        # Calculate illumination angle
+        illumination_cos = empty_image.expression(
+            'cos(sz) * cos(ps) + sin(sz) * sin(ps) * cos(sa - pa)',
+            {
+                'sz': ee.Number(meanZenith).multiply(np.pi).divide(180),  # Convert solar zenith to radians
+                'sa': ee.Number(meanAzimuth).multiply(np.pi).divide(180),  # Convert solar azimuth to radians
+                'ps': slope,
+                'pa': aspect
+            }
+        )
+        # The result is the cosine of the illumination angle
+        # To get the angle itself -> acos
+        illumination_angle_r = illumination_cos.acos()
+        illumination_angle = illumination_angle_r.multiply(180).divide(np.pi)
+
+        # Round to full numbers, convert to int, and cap at 90
+        illumination_angle = illumination_angle.round().toInt().clamp(0, 90).rename('terrainShadowMask')
+
+        # add the additonal terrainShadow band
+        image = image.addBands(illumination_angle)
+
+        return image
+
+    # This function detects and updates terrain shadows
     def addTerrainShadow(image):
         # get the solar position
         meanAzimuth = image.get('MEAN_SOLAR_AZIMUTH_ANGLE')
@@ -382,15 +421,17 @@ def addTerrainShadow(image):
         # Terrain shadow
         terrainShadow = ee.Terrain.hillShadow(
             DEM_sa3d, meanAzimuth, meanZenith, 100, True)
-        terrainShadow = terrainShadow.Not().rename(
-            'terrainShadowMask')  # invert the binaries
+        terrainShadow = terrainShadow.Not() # invert the binaries
 
-        # add the additonal terrainShadow band
-        image = image.addBands(terrainShadow)
+        # Update the existing terrainShadowMask band
+        updatedMask = image.select('terrainShadowMask').where(terrainShadow, 100)
+
+        # Replace the existing terrainShadowMask band
+        image = image.addBands(updatedMask, ['terrainShadowMask'], True)
 
         return image
 
-    # This  adds terrain shadows from precalcuated terrain
+    # This updates terrain shadows from precalcuated terrain
     def addTerrainShadow_predefined(image, start_date, terrain_shadow_collection, S2_sr):
 
         # Define the day of year
@@ -433,18 +474,20 @@ def find_closest_band(current, previous):
 
         band_image = terrain_shadow_asset.select(
             'shadow_' + closest_band_name.getInfo())
-        band_image = band_image.rename('terrainShadowMask')
+
+        # Update the existing terrainShadowMask band
+        updatedMask = image.select('terrainShadowMask').where(band_image, 100)
 
-        # Add the terrain shadow mask band to the input image
-        image = image.addBands(band_image)
+        # Replace the existing terrainShadowMask band
+        image = image.addBands(updatedMask, ['terrainShadowMask'], True)
 
         return image
 
     # This function adds the masked-pixel-percentage (clouds, cloud shadows, QA masks) as a property to each image
     def addMaskedPixelCount(image):
         # counter the umber of pixel that are masked by cloud or shadows
         image_mask = image.select('cloudAndCloudShadowMask').gt(
-            0).Or(image.select('terrainShadowMask').gt(0))
+            0).Or(image.select('terrainShadowMask').gt(99))
         statsMasked = image_mask.reduceRegion(
             reducer=ee.Reducer.sum(),
             geometry=image.geometry().intersection(aoi_CH_simplified),
@@ -519,6 +562,9 @@ def set_date(img):
             S2_sr = ee.ImageCollection(
                 S2_sr).map(maskCloudsAndShadowsSTwoCloudless)
 
+    # Add the illumination angle as terrainShadowMask band
+    S2_sr = S2_sr.map(addIlluminationAngel)
+
     # SWITCH
     if terrainShadowDetection is True:
         print('--- Terrain shadow detection applied ---')

diff --git a/step1_processors/step1_processor_vhi.py b/step1_processors/step1_processor_vhi.py
@@ -73,7 +73,7 @@ def loadNdviCurrentData(image):
     """
     # Apply the cloud and terrain shadow mask within the S2 image collection
     def applyMasks(image):
-        image = image.updateMask(image.select('terrainShadowMask').eq(0))
+        image = image.updateMask(image.select('terrainShadowMask').lt(65))
         image = image.updateMask(image.select('cloudAndCloudShadowMask').eq(0))
         return image
     S2_col_masked = image.map(applyMasks)