From 9d47e25123a43f96ec010038572948e34d9d1980 Mon Sep 17 00:00:00 2001 From: davidanagy Date: Mon, 16 Nov 2020 15:49:32 -0800 Subject: [PATCH] Project Canopy pipeline v1 --- cloud_score.py | 11 +++- gee_task_manager.py | 7 ++- helpers.py | 126 +++++++++++++++++++++++++++++++++++++++++++- project_shadows.py | 95 +++++++++++++++++++++------------ 4 files changed, 199 insertions(+), 40 deletions(-) diff --git a/cloud_score.py b/cloud_score.py index 00c7912..7f6ede2 100644 --- a/cloud_score.py +++ b/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')) diff --git a/gee_task_manager.py b/gee_task_manager.py index 721fd8e..58a0623 100644 --- a/gee_task_manager.py +++ b/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 diff --git a/helpers.py b/helpers.py index 1a102f6..fe38613 100644 --- a/helpers.py +++ b/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,6 +36,43 @@ 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() @@ -35,12 +80,14 @@ def exportImageToGCS(img=None, roi=None, bucket=None, filename=None, dest_path=N 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 diff --git a/project_shadows.py b/project_shadows.py index 1ae092a..6e607bf 100644 --- a/project_shadows.py +++ b/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