import pyotb
import rasterio
from scipy import ndimage

# Pansharpening
pxs = pyotb.BundleToPerfectSensor(inp='pan.vrt', inxs='xs.tif') 

# Computing NDVI
ndvi = (pxs[:, :, -1] - pxs[:, :, 0]) / (pxs[:, :, -1] + pxs[:, :, 0])

# Computing a boolean raster, describing the bare soils i.e. pixels without vegetation
bare_soils = (ndvi < 0.3)

# Creating a boolean cloud mask from the GML vector
cloud_mask = pyotb.Rasterization('cloud_mask.GML', im=bare_soils)

# Masking clouds (i.e. assigning to 0) on the result
bare_soils_masked = pyotb.where(cloud_mask == 1, 0, bare_soils)

# Exporting to rasterio
bare_soils_arr, profile = bare_soils_masked.to_rasterio()  # warning, this loads the whole array in memory!

# Gathering the contiguous groups of pixels with an integer index
labeled_bare_soils, _ = ndimage.label(bare_soils_arr)

# Writing the result to disk
with rasterio.open('labeled_bare_soil.tif', 'w', **profile) as f:
    f.write(labeled_bare_soils)