How to find how many times each pixel is covered in several GeoTIFFs?

Question

I have a few GeoTIFFs which might have overlap, I want to find how many times each pixel is covered in several GeoTIFFs. I thought it is best to create a heat map in Python (using Gdal or Rasterio) showing the number or count of available data for each pixel, in other words, if a pixel is covered by one of the GeoTIFFs put 1, if it is covered by two of them put 2, if three times put 3 and if there is no coverage at all by none put NaN. But how? I worked on merging the GeoTIFFs (called merged.tiff) and then creating a heatmap using the below code, but it is not showing what I want.

with rasterio.open('merged.tiff') as src:
    data = src.read(1)
mask = (data > 0).astype(np.uint8)
pixel_count = np.sum(mask)
heatmap = pixel_count * mask
# Plot the heatmap
plt.imshow(heatmap, cmap='hot', interpolation='nearest')
plt.colorbar(label='Number of Pixels Covered')
plt.title('Pixel Count Heat Map')
plt.show()

Answer 1

It's not clear what you're trying to do. What means covered vs not-covered for example, is the value of 0 what determines that? The fact that you can stack them into a single file by definition means that every pixel is covered. Do you want for example want to count the amount of non-nodata values over all bands/layers?

from osgeo import gdal

ds = gdal.OpenEx("merged.tiff")
data = ds.ReadAsArray()
ds = None

To count the amount of non-zero pixel values over the bands you'll need to specify the axis, otherwise np.sum will return just the total sum over all dimensions, returning a single value.

heatmap = (data > 0).sum(axis=0)

For the special case of 0 using the dedicated Numpy function is probably a little faster but less flexible: heatmap = np.count_nonzero(data, axis=0).

Using the nodata value as defined in the metadata of the file can be done by reading those values before closing the gdal dataset (before ds = None) with:

nodata = np.asarray(list(map(
    lambda i: ds.GetRasterBand(i).GetNoDataValue(), 
    range(1,ds.RasterCount+1),
)))

After which a count can be done with:

assert np.isfinite(nodata).all()
heatmap = (data != nodata.reshape(-1,1,1)).sum(axis=0)

Such a comparison only works if the nodata values aren't np.nan or np.inf, if that's the case a dedicated check can be done with ~np.isnan or np.isfinite.

Finally show the result with something like:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(5,6), facecolor="w", layout="compressed")
im = ax.imshow(heatmap, cmap="hot", interpolation="nearest")
cb = fig.colorbar(
    im, ax=ax, shrink=0.5, orientation="horizontal", 
    label="Number of Pixels Covered",
)
ax.set(title="Pixel Count Heat Map")
plt.show()