Stacking two TIFF images when they have different projections, area covered, resolution

18 Views Asked by At

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. enter image description here

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... enter image description here

0

There are 0 best solutions below