with_rasterio.py

import rasterio
from functions_for_rasterio import rcs_pansharpen, pansharpen, dummy_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()
    profile = pan_src.profile

print(xs.shape)
print(pan.shape)

# Pansharpening
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')
with rasterio.open('pansharpened_rasterio.tif', 'w', **profile) as out_ds:
    out_ds.write(pxs)