I have two images with different projections, sizes, resolutions and area covered. Let's call these two images 'MPF' - which covers a huge part of the Arctic reiong, and 'S1' (sentinel-1) which covers a much smaller area, north of Greenland. I want to stack them making sure they are overlapping correctly. When I do the following way, the result are two images stacked but looks like it simply stacks them without respecting the actual geo-location or coordinates. I tried to ways, both failing:
FIRST TRIAL
import rasterio
import os
file_list = [MPF, s1]
# Read metadata of the first file
with rasterio.open(file_list[0]) as src0:
meta = src0.meta
# Update meta to reflect the number of layers
meta.update(count=len(file_list) + 1) # +1 to include the first band
# Specify the output directory and file name
output_directory = 'outputpath'
output_tiff_stack_rasterio = os.path.join(output_directory, 'rasterioStackTrial1.tif')
# Create the output directory if it doesn't exist
os.makedirs(output_directory, exist_ok=True)
# Create a new multi-band TIFF to stack the files
with rasterio.open(out_tiff, 'w', **meta) as dst:
for id, layer in enumerate(file_list, start=1):
with rasterio.open(layer) as src1:
dst.write_band(id, src1.read(1))
print(f"Stacked TIFF file saved as: {output_tiff_stack_rasterio}")
import rasterio
import matplotlib.pyplot as plt
# Path to the stacked TIFF file
stacked_tiff = output_tiff_stack_rasterio
# Open the stacked TIFF file
with rasterio.open(stacked_tiff) as src:
# Read both bands
band1 = src.read(1)
band2 = src.read(2)
# Create a figure with two subplots for the bands
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
# Display band 1
im1 = ax1.imshow(band1, cmap='viridis')
ax1.set_title('MPF - 2019/19/09')
# Display band 2
im2 = ax2.imshow(band2, cmap='viridis')
ax2.set_title('Sentinel-1 - 2019/19/09')
# Add colorbars for both subplots
cbar1 = fig.colorbar(im1, ax=ax1)
cbar2 = fig.colorbar(im2, ax=ax2)
plt.show()
print('Name of output is:', output_tiff_stack_rasterio)
Resulting image is below. The image on the left covers the whole Arctic, the one on the right is just a tiny portion of it.

SECOND TRIAL:
########## TRYING TO STACK WITH RASTERIO TRIAL 2########################
import rasterio
from rasterio.warp import reproject, Resampling
import numpy as np
# Load the two images
with rasterio.open(image1_path) as src1, rasterio.open(image2_path) as src2:
# Check CRS (coordinate reference system) of each image
print("Image 1 CRS:", src1.crs)
print("Image 2 CRS:", src2.crs)
# Reproject the second image to match the CRS of the first image
src2 = src2.read(
out_shape=src1.shape,
resampling=Resampling.bilinear
)
# Create an output file with two bands
output_path = '/home/b1782dd2-17b6-46da-a896-53dbedd4b60a/NCfiles_LEE/rasterioStackTrial3'
with rasterio.open(output_path, 'w', driver='GTiff', width=src1.width, height=src1.height, count=2, dtype=src1.dtypes[0]) as dst:
# Write the data from each source image to the corresponding band in the output
dst.write(src1.read(1), 1)
dst.write(src2[0], 2) # Assuming src2 has shape (1, height, width)
# Update the metadata for the output image
dst.crs = src1.crs
dst.transform = src1.transform
print("Merged image is saved to", output_path)
rasterio2 = output_path
# Open the stacked TIFF file
with rasterio.open(rasterio2) as src:
# Read both bands
band1 = src.read(1)
band2 = src.read(2)
# Create a figure with two subplots for the bands
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
# Display band 1
im1 = ax1.imshow(band1, cmap='viridis')
ax1.set_title('MPF - 2019/19/09')
# Display band 2
im2 = ax2.imshow(band2, cmap='viridis')
ax2.set_title('Sentinel-1 - 2019/19/09')
# Add colorbars for both subplots
cbar1 = fig.colorbar(im1, ax=ax1)
cbar2 = fig.colorbar(im2, ax=ax2)
The resulting image, for some reason the band refering to the first image goes blank...
