Merge rasters together with overlapping extents without losing the pixels in the overlap area

Question

I am trying to merge 4 rasters together using the rasterio.merge function. All four rasters have the same CRS and resolution, but they all have different, but overlapping extents. The extent of all four rasters overlap each other, but there are no pixels with positive values that overlap from the different rasters. The merge was partially successful as a new raster was formed; however, the merged raster only shows nodata values for pixels that exist in the area where the extents overlap.

Merged raster:

Four input rasters displayed together (north, east, west, and south EU)

I have tried to use the rasterio.merge (as merge) function:

mosaic, output = merge(raster_to_mosaic, method = 'first', nodata = -9999)

and I tried all merge methods (first, last, min, max, count, sum)

I expected to get a raster back which overrided pixels nodata that were overlapping with pixels with real values.

Whole script:

import rasterio as rio
from rasterio.merge import merge
from pathlib import Path


east_rent_filepath = "/merged/east_SSP1_2020_rent_hectare_filtered.tif"
north_rent_filepath = "/merged/north_SSP1_2020_rent_hectare_filtered.tif"
west_rent_filepath = "/merged/west_SSP1_2020_rent_hectare_filtered.tif"
south_rent_filepath = "/merged/south_SSP1_2020_rent_hectare_filtered.tif"
out_path = "/merged/merged_count.tif"


# Open the raster file for reading and writing
with rasterio.open(east_rent_filepath, "r+") as src:
    # Read the data into a numpy array
    data = src.read(1)
    
    # Set the NoData value
    data[data == -9999] = src.nodata
    
    # Update the raster with the new data
    src.write(data, 1)

# Open the raster file for reading and writing
with rasterio.open(north_rent_filepath, "r+") as src:
    # Read the data into a numpy array
    data = src.read(1)
    
    # Set the NoData value
    data[data == -9999] = src.nodata
    
    # Update the raster with the new data
    src.write(data, 1)

# Open the raster file for reading and writing
with rasterio.open(west_rent_filepath, "r+") as src:
    # Read the data into a numpy array
    data = src.read(1)
    
    # Set the NoData value
    data[data == -9999] = src.nodata
    
    # Update the raster with the new data
    src.write(data, 1)

# Open the raster file for reading and writing
with rasterio.open(south_rent_filepath, "r+") as src:
    # Read the data into a numpy array
    data = src.read(1)
    
    # Set the NoData value
    data[data == -9999] = src.nodata
    
    # Update the raster with the new data
    src.write(data, 1)    
print("updated nodata to -9999")    


path = Path('/merged/')
Path('output').mkdir(parents=True, exist_ok=True)
raster_files = list(path.iterdir())
raster_to_mosaic = []

for p in raster_files:
    raster = rio.open(p)
    raster_to_mosaic.append(raster)

mosaic, output = merge(raster_to_mosaic, method = 'first', nodata = -9999)

output_meta = raster.meta.copy()
output_meta.update(
    {"driver": "GTiff",
        "height": mosaic.shape[1],
        "width": mosaic.shape[2],
        "transform": output,
    }
)
with rio.open(out_path, "w", **output_meta) as m:
    m.write(mosaic)

Answer 1

Resolved this issue by a work around. I first clipped all rasters to a common extent in QGIS using the "Clip raster by extent" function (used a raster layer extent which included all the rasters) and then I used the following code to stack the rasters on top of each other and ignoring any NaN values, and found the maximum value across the four rasters:

import rasterio
import numpy as np

# Open invoerrasters
east_extent = r"/east.tif"
west_extent = r"/west.tif"
north_extent = r"/north.tif"
south_extent = r"/south.tif"
combined_path = "C:combined.tif"


with rasterio.open(east_extent) as src_East, rasterio.open(west_extent) as src_West, rasterio.open(north_extent) as src_North, rasterio.open(south_extent) as src_South:
    # Lees metadata van één invoerraster (ervan uitgaande dat beide rasters dezelfde metadata hebben)
    meta = src_East.meta.copy()

    # Werk metadata bij voor het uitvoerraster
    meta.update(driver='GTiff', dtype='float32')

    # Maak uitvoerrasterbestand aan
    out_path = combined_path
    with rasterio.open(out_path, 'w', **meta) as dst:
        # Voer de vermenigvuldiging uit op de invoerrasters
        east = src_East.read(1)  # Lees eerste band van raster1
        west = src_West.read(1)
        north = src_North.read(1)
        south = src_South.read(1)
        
        arrays = [east, south, north, west]

        combined_array = np.stack(arrays)

        res = np.nanmax(combined_array, axis = 0)

        dst.write(res, 1)  # Schrijf het resultaat naar de eerste band

print(f"Uitvoerraster is succesvol geschreven. Hoort bij bestandspad {uitvoer_pad}")