SCRIPTS: some modifs

scripts/functions_for_rasterio.py

@@ -21,8 +21,10 @@ def rcs_pansharpen(xs, pan):
     print(0)
     xs_zoom = ndimage.zoom(xs, zoom=(1, 4, 4), order=0)
     print(1)
+    print(xs_zoom.shape)
     pan_low_pass = ndimage.gaussian_filter(pan, sigma=1)
     print(2)
+    print(pan_low_pass.shape)
 
     pansharpened = xs_zoom * (pan / pan_low_pass)
     print(3)
@@ -86,14 +88,15 @@ def pansharpen(m, pan, method='browley', W=0.1, all_data=False):
     image = None
 
     if method == 'simple_browley':
-        all_in = R+G+B
-        prod = np.multiply(all_in, pan)
+        all_in = R+G+B+I
+        prod = np.multiply(all_in, np.squeeze(pan))
 
-        r = np.multiply(R, pan/all_in)[:, :, np.newaxis]
-        g = np.multiply(G, pan/all_in)[:, :, np.newaxis]
-        b = np.multiply(B, pan/all_in)[:, :, np.newaxis]
+        r = np.multiply(R, np.squeeze(pan)/all_in)[:, :, np.newaxis]
+        g = np.multiply(G, np.squeeze(pan)/all_in)[:, :, np.newaxis]
+        b = np.multiply(B, np.squeeze(pan)/all_in)[:, :, np.newaxis]
+        i = np.multiply(I, np.squeeze(pan)/all_in)[:, :, np.newaxis]
 
-        image = np.concatenate([r,g,b], axis=2)
+        image = np.concatenate([r,g,b,i], axis=2)
 
     if method == 'sample_mean':
         r = 0.5 * (R + pan)[:, :, np.newaxis]
@@ -103,7 +106,7 @@ def pansharpen(m, pan, method='browley', W=0.1, all_data=False):
         image = np.concatenate([r,g,b], axis=2)
 
     if method == 'esri':
-        ADJ = pan-rgbn_scaled.mean(axis=2)
+        ADJ = np.squeeze(pan)-rgbn_scaled.mean(axis=2)
         r = (R + ADJ)[:, :, np.newaxis]
         g = (G + ADJ)[:, :, np.newaxis]
         b = (B + ADJ)[:, :, np.newaxis]

scripts/pyotb_usecase.py

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

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'
-
-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)
+import time
 
+start = time.time()
+pxs = pyotb.BundleToPerfectSensor(inp='pan.vrt', inxs='xs.tif', ram=190000) 
+pxs.write('pxs_pyotb.tif', pixel_type='uint16')
 
+print('ELAPSED TIME:', time.time() - start)

scripts/with_rasterio.py

 import rasterio
 from functions_for_rasterio import rcs_pansharpen, pansharpen, dummy_pansharpen
+from scipy import ndimage
+
+
+import time
+
+start = time.time()
 
 xs_path = 'xs.tif'
 pan_path = 'pan.vrt'
@@ -17,10 +23,24 @@ with rasterio.open(pan_path) as pan_src:
 print(xs.shape)
 print(pan.shape)
 
+""" Just to try the gaussian filter
+pan_low_pass = ndimage.gaussian_filter(pan, sigma=1)
+
+profile.update(driver='GTiff')
+with rasterio.open('pan_gaussian.tif', 'w', **profile) as out_ds:
+    out_ds.write(pan_low_pass)
+"""
+
+print('Read time', time.time() - start)
+
 # Pansharpening
-pxs = pansharpen(xs[:, :1000,:1000], pan[:, :4000,:4000])
+pxs = rcs_pansharpen(xs, pan)
+
+print('After pan sharpen', time.time() - start)
 
 # Writing pansharpened output
 profile.update(count=4, height=pxs.shape[1], width=pxs.shape[2], driver='GTiff')
-with rasterio.open('/tmp/pansharpened_rasterio.tif', 'w', **profile) as out_ds:
+with rasterio.open('pansharpened_rasterio_rcs.tif', 'w', **profile) as out_ds:
     out_ds.write(pxs)
+
+print('ELAPSED TIME', time.time() - start)