ENH: dummy pansharpening

functions_for_rasterio.py

 from scipy import ndimage
+import numpy as np
 
 
+def dummy_pansharpen(xs, pan):
+    xs_zoom = ndimage.zoom(xs, zoom=(1, 4, 4), order=0)
+
+    total = sum(xs_zoom[i, :, :] for i in range(4))
+
+    pxs = (xs_zoom / total) * pan
+    return pxs
+
+
+#### DO NO WORK ####
 def rcs_pansharpen(xs, pan):
     """
     Implemented by me
     cf https://www.sfpt.fr/hyperspectral/wp-content/uploads/2018/01/May_panorama_pansharpening.pdf
     """
-
-    xs_zoom = ndimage.zoom(xs, zoom=4)
+    print(0)
+    xs_zoom = ndimage.zoom(xs, zoom=(1, 4, 4), order=0)
+    print(1)
     pan_low_pass = ndimage.gaussian_filter(pan, sigma=1)
+    print(2)
 
     pansharpened = xs_zoom * (pan / pan_low_pass)
+    print(3)
     return pansharpened
 
 
 ## NOT USED ANYMORE (can not compare speed with OTB)
-"""
 # functions
 def stretch(bands, lower_percent=2, higher_percent=98):
     out = np.zeros_like(bands)
@@ -110,4 +123,4 @@ def pansharpen(m, pan, method='browley', W=0.1, all_data=False):
         return rgbn_scaled, image, I
     else:
         return image
-"""
\ No newline at end of file
+

with_rasterio.py

 import rasterio
-from functions_for_rasterio import rcs_pansharpen
+from functions_for_rasterio import rcs_pansharpen, pansharpen, dummy_pansharpen
 
 xs_path = 'xs.tif'
 pan_path = 'pan.vrt'
@@ -12,10 +12,15 @@ with rasterio.open(xs_path) as xs_src:
 # Reading panchromatic
 with rasterio.open(pan_path) as pan_src:
     pan = pan_src.read()
+    profile = pan_src.profile
 
 print(xs.shape)
 print(pan.shape)
 
 # Pansharpening
-pxs = rcs_pansharpen(xs, pan)
+pxs = dummy_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:
+    out_ds.write(pxs)