Project Canopy Proposed Changes (v1)

cloud_score.py

@@ -1,6 +1,6 @@
 import ee
 
-from helpers import rescale, dilatedErossion
+from .helpers import rescale, dilatedErossion
 
 def sentinel2CloudScore(img):
     toa = img.select(['B1','B2','B3','B4','B5','B6','B7','B8','B8A', 'B9','B10', 'B11','B12']) \
@@ -34,11 +34,18 @@ def sentinel2CloudScore(img):
     score = score.max(ee.Image(0.001))
 
     # Remove small regions and clip the upper bound
-    dilated = dilatedErossion(score).min(ee.Image(1.0))
+    ## commenting this out--pointless to do this unless we do the
+    ## following line (42), which is also commented out in the JS
+    ## dilated = dilatedErossion(score).min(ee.Image(1.0))
+
+    ## This was commented out in the javascript--why?
+    ## score = score.multiply(dilated)
 
     score = score.reduceNeighborhood(
         reducer=ee.Reducer.mean(),
         kernel=ee.Kernel.square(5)
       )
+
+    print("Generated cloud score")
 
     return img.addBands(score.rename('cloudScore'))

gee_task_manager.py

@@ -2,12 +2,15 @@
 import ee
 import dill
 import os
+import sys
 import time
 import traceback
 import json
 import copy
-from gevent import monkey
-monkey.patch_all()
+
+# will need when using outside of 
+# from gevent import monkey
+# monkey.patch_all()
 
 from gevent.fileobject import FileObjectThread
 from gevent.queue import Queue, Empty

helpers.py

@@ -17,7 +17,15 @@ def exportImageCollectionToGCS(imgC, bucket=None, resolution=10, start=False):
 
     return(task_ids)
 
-def exportImageToGCS(img=None, roi=None, bucket=None, filename=None, dest_path=None, resolution=10, start=False):
+def exportImageToGCS(img=None, roi=None, bucket=None, filename=None, dest_path=None, resolution=10, start=True, sensor_name=None):
+    ## same as in the JS version
+
+    if sensor_name == 'copernicus/s2':
+        img = img.select(['B4', 'B3', 'B2'])
+    elif sensor_name == 'copernicus/s2_sr':
+        img = img.select(['TCI_R', 'TCI_G', 'TCI_B'])
+
+#     print(img.getInfo())
 
     export = ee.batch.Export.image.toCloudStorage(
       image=img,
@@ -28,19 +36,58 @@ def exportImageToGCS(img=None, roi=None, bucket=None, filename=None, dest_path=N
       bucket=bucket,
       maxPixels=1e13
     )
+
+#     print()
+
+    if start:
+        export.start()
+
+
+
+#     print("exporting final image")
+    return export
+
+
+def exportImageToGDrive(img=None, roi=None, drive_folder=None, filename=None, dest_path=None, resolution=10, start=True):
+
+#downConfig = {
+    #'scale': 10, 
+    #"maxPixels": 1.0E13, 
+    #'driveFolder': 'image3',
+    #"driveFileNamePrefix":str(i)
+#}  # scale means resolution.
+    #image_to_dl = ee.Image(image_dl_list.get(i))
+    # img_2 = image_to_dl.select('B.+')
+    #img_2 = image_to_dl.select(RGB)
+    #name = img_2.getInfo()["id"].split("/")[-1]
+    #print("Image to Download:",name)
+    #task = ee.batch.Export.image(img_2, name, downConfig)
+    #task.start()
+
+    downConfig = {
+        'scale': resolution,
+        'region': roi,
+        'driveFileNamePrefix': dest_path,
+        'driveFolder': drive_folder,
+        'maxPixels': 1e13,
+    }
+
+    export = ee.batch.Export.image(img, filename, downConfig)
 
     if start:
         export.start()
 
     return(export)
 
 def rescale(img, exp, thresholds):
+    ## identical to the javascript
     return img.expression(exp, {"img": img}) \
               .subtract(thresholds[0]) \
               .divide(thresholds[1] - thresholds[0])
 
 def dilatedErossion(score, dilationPixels=3, erodePixels=1.5):
   # Perform opening on the cloud scores
+  ## same as the JS version
   score = score \
             .reproject('EPSG:4326', None, 20) \
             .focal_min(radius=erodePixels, kernelType='circle', iterations=3) \
@@ -53,13 +100,18 @@ def dilatedErossion(score, dilationPixels=3, erodePixels=1.5):
 #        input: imgC - Image collection including "cloudScore" band for each image
 #              threshBest - A threshold percentage to select the best image. This image is used directly as "cloudFree" if one exists.
 #        output: A single cloud free mosaic for the region of interest
+
 def mergeCollection(imgC, keepThresh=5, filterBy='CLOUDY_PERCENTAGE', filterType='less_than', mosaicBy='cloudShadowScore'):
     # Select the best images, which are below the cloud free threshold, sort them in reverse order (worst on top) for mosaicing
+    ## same as the JS version
     best = imgC.filterMetadata(filterBy, filterType, keepThresh).sort(filterBy, False)
+    #print('Info on first image of collection:', imgC.first().getInfo())
     filtered = imgC.qualityMosaic(mosaicBy)
 
     # Add the quality mosaic to fill in any missing areas of the ROI which aren't covered by good images
     newC = ee.ImageCollection.fromImages( [filtered, best.mosaic()] )
+
+    print("collection merged")
 
     return ee.Image(newC.mosaic())
 
@@ -76,6 +128,8 @@ def calcCloudCoverage(img, cloudThresh=0.2):
               )
 
     roi = ee.Geometry(img.get('ROI'))
