Commit b2fe0f84 authored by Narcon Nicolas's avatar Narcon Nicolas
Browse files

First scripts for comparison of pyotb vs rasterio for pansharpening

parents
No related merge requests found
Showing with 133 additions and 0 deletions
+133 -0
import numpy as np
from skimage.transform import rescale
import skimage.color as color
# functions
def stretch(bands, lower_percent=2, higher_percent=98):
out = np.zeros_like(bands)
for i in range(3):
a = 0
b = 255
c = np.percentile(bands[:,:,i], lower_percent)
d = np.percentile(bands[:,:,i], higher_percent)
t = a + (bands[:,:,i] - c) * (b - a) / (d - c)
t[t<a] = a
t[t>b] = b
out[:,:,i] =t
return out.astype(np.uint8)
def pansharpen(m, pan, method='browley', W=0.1, all_data=False):
# 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
# scaled them
rgbn_scaled = np.empty((m.shape[1]*4, m.shape[2]*4, 4))
for i in range(4):
img = rgbn[:,:,i]
scaled = rescale(img, (4,4))
rgbn_scaled[:,:,i] = scaled
# check size and crop for pan band
if pan.shape[0] < rgbn_scaled.shape[0]:
rgbn_scaled = rgbn_scaled[:pan.shape[0],:, :]
else:
pan = pan[:rgbn_scaled.shape[0], :]
if pan.shape[1] < rgbn_scaled.shape[1]:
rgbn_scaled = rgbn_scaled[:,:pan.shape[1], :]
else:
pan = pan[:,:rgbn_scaled.shape[1]]
R = rgbn_scaled[:,:,0]
G = rgbn_scaled[:,:,1]
B = rgbn_scaled[:,:,2]
I = rgbn_scaled[:,:,3]
image = None
if method == 'simple_browley':
all_in = R+G+B
prod = np.multiply(all_in, 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]
image = np.concatenate([r,g,b], axis=2)
if method == 'sample_mean':
r = 0.5 * (R + pan)[:, :, np.newaxis]
g = 0.5 * (G + pan)[:, :, np.newaxis]
b = 0.5 * (B + pan)[:, :, np.newaxis]
image = np.concatenate([r,g,b], axis=2)
if method == 'esri':
ADJ = pan-rgbn_scaled.mean(axis=2)
r = (R + ADJ)[:, :, np.newaxis]
g = (G + ADJ)[:, :, np.newaxis]
b = (B + ADJ)[:, :, np.newaxis]
i = (I + ADJ)[:, :, np.newaxis]
image = np.concatenate([r,g,b,i], axis=2)
if method == 'browley':
DNF = (pan - W*I)/(W*R+W*G+W*B)
r = (R * DNF)[:, :, np.newaxis]
g = (G * DNF)[:, :, np.newaxis]
b = (B * DNF)[:, :, np.newaxis]
i = (I * DNF)[:, :, np.newaxis]
image = np.concatenate([r,g,b,i], axis=2)
if method == 'hsv':
hsv = color.rgb2hsv(rgbn_scaled[:,:,:3])
hsv[:,:,2] = pan - I*W
image = color.hsv2rgb(hsv)
if all_data:
return rgbn_scaled, image, I
else:
return image
\ No newline at end of file
with_pyotb.py 0 → 100644
import pyotb
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')
cloud_mask = pyotb.Rasterization(cloud_path, im=pan_path)
import rasterio
from functions_for_rasterio import pansharpen
xs_path = 'xs.tif'
pan_path = 'pan.vrt'
cloud_path = 'cloud_mask.GML'
# Reading multispectral
with rasterio.open(xs_path) as xs_src:
xs = xs_src.read()
# Reading panchromatic
with rasterio.open(pan_path) as pan_src:
pan = pan_src.read()
print(xs.shape)
print(pan.shape)
# Pansharpening
pxs = pansharpen(xs, pan)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment