Cloud mask to calculate ndvi

development python
#1

I want to calculate ndvi whith cloud mask using snappy in python. I have tried to use product.getMaskGroup().get(‘opaque_clouds_10m’)
But I don’t know how to apply

mask_group = product.getMaskGroup()
mask_names = mask_group.getNodeNames()

#Select masks
mask_scl =  mask_group.get('opaque_clouds_10m') # ok

#Select Red and NIR bands
red_band = product.getBand('B4')
nir_band = product.getBand('B8')


#output product (ndvi) & new band
output_product = Product('NDVI', 'NDVI', width, height)
ProductUtils.copyGeoCoding(product, output_product)
output_band = output_product.addBand('ndvi', ProductData.TYPE_FLOAT32)

#output writer
output_product_writer = ProductIO.getProductWriter('GeoTiff')
output_product.setProductWriter(output_product_writer)
output_product.writeHeader('ndvi_output2.tif') #ok
output_product.writeHeader(name+'ndvi.tif')


#compute & save ndvi line by line
red_row = numpy.zeros(width, dtype=numpy.float32)
nir_row = numpy.zeros(width, dtype=numpy.float32)

mask_scl2 = jpy.cast(mask_scl, Mask)
slc_data = numpy.zeros(width * height, numpy.uint32)
mask_scl2.readPixels(0, 0, width, height,slc_data)

for y in range(height):
    print("processing line ", y, " of ", height)

    red_row = red_band.readPixels(0, y, width, 1, red_row)
    nir_row = nir_band.readPixels(0, y, width, 1, nir_row)


    red_band.readPixels(0, 0, w, h, slc_data)
    nir_band.readPixels(0, 0, w, h, slc_data)

I need to create a numpy maks condition?

    #create a numpy mask condition
    invalid_mask = numpy.where(mask_array == 0, 1, 0)
    masked_data = numpy.ma.array(data_array, mask=invalid_mask, fill_value=numpy.nan)

    ndvi = numpy.divide((nir_row - red_row),(nir_row + red_row))
    output_band.writePixels(0, y, width, 1, ndvi)


output_product.closeIO()

Can you help me?
Thank you

#2

I think you found already this example, snappy_ndvi_with_masks.py (github.com). But you want to use the scl mask.

I see some problems in your code.

    red_row = red_band.readPixels(0, y, width, 1, red_row)
    nir_row = nir_band.readPixels(0, y, width, 1, nir_row)


    red_band.readPixels(0, 0, w, h, slc_data)
    nir_band.readPixels(0, 0, w, h, slc_data)

Here you read rows AND the full scene. Do wither one or the other.
Also, you read the data of the red and nir band into slc_data.

Here a code snippet which might work. I haven’t tested. But it might help anyway.

mask_scl2 = jpy.cast(mask_scl, Mask)

red_row = numpy.zeros(width, dtype=numpy.float32)
nir_row = numpy.zeros(width, dtype=numpy.float32)
slc_data = numpy.zeros(width, numpy.uint8)

for y in range(height):
    print("processing line ", y, " of ", height)

    red_row = red_band.readPixels(0, y, width, 1, red_row)
    nir_row = nir_band.readPixels(0, y, width, 1, nir_row)
    mask_scl2.readPixels(0, 0, width, 1,slc_data)

    masked_red= numpy.ma.array(r7, mask=slc_data, fill_value=numpy.nan)
    masked_nir = numpy.ma.array(r10, mask=slc_data, fill_value=numpy.nan)

    ndvi = (masked_nir - masked_red) / (masked_nir + masked_red )
    output_band.writePixels(0, y, width, 1, ndvi.filled(numpy.nan))

output_product.closeIO()