An error occurred while loading the file. Please try again.
-
Guillaume Perréal authored47d6b51c
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)