diff --git a/functions_for_rasterio.py b/functions_for_rasterio.py index 304bdc76c1e203d3e719ed47c09641e635e21cdb..ddc1503f320b5d1910e1dc6717ffe7a8ca29daf8 100644 --- a/functions_for_rasterio.py +++ b/functions_for_rasterio.py @@ -1,4 +1,5 @@ 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: diff --git a/with_rasterio.py b/with_rasterio.py index 0c979e4e91cc80496de010c3108bdbed875f21bd..4d3ea871eb7e6db74aac0ec1d92e144be9db562c 100644 --- a/with_rasterio.py +++ b/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')