diff --git a/scripts/functions_for_rasterio.py b/scripts/functions_for_rasterio.py
index ddc1503f320b5d1910e1dc6717ffe7a8ca29daf8..5ed85cb8f5cf2ada01f1ea7d49484253582e94cb 100644
--- a/scripts/functions_for_rasterio.py
+++ b/scripts/functions_for_rasterio.py
@@ -53,12 +53,13 @@ def pansharpen(m, pan, method='browley', W=0.1, all_data=False):
 
     # get m_bands
     rgbn = np.empty((m.shape[0], m.shape[1], 4))
+    print(1)
     rgbn[:,:,0] = m[:,:, 0] # red
     rgbn[:,:,1] = m[:,:, 1] # green
     rgbn[:,:,2] = m[:,:, 2] # blue
     rgbn[:,:,3] = m[:,:, 3] # NIR-1
 
-    # scaled them
+    # scale them
     rgbn_scaled = np.empty((m.shape[0]*4, m.shape[1]*4, 4))
 
     for i in range(4):
diff --git a/scripts/with_pyotb.py b/scripts/with_pyotb.py
index 80627768340629472a80bd90ff0f9b0f4ebc4f5d..8d99260c83b292310f5b1a6fbcc1c63c8f71dcd8 100644
--- a/scripts/with_pyotb.py
+++ b/scripts/with_pyotb.py
@@ -1,12 +1,27 @@
 import pyotb
+import rasterio
+from scipy import ndimage
 
 xs_path = 'xs.tif'
 pan_path = 'pan.vrt'
 cloud_path = 'cloud_mask.GML'
 
-# Pansharpening
-pxs = pyotb.BundleToPerfectSensor(inp=pan_path, inxs=xs_path, method='rcs') 
-pxs.write('pxs.tif')
+pxs = pyotb.BundleToPerfectSensor(inp='pan.vrt', inxs='xs.tif') 
+
+ndvi = (pxs[:, :, -1] - pxs[:, :, 0]) / (pxs[:, :, -1] + pxs[:, :, 0])
+
+bare_soils = (ndvi < 0.3)
+
+cloud_mask = pyotb.Rasterization('cloud_mask.GML', im=bare_soils)
+
+bare_soils_masked = pyotb.where(cloud_mask == 1, 0, bare_soils)
+
+bare_soils_arr, profile = bare_soils_masked.to_rasterio()  # warning, this loads the whole array in memory!
+
+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)
 
 
-cloud_mask = pyotb.Rasterization(cloud_path, im=pan_path)
diff --git a/scripts/with_rasterio.py b/scripts/with_rasterio.py
index 4d3ea871eb7e6db74aac0ec1d92e144be9db562c..30bf471b1ead83068ec647bde57434a4522211e7 100644
--- a/scripts/with_rasterio.py
+++ b/scripts/with_rasterio.py
@@ -22,5 +22,5 @@ pxs = pansharpen(xs[:, :1000,:1000], pan[:, :4000,:4000])
 
 # Writing pansharpened output
 profile.update(count=4, height=pxs.shape[1], width=pxs.shape[2], driver='GTiff')
-with rasterio.open('pansharpened_rasterio.tif', 'w', **profile) as out_ds:
+with rasterio.open('/tmp/pansharpened_rasterio.tif', 'w', **profile) as out_ds:
     out_ds.write(pxs)