diff --git a/cronModel.py b/cronModel.py new file mode 100644 index 0000000..a325064 --- /dev/null +++ b/cronModel.py @@ -0,0 +1,330 @@ +import ee +import sys +from datetime import date +import math, time + +# initialize EE +ee.Initialize() + +def submitJob(item,mpl): + + # set some ee.Model parameters + bands = ['VH_after0','VH_after1','VH_before0', 'VH_before1','VH_before2','VV_after0','VV_after1','VV_before0', 'VV_before1', 'VV_before2'] + PROJECT = 'cipalawan'; + MODEL_NAME = 'alerts'; + VERSION_NAME = 'sarModelv18' + + # Load the trained model and use it for prediction. + model = ee.Model.fromAiPlatformPredictor( + projectName= PROJECT, + modelName= MODEL_NAME, + version= VERSION_NAME, + inputTileSize= [128,128], + inputOverlapSize= [16,16], + proj= ee.Projection('EPSG:4326').atScale(10), + fixInputProj= True, + outputBands= {'landclass': {'type': ee.PixelType.float(),'dimensions': 1 }}) + MODE = 'DESCENDING' + + # Get the projection that is needed for the study area + projection = ee.Projection('EPSG:32650') + + + # Import Sentinel-1 Collection + s1 = ee.ImageCollection('COPERNICUS/S1_GRD')\ + .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))\ + .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))\ + .filter(ee.Filter.eq('orbitProperties_pass', MODE))\ + .filter(ee.Filter.eq('instrumentMode', 'IW'))\ + .map(erodeGeometry)\ + .map(terrainCorrection)\ + .map(applySpeckleFilter)\ + .map(addRatio) + + # get the image that hasnt been processed + img = ee.Image(s1.filter(ee.Filter.eq("system:index",item)).first()) + # get the geometry + geom = img.geometry() + # get the properties + prop = img.toDictionary() + # get the name of the image + name = img.get("system:index").getInfo() + # get the timestamp + timeStamp = img.get("system:time_start") + date = ee.Date(timeStamp) + + # filter for area of the image + s1Item = s1.filterBounds(geom) + # get the before images + beforeSeries = createSeriesBefore(s1Item,date.advance(-1,"days")) + # get the after image + afterSeries = createSeriesAfter(s1Item,date.advance(1,"days")).select(["VV_after1","VH_after1"]) + # this is the current image + after = img.select(["VV","VH"],["VV_after0","VH_after0"]) + + # combine images into a single image + image = beforeSeries.addBands(after).addBands(afterSeries).unmask(0,False) + # convert to float for include_preprocessing + image = image.select(bands).toFloat() + + # run the cnn + prediction = ee.Image(model.predictImage(image.toArray()).arrayFlatten([["landclass","other"]]).toFloat()) + prediction = prediction.select("landclass").multiply(100).toInt().clip(mpl) + + # set the export location + outputName = "projects/cipalawan/assets/alertsv4/" + name + + geom = geom.transform("EPSG:4326",0.01) + output = ee.Image(prediction).clip(geom) + output = clipEdge(output) + output = output.set(prop).set("system:time_start",timeStamp) + + # create the export task + export_task = ee.batch.Export.image.toAsset(image=output, description="alerts "+ name, assetId=outputName,region=geom.getInfo()['coordinates'], maxPixels=1e13,scale=10 ) + # run the export task + export_task.start() + + + +def createSeriesBefore(collection,date,iters=3,nday =12): + + iterations = list(range(1,iters*nday,nday)) + names = list(["_before{:01d}".format(x) for x in range(0,iters,1)]) + #print(iterations) + #print(names) + + def returnCollection(day,name): + #print(day,name) + start = ee.Date(date).advance(-day,"days").advance(-nday,"days") + end = ee.Date(date).advance(-day,"days") + bandNames = ["VV"+name,"VH"+name] + #print(bandNames) + return ee.Image(collection.filterDate(start,end).mean())\ + .select(["VV","VH"],bandNames)\ + .set("system:time_start",start) + + return toBands(ee.ImageCollection.fromImages(list(map(returnCollection,iterations,names)))) + +def createSeriesAfter(collection,date,iters=2,nday =12): + + + iterations = list(range(1,iters*nday,nday)) + names = list(["_after{:01d}".format(x) for x in range(0,iters,1)]) + + def returnCollection(day,name): + start = ee.Date(date).advance(day,"days") + end = ee.Date(date).advance(day,"days").advance(nday,"days") + bandNames = ["VV"+name,"VH"+name] + + return ee.Image(collection.filterDate(start,end).mean())\ + .select(["VV","VH"],bandNames)\ + .set("system:time_start",start) + + return toBands(ee.ImageCollection.fromImages(list(map(returnCollection,iterations,names)))) +# +# * Clips 5km edges + +def erodeGeometry(image): + return image.clip(image.geometry().buffer(-1000)) + +def clipEdge(image): + return image.clip(image.geometry().buffer(-500)) + + +def applySpeckleFilter(img): + + vv = img.select('VV') + vh = img.select('VH') + vv = speckleFilter(vv).rename('VV'); + vh = speckleFilter(vh).rename('VH'); + return ee.Image(ee.Image.cat(vv,vh).copyProperties(img,['system:time_start'])).clip(img.geometry()).copyProperties(img); + + +def speckleFilter(image): + """ apply the speckle filter """ + ksize = 3 + enl = 7; + + geom = image.geometry() + + # Convert image from dB to natural values + nat_img = toNatural(image); + + # Square kernel, ksize should be odd (typically 3, 5 or 7) + weights = ee.List.repeat(ee.List.repeat(1,ksize),ksize); + + # ~~(ksize/2) does integer division in JavaScript + kernel = ee.Kernel.fixed(ksize,ksize, weights, (ksize//2), (ksize//2), False); + + # Get mean and variance + mean = nat_img.reduceNeighborhood(ee.Reducer.mean(), kernel); + variance = nat_img.reduceNeighborhood(ee.Reducer.variance(), kernel); + + # "Pure speckle" threshold + ci = variance.sqrt().divide(mean);# square root of inverse of enl + + # If ci <= cu, the kernel lies in a "pure speckle" area -> return simple mean + cu = 1.0/math.sqrt(enl); + + # If cu < ci < cmax the kernel lies in the low textured speckle area + # -> return the filtered value + cmax = math.sqrt(2.0) * cu; + + alpha = ee.Image(1.0 + cu*cu).divide(ci.multiply(ci).subtract(cu*cu)); + b = alpha.subtract(enl + 1.0); + d = mean.multiply(mean).multiply(b).multiply(b).add(alpha.multiply(mean).multiply(nat_img).multiply(4.0*enl)); + f = b.multiply(mean).add(d.sqrt()).divide(alpha.multiply(2.0)); + + # If ci > cmax do not filter at all (i.e. we don't do anything, other then masking) + # Compose a 3 band image with the mean filtered "pure speckle", + # the "low textured" filtered and the unfiltered portions + out = ee.Image.cat(toDB(mean.updateMask(ci.lte(cu))),toDB(f.updateMask(ci.gt(cu)).updateMask(ci.lt(cmax))),image.updateMask(ci.gte(cmax))); + + return out.reduce(ee.Reducer.sum()).clip(geom); + + +def toNatural(img): + """Function to convert from dB to natural""" + return ee.Image(10.0).pow(img.select(0).divide(10.0)); + +def toDB(img): + """ Function to convert from natural to dB """ + return ee.Image(img).log10().multiply(10.0); + + +def toBands(collection): + + def createStack(img,prev): + return ee.Image(prev).addBands(img) + + stack = ee.Image(collection.iterate(createStack,ee.Image(1))) + stack = stack.select(ee.List.sequence(1, stack.bandNames().size().subtract(1))); + return stack; + +# Implementation by Andreas Vollrath (ESA), inspired by Johannes Reiche (Wageningen) +def terrainCorrection(image): + date = ee.Date(image.get('system:time_start')) + imgGeom = image.geometry() + srtm = ee.Image('USGS/SRTMGL1_003').clip(imgGeom) # 30m srtm + sigma0Pow = ee.Image.constant(10).pow(image.divide(10.0)) + + #Article ( numbers relate to chapters) + #2.1.1 Radar geometry + theta_i = image.select('angle') + phi_i = ee.Terrain.aspect(theta_i).reduceRegion(ee.Reducer.mean(), theta_i.get('system:footprint'), 1000).get('aspect') + + #2.1.2 Terrain geometry + alpha_s = ee.Terrain.slope(srtm).select('slope') + phi_s = ee.Terrain.aspect(srtm).select('aspect') + + # 2.1.3 Model geometry + # reduce to 3 angle + phi_r = ee.Image.constant(phi_i).subtract(phi_s) + + #convert all to radians + phi_rRad = phi_r.multiply(math.pi / 180) + alpha_sRad = alpha_s.multiply(math.pi / 180) + theta_iRad = theta_i.multiply(math.pi / 180) + ninetyRad = ee.Image.constant(90).multiply(math.pi / 180) + + # slope steepness in range (eq. 2) + alpha_r = (alpha_sRad.tan().multiply(phi_rRad.cos())).atan() + + # slope steepness in azimuth (eq 3) + alpha_az = (alpha_sRad.tan().multiply(phi_rRad.sin())).atan() + + # local incidence angle (eq. 4) + theta_lia = (alpha_az.cos().multiply((theta_iRad.subtract(alpha_r)).cos())).acos() + theta_liaDeg = theta_lia.multiply(180 / math.pi) + + # 2.2 + # Gamma_nought_flat + gamma0 = sigma0Pow.divide(theta_iRad.cos()) + gamma0dB = ee.Image.constant(10).multiply(gamma0.log10()) + ratio_1 = gamma0dB.select('VV').subtract(gamma0dB.select('VH')) + + # Volumetric Model + nominator = (ninetyRad.subtract(theta_iRad).add(alpha_r)).tan() + denominator = (ninetyRad.subtract(theta_iRad)).tan() + volModel = (nominator.divide(denominator)).abs() + + # apply model + gamma0_Volume = gamma0.divide(volModel) + gamma0_VolumeDB = ee.Image.constant(10).multiply(gamma0_Volume.log10()) + + # we add a layover/shadow maskto the original implmentation + # layover, where slope > radar viewing angle + alpha_rDeg = alpha_r.multiply(180 / math.pi) + layover = alpha_rDeg.lt(theta_i); + + # shadow where LIA > 90 + shadow = theta_liaDeg.lt(85) + + # calculate the ratio for RGB vis + ratio = gamma0_VolumeDB.select('VV').subtract(gamma0_VolumeDB.select('VH')) + + output = gamma0_VolumeDB.addBands(ratio).addBands(alpha_r).addBands(phi_s).addBands(theta_iRad)\ + .addBands(layover).addBands(shadow).addBands(gamma0dB).addBands(ratio_1) + + return output.select(['VV', 'VH'], ['VV', 'VH']).set("system:time_start",date).clip(imgGeom ).copyProperties(image) + + +# add the ratio band to the image +def addRatio(img): + geom = img.geometry() + vv = toNatural(img.select(['VV'])).rename(['VV']); + vh = toNatural(img.select(['VH'])).rename(['VH']); + ratio = vh.divide(vv).rename(['ratio']); + return ee.Image(ee.Image.cat(vv,vh,ratio).copyProperties(img,['system:time_start'])).clip(geom).copyProperties(img); + + + +# import the administrative boundaries +WDPAareas = ee.FeatureCollection("WCMC/WDPA/current/polygons"); +mpl = WDPAareas.filter(ee.Filter.eq("ORIG_NAME","Mt. Mantalingahan Protected Landscape")).geometry().buffer(1000); +mpl= ee.FeatureCollection("projects/cipalawan/assets/palawanIsland").geometry().buffer(1000); + +# we only have descening imagery +MODE = 'DESCENDING' + +region = mpl + +# Declare start and end period +cdate = date.today() +start = "2022-10-01" +end = "2022-11-15" + +s1Collection = ee.ImageCollection("COPERNICUS/S1_GRD").filterDate(start,end)\ + .filterBounds(region)\ + .sort("system:time_start")\ + .filter(ee.Filter.eq('orbitProperties_pass', MODE))\ + .sort("system:time_start",True)\ + .aggregate_histogram("system:index").getInfo() + + +gplCollection = ee.ImageCollection("projects/cipalawan/assets/alertsv4").aggregate_histogram("system:index").getInfo() + + +# add s1 filenames in colelction to two different lists +s1List1 = [] +s1List2 = [] +for item in s1Collection: + s1List1.append(item) + s1List2.append(item) + +# add the processed imagery to a list +gplList = [] +for item in gplCollection: + gplList.append(item) + +# remove images that have already been processed from the list +for x in s1List1: + for z in gplList: + if x == z: + s1List2.remove(x) + + +# process imagery that has not been processed +for item in s1List2: + print(counter,item) + submitJob(item,mpl) diff --git a/mergeVectorAlerts.py b/mergeVectorAlerts.py new file mode 100644 index 0000000..826516a --- /dev/null +++ b/mergeVectorAlerts.py @@ -0,0 +1,34 @@ +import ee, os +import subprocess +import datetime + + +ee.Initialize() + +def system_call(command): + p = subprocess.Popen([command], stdout=subprocess.PIPE, shell=True,text=True) + return p.stdout.read() + +batcmd="earthengine ls projects/cipalawan/assets/palawanAlerts/" + +result = system_call(batcmd).split("\n") + +ft = ee.FeatureCollection(result[0]) + +for i in range(1,len(result),1): + if len(result[i]) > 0: + fc = ee.FeatureCollection(result[i]) + ft = ft.merge(fc) + + +# get the date for today +date = datetime.datetime.now() +year = date.year +month = date.month +day = date.day +name = str(year)+str(month)+str(day) + +# create the export task +exportTask = ee.batch.Export.table.toAsset(collection= ee.FeatureCollection(ft),description= name,assetId="projects/cipalawan/assets/palawanAlerts/alerts_"+name) +# run the export task +exportTask.start() diff --git a/setDates.py b/setDates.py new file mode 100644 index 0000000..0531649 --- /dev/null +++ b/setDates.py @@ -0,0 +1,21 @@ +import os, subprocess + + +output = subprocess.check_output("earthengine ls projects/cipalawan/assets/alertsv4",shell=True) + +myList = output.decode('utf8').split("\n") + +counter = 0 +for items in myList: + try: + myName = items.split("/")[4] + year = myName[17:21] + month = myName[21:23] + day = myName[23:25] + #print(items, counter, year, month, day) + + time = str(year) + "-" + str(month) + "-" + str(day) + "T00:00:00" + counter+=1 + os.system("earthengine asset set --time_start " + time + " " + items) + except: + pass diff --git a/vectorizeData.py b/vectorizeData.py new file mode 100644 index 0000000..4bf2625 --- /dev/null +++ b/vectorizeData.py @@ -0,0 +1,112 @@ +import ee, subprocess + +ee.Initialize() + + +# set projection +PRJ = "EPSG:32650" + +# function to calculate minimum mapping unit using connected pixel count. +def apply_mmu (binary, band_name, mmu_value): + + img = binary; + binary = binary.updateMask(binary); + + #// Get the connected pixel count + mmu = binary.connectedPixelCount(mmu_value); + + # Get the connected pixel count + mmu = binary.connectedPixelCount(mmu_value, True).reproject(PRJ, None, 10).gte(mmu_value).toInt8(); + + # Create the cleaned layer + mmu_applied = img.multiply(mmu).unmask(0); + + return mmu_applied.rename([band_name]).toInt8(); + +# system call to get the list of feature collections +def system_call(command): + p = subprocess.Popen([command], stdout=subprocess.PIPE, shell=True,text=True) + return p.stdout.read() + +# import project shapefiles +projectArea = ee.FeatureCollection("projects/cipalawan/assets/Project_Area_Footprint"); +forest = ee.FeatureCollection("projects/cipalawan/assets/LAWIN/FC_PALAWAN2015").filterBounds(projectArea); +ancestraldomains = ee.FeatureCollection("projects/cipalawan/assets/Priority_ancestraldomains"); + +# import the alerts and filter for aoi +gplCollection = ee.ImageCollection("projects/cipalawan/assets/alertsv4").filterBounds(projectArea).filterBounds(forest).aggregate_histogram("system:index").getInfo() + +# command to get list of alerts +batcmd="earthengine ls projects/cipalawan/assets/palawanAlerts/" + +# get the current file list +result = system_call(batcmd).split("\n") + +# add all filenames to a list +ftFiles = [] +for item in result: + try: + ftFiles.append(item.split("/")[4]) + except: + pass + +# add the processed imagery to a list +gplList1 = [] +gplList2 = [] +for item in gplCollection: + gplList1.append(item) + gplList2.append(item) + +# remove images that have already been processed from the list +for x in ftFiles: + for z in gplList2: + if x == z: + gplList1.remove(x) + + +# process imagery that was not processed yet +for item in gplList1: + # get the image + img = ee.Image("projects/cipalawan/assets/alertsv4/"+item) + # get the name + name = img.get("system:index").getInfo() + + # get date information + date = ee.Date(img.get("system:time_start")) + year = date.get("year") + month = date.get("month") + day= date.get("day") + doy = date.format("D") + + # apply a theshold and mmu + img = img.gt(50).selfMask(); + mmu = apply_mmu (img, "mmu", 50); + img = img.updateMask(mmu).clip(forest).clip(projectArea); + + # calculate area for each alert + def calcArea(feat): + area = feat.area(1) + return feat.set("size",area).set("year",year).set("month",month).set("day",day).set("doy",doy) + + # vectorize the data + alertshp = img.reduceToVectors(geometry=projectArea,scale=10,maxPixels=1e13,tileScale=16); + # calculate the area of each polygon + alertshp = alertshp.map(calcArea) + + # create a list to set unique ID to each feature + feature_list = alertshp.toList(alertshp.size()); + indexes = ee.List.sequence(1, alertshp.size()); + + def setId(t): + t = ee.List(t); + f = ee.Feature(t.get(0)); + idx = ee.Number(t.get(1)); + return f.set("ID", idx); + + with_ids = feature_list.zip(indexes).map(setId) + alertshp = ee.FeatureCollection(with_ids); + + # create export task + exportTask = ee.batch.Export.table.toAsset(collection= ee.FeatureCollection(alertshp),description= name,assetId="projects/cipalawan/assets/palawanAlerts/"+name) + # start the task + exportTask.start()