+    #line below to used debug issue with export tile pipeline
+#     roi = img.geometry()
 
     intersection = roi.intersection(imgPoly, ee.ErrorMargin(0.5))
     cloudMask = img.select(['cloudScore']).gt(cloudThresh).clip(roi).rename('cloudMask')
@@ -87,15 +141,50 @@ def calcCloudCoverage(img, cloudThresh=0.2):
       geometry=roi,
       scale=10,
       maxPixels=1e12,
+      ## bottom two not in the javascript version
       bestEffort=True,
       tileScale=16
     )
 
+    ## maxAreaError not in the javascript version, which uses the default
+    ## for the .area function calls
     maxAreaError = 10
     cloudPercent = ee.Number(stats.get('cloudMask')).divide(imgPoly.area(maxAreaError)).multiply(100)
     coveragePercent = ee.Number(intersection.area(maxAreaError)).divide(roi.area(maxAreaError)).multiply(100)
     cloudPercentROI = ee.Number(stats.get('cloudMask')).divide(roi.area(maxAreaError)).multiply(100)
 
+    img = img.set('CLOUDY_PERCENTAGE', cloudPercent)
+    img = img.set('ROI_COVERAGE_PERCENT', coveragePercent)
+    img = img.set('CLOUDY_PERCENTAGE_ROI', cloudPercentROI)
+
+    print("calculated cloud coverage values")
+
+    return img
+
+def calcCloudCoverage_java(img, cloudThresh=0.2):
+    imgPoly = ee.Algorithms.GeometryConstructors.Polygon(
+              ee.Geometry( img.get('system:footprint') ).coordinates()
+              )
+
+    roi = ee.Geometry(img.get('ROI'))
+
+
+    intersection = roi.intersection(imgPoly, ee.ErrorMargin(0.5))
+    cloudMask = img.select(['cloudScore']).gt(cloudThresh).clip(roi).rename('cloudMask')
+
+    cloudAreaImg = cloudMask.multiply(ee.Image.pixelArea())
+
+    stats = cloudAreaImg.reduceRegion(
+      reducer=ee.Reducer.sum(),
+      geometry=roi,
+      scale=10,
+      maxPixels=1e12,
+    )
+
+    cloudPercent = ee.Number(stats.get('cloudMask')).divide(imgPoly.area()).multiply(100)
+    coveragePercent = ee.Number(intersection.area(maxAreaError)).divide(roi.area()).multiply(100)
+    cloudPercentROI = ee.Number(stats.get('cloudMask')).divide(roi.area()).multiply(100)
+
     img = img.set('CLOUDY_PERCENTAGE', cloudPercent)
     img = img.set('ROI_COVERAGE_PERCENT', coveragePercent)
     img = img.set('CLOUDY_PERCENTAGE_ROI', cloudPercentROI)
@@ -104,3 +193,38 @@ def calcCloudCoverage(img, cloudThresh=0.2):
 
 def clipToROI(x, roi):
     return x.clip(roi).set('ROI', roi)
+
+def inject_B10(img):
+    ee_img_ind = img.get('system:index')
+    coll = ee.ImageCollection('COPERNICUS/S2')\
+        .filterMetadata('system:index', 'equals', ee_img_ind)
+
+    L1C_img = coll.first()
+
+    B10 = L1C_img.select('B10')
+
+    return img.addBands(B10)
+
+## new function found in the javascript
+def computeQualityScore(img):
+    score = img.select(['cloudScore']).max(img.select(['shadowScore']))
+
+    score = score.reproject('EPSG:4326', None, 20).reduceNeighborhood(
+        reducer=ee.Reducer.mean(),
+        kernel=ee.Kernel.square(5)
+    )
+
+    score = score.multiply(-1)
+
+    print("computed quality score")
+
+    return img.addBands(score.rename('cloudShadowScore'))
+
+## new function found in the javascript
+def uniqueValues(collection, field):
+    values = ee.Dictionary(
+        collection.reduceColumns(ee.Reducer.frequencyHistogram(), [field])\
+            .get('histogram')
+    ).keys()
+
+    return values

project_shadows.py

@@ -1,61 +1,86 @@
 import ee
+# original code did not import math
+import math
 
-from helpers import dilatedErossion
+from .helpers import dilatedErossion
 
 # Implementation of Basic cloud shadow shift
 # Author: Gennadii Donchyts
 # License: Apache 2.0
 # Modified by Lloyd Hughes to reduce spurious cloud shadow masks
