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'
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()
    profile = pan_src.profile

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 = 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('pansharpened_rasterio_rcs.tif', 'w', **profile) as out_ds:
    out_ds.write(pxs)

print('ELAPSED TIME', time.time() - start)