1

With the code below, I computed the NDVI, NDWI and NDWIMOD indices of a R.O.I. in order to detect water areas in that R.O.I.

Now, my goal concerns to compare the water areas with a given shapefile in order to find which areas are correctly detected by computing the indices. I want to do it for each index separately.

Do you have any advice about how I can carry out it? Moreover, How can I save the new areas detected?

The code is the following:

import ee
import ee.mapclient

ee.Initialize()
ee.mapclient.centerMap(15.88, 41.78, 10)

def maskS2clouds(collection): 
    qa = collection.select('QA60')
    cloudBitMask = 1 << 10
    cirrusBitMask = 1 << 11
    mask = qa.bitwiseAnd(cloudBitMask).eq(0)and(qa.bitwiseAnd(cirrusBitMask).eq(0))
    return collection.updateMask(mask).divide(10000)

polygon = ee.Geometry.Polygon([[
[16.440171600655958,42.033470726105556], 
[15.511826874093458,42.12316884368949],
[14.621934295968458,41.689821328579825],
[15.441789032296583,41.00320702735848],
[16.246537567452833,41.36906568596335],
[16.440171600655958,42.033470726105556]]])

collection = ee.ImageCollection('COPERNICUS/S2').filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)).filterDate('2018-01-1','2018-06-30').map(maskS2clouds).filterBounds(polygon)

medianpixel = collection.median()
medianpixelclipped = medianpixel.clip(polygon)

ndvi = collection.mean().normalizedDifference(['B8','B4'])
ndviViz = {'min':0.6, 'max': 1, 'palette': ['0000FF', 'FF0000']}
ndvimasked = ndvi.updateMask(ndvi.lt(0.0))
ee.mapclient.addToMap(ndvimasked, ndviViz, 'Sentinel_2 mosaic')

ndwi2 =  collection.mean().normalizedDifference(['B3','B8'])
ndwi2Viz = {'min': -0.4, 'max': 0.0, 'palette': ['FF0000', '00FF00']}
ndwi2masked = ndwi2.updateMask(ndwi2.gt(-0.125))
ee.mapclient.addToMap(ndwi2masked, ndwi2Viz, 'Sentinel_2 mosaic')

ndwi2mod = collection.mean().normalizedDifference(['B8','B3'])
ndwi2modViz = {'min': 0.0, 'max': 0.4, 'palette': ['FFFFFF', 'FF0000']}
ndwi2modmasked = ndwi2mod.updateMask(ndwi2mod.lt(0.0))
ee.mapclient.addToMap(ndwi2modmasked, ndwi2modViz, 'Sentinel_2 mosaic')

capitanata = ee.FeatureCollection('users/vpaparella96/vasche_puglia')
ee.mapclient.addToMap(capitanata)

I shared to everyone the shapefile.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.