-def sentinel2ProjectShadows(image, cloudHeights=list(range(200, 10000, 250)), cloudThresh=0.2, irSumThresh=0.3, ndviThresh=-0.1):
-  meanAzimuth = image.get('MEAN_SOLAR_AZIMUTH_ANGLE')
-  meanZenith = image.get('MEAN_SOLAR_ZENITH_ANGLE')
 
-  cloudHeights = ee.List(cloudHeights)
 
-  cloudMask = image.select(['cloudScore']).gt(cloudThresh)
-
-  #Find dark pixels
-  darkPixelsImg = image.select(['B8','B11','B12']) \
+def sentinel2ProjectShadows(
+    image,
+    cloudHeights=list(
+        range(
+            200,
+            10000,
+            250)),
+            cloudThresh=0.2,
+            irSumThresh=0.3,
+            ndviThresh=-
+            0.1):
+      # this function was not used anywhere else in the code.
+      # otherwise, same as the JS version
+
+    meanAzimuth = image.get('MEAN_SOLAR_AZIMUTH_ANGLE')
+    meanZenith = image.get('MEAN_SOLAR_ZENITH_ANGLE')
+
+    cloudHeights = ee.List(cloudHeights)
+    cloudMask = image.select(['cloudScore']).gt(cloudThresh)
+
+      # Find dark pixels
+    darkPixelsImg = image.select(['B8', 'B11', 'B12']) \
                     .divide(10000) \
                     .reduce(ee.Reducer.sum())
 
-  ndvi = image.normalizedDifference(['B8','B4'])
-  waterMask = ndvi.lt(ndviThresh)
+    ndvi = image.normalizedDifference(['B8', 'B4'])
+    waterMask = ndvi.lt(ndviThresh)
 
-  darkPixels = darkPixelsImg.lt(irSumThresh)
+    darkPixels = darkPixelsImg.lt(irSumThresh)
 
-  # Get the mask of pixels which might be shadows excluding water
-  darkPixelMask = darkPixels.And(waterMask.Not())
-  darkPixelMask = darkPixelMask.And(cloudMask.Not())
+    # Get the mask of pixels which might be shadows excluding water
+    darkPixelMask = darkPixels.And(waterMask.Not())
+    darkPixelMask = darkPixelMask.And(cloudMask.Not())
 
-  #Find where cloud shadows should be based on solar geometry
-  #Convert to radians
-  azR = ee.Number(meanAzimuth).add(180).multiply(math.pi).divide(180.0)
-  zenR = ee.Number(meanZenith).multiply(math.pi).divide(180.0)
+    # Find where cloud shadows should be based on solar geometry
+    # Convert to radians
+    azR = ee.Number(meanAzimuth).add(180).multiply(math.pi).divide(180.0)
+    zenR = ee.Number(meanZenith).multiply(math.pi).divide(180.0)
 
-  def findShadows(cloudHeight):
-    cloudHeight = ee.Number(cloudHeight)
+    def findShadows(cloudHeight):
+        cloudHeight = ee.Number(cloudHeight)
 
-    shadowCastedDistance = zenR.tan().multiply(cloudHeight)#Distance shadow is cast
-    x = azR.sin().multiply(shadowCastedDistance).multiply(-1)#.divide(nominalScale)#X distance of shadow
-    y = azR.cos().multiply(shadowCastedDistance).multiply(-1)#Y distance of shadow
-    #return cloudMask.changeProj(cloudMask.projection(), cloudMask.projection().translate(x, y))
-    return image.select(['cloudScore']).displace(ee.Image.constant(x).addBands(ee.Image.constant(y)))
+        shadowCastedDistance = zenR.tan().multiply(
+            cloudHeight)  # Distance shadow is cast
+        # .divide(nominalScale)#X distance of shadow
+        x = azR.sin().multiply(shadowCastedDistance).multiply(-1)
+        y = azR.cos().multiply(shadowCastedDistance).multiply(-1)  # Y distance of shadow
+        # return cloudMask.changeProj(cloudMask.projection(),
+        # cloudMask.projection().translate(x, y))
+        return image.select(
+        ['cloudScore']).displace(
+        ee.Image.constant(x).addBands(
+            ee.Image.constant(y)))
 
-  #Find the shadows
-  shadows = cloudHeights.map( findShadows )
+    # Find the shadows
+    shadows = cloudHeights.map(findShadows)
 
-  shadowMasks = ee.ImageCollection.fromImages(shadows)
-  shadowMask = shadowMasks.mean()
+    shadowMasks = ee.ImageCollection.fromImages(shadows)
+    shadowMask = shadowMasks.mean()
 
-  # #Create shadow mask
-  shadowMask = dilatedErossion(shadowMask.multiply(darkPixelMask))
+    # #Create shadow mask
+    shadowMask = dilatedErossion(shadowMask.multiply(darkPixelMask))
 
-  shadowScore = shadowMask.reduceNeighborhood(
+    shadowScore = shadowMask.reduceNeighborhood(
         reducer=ee.Reducer.max(),
         kernel=ee.Kernel.square(1)
       )
 
-  image = image.addBands(shadowScore.rename(['shadowScore']))
+    image = image.addBands(shadowScore.rename(['shadowScore']))
+
+
+    print("completed shadow score")
 
-  return image
+    return image