ENH: less dummy pansharpening with rasterio

functions_for_rasterio.py

 from scipy import ndimage
+from skimage.transform import rescale
 import numpy as np
 
 
@@ -45,15 +46,20 @@ def stretch(bands, lower_percent=2, higher_percent=98):
 
 def pansharpen(m, pan, method='browley', W=0.1, all_data=False):
 
+    # change convention of pan and ms
+    pan = np.moveaxis(pan, 0, -1)
+    m = np.moveaxis(m, 0, -1)
+
+
     # get m_bands
-    rgbn = np.empty((m.shape[1], m.shape[2], 4))
-    rgbn[:,:,0] = m[0,:,:] # red
-    rgbn[:,:,1] = m[1,:,:] # green
-    rgbn[:,:,2] = m[2,:,:] # blue
-    rgbn[:,:,3] = m[3,:,:] # NIR-1
+    rgbn = np.empty((m.shape[0], m.shape[1], 4))
+    rgbn[:,:,0] = m[:,:, 0] # red
+    rgbn[:,:,1] = m[:,:, 1] # green
+    rgbn[:,:,2] = m[:,:, 2] # blue
+    rgbn[:,:,3] = m[:,:, 3] # NIR-1
 
     # scaled them
-    rgbn_scaled = np.empty((m.shape[1]*4, m.shape[2]*4, 4))
+    rgbn_scaled = np.empty((m.shape[0]*4, m.shape[1]*4, 4))
 
     for i in range(4):
         img = rgbn[:,:,i]
@@ -62,19 +68,19 @@ def pansharpen(m, pan, method='browley', W=0.1, all_data=False):
 
     # check size and crop for pan band
     if pan.shape[0] < rgbn_scaled.shape[0]:
-        rgbn_scaled = rgbn_scaled[:pan.shape[0],:, :]
+        rgbn_scaled = rgbn_scaled[:pan.shape[0], :, :]
     else:
-        pan = pan[:rgbn_scaled.shape[0], :]
+        pan = pan[:rgbn_scaled.shape[0], :, :]
 
     if pan.shape[1] < rgbn_scaled.shape[1]:
-        rgbn_scaled = rgbn_scaled[:,:pan.shape[1], :]
+        rgbn_scaled = rgbn_scaled[:, :pan.shape[1], :]
     else:
-        pan = pan[:,:rgbn_scaled.shape[1]]
+        pan = pan[:, :rgbn_scaled.shape[1], :]
 
-    R = rgbn_scaled[:,:,0]
-    G = rgbn_scaled[:,:,1]
-    B = rgbn_scaled[:,:,2]
-    I = rgbn_scaled[:,:,3]
+    R = rgbn_scaled[:, :, 0]
+    G = rgbn_scaled[:, :, 1]
+    B = rgbn_scaled[:, :, 2]
+    I = rgbn_scaled[:, :, 3]
 
     image = None
 
@@ -105,7 +111,13 @@ def pansharpen(m, pan, method='browley', W=0.1, all_data=False):
         image = np.concatenate([r,g,b,i], axis=2)
 
     if method == 'browley':
-        DNF = (pan - W*I)/(W*R+W*G+W*B)
+        print('pan', pan.shape)
+        print('W', W)
+        print('I', I.shape)
+        print('R', R.shape)
+        print('G', G.shape)
+        print('B', B.shape)
+        DNF = (np.squeeze(pan) - W*I)/(W*R+W*G+W*B)
 
         r = (R * DNF)[:, :, np.newaxis]
         g = (G * DNF)[:, :, np.newaxis]
@@ -119,6 +131,9 @@ def pansharpen(m, pan, method='browley', W=0.1, all_data=False):
         hsv[:,:,2] = pan - I*W
         image = color.hsv2rgb(hsv)
 
+
+    image = np.moveaxis(image, -1, 0)
+
     if all_data:
         return rgbn_scaled, image, I
     else:

with_rasterio.py

@@ -18,7 +18,7 @@ print(xs.shape)
 print(pan.shape)
 
 # Pansharpening
-pxs = dummy_pansharpen(xs[:, :1000,:1000], pan[:, :4000,:4000])
+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')