ENH: more scripts

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):

scripts/with_pyotb.py

 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)

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)