Batch processing geoTIFFs: for-loop only exports last dataset

Question

I am working on a batch processing script that goes through a folder of satellite images (geoTIFF format) and recognizes when two or more images were captured on the same date. Then it will interpolate the pixel values onto a standardized 1000m resolution grid (the images all have slightly different pixel sizes and this "regridding" helps to fix any pixel offset issues) and merge the images together. The merged images do not overlap. I am basically just trying to stitch them together at the edges. In the end the script should export all images as geoTIFFs into a different folder.

Here are the main parts of my code. First I define x and y components of my grid for later use:

# creating standardized grid with 1000m pixel resolution:
# calculate 1000m size pixels in degrees lat and lon:
a = 6378137  # WGS84 large half axis (latitude) in meters
b = 6356752  # WGS84 small half axis (longitude) in meters

earth_circum_lat = 2 * math.pi * a  # latitude circumference of the earth
earth_circum_lon = 2 * math.pi * b  # longitude circumference of the earth

deg_lat = 360 / earth_circum_lat * 1000
deg_lon = 360 / earth_circum_lon * 1000

# create x and y components for the grid:
y = np.arange(5., -5., -deg_lat)
x = np.arange(-95., -85., deg_lon)

Then I open the geoTIFFs and and find out, how many different dates have data:

# open all MODIS geoTIFF-files:
folder = "E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/batch"
list_of_paths = glob.glob(folder + '/*.tif', recursive=True)

# make list with all different filenames (dates) in this folder:
modis = []  # initialize empty list for all file names
for i in range(0, np.size(list_of_paths)):
    # files naming convention "cloud_effective_radius_YYYYMMDD_HHMMSS.tif":
    modis.append(list_of_paths[i].split('ius_')[1][0:8])

# find out how many dates are in the folder:
modis = np.unique(modis)  # remove duplicates from array
print(modis)
print('
data from {} different dates in this folder
'.format(np.size(modis)))

Lastly, I interpolate the images to the standardized grid, merge them together and then export them as geoTIFFs into my destination folder:

# batch process for interpolating and merging geoTIFFs:
for i in range(0, np.size(modis)):
    # all files names of this day in one list
    list_files_date = glob.glob(os.path.join(folder, 'cloud_effective_radius_{0}*.tif'.format(modis[i])))

    # if more than one file for one date -> merge them together
    if len(list_files_date) > 1:
        ds = rioxarray.open_rasterio(list_files_date[0], engine='rasterio')
        ds1 = rioxarray.open_rasterio(list_files_date[1], engine='rasterio')

        # regrid data -> interpolate x and y so that all data are on the same grid
        ds_interp = ds.interp(y=y, x=x, method="nearest")
        ds1_interp = ds1.interp(y=y, x=x, method="nearest")

        # merge every data variable one by one
        ds_merged = ds_interp.combine_first(ds1_interp)
        # it combines both datasets but in case both have a value != nan it uses the first value
        # if only one is nan it uses the value != nan

        if np.size(list_files_date) > 2:
            for j in range(2, np.size(list_files_date)):
                ds = rioxarray.open_rasterio(list_files_date[j], engine='rasterio')
                ds_interp = ds.interp(y=y, x=x, method="nearest")

                ds_merged = ds_merged.combine_first(ds_interp)

    # if just one file for one date -> just regrid it
    else:
        list_files_date = (glob.glob(os.path.join(folder, 'cloud_effective_radius_{0}*.tif'.format(modis[i]))))
        ds = rioxarray.open_rasterio(list_files_date[0], engine='rasterio')
        ds_merged = ds.interp(y=y, x=x, method="nearest")

    # export raster as geoTIFFs:
    img_number = 1
    ds_merged.rio.to_raster("E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/2021_06_merged/ds_merged" +
                            str(img_number) + ".tif",
                            driver="GTiff")
    img_number += 1

When I run the script, only one merged geoTIFF is created in the destination folder. This is always the one from the last day, so I assume that the rest gets overwritten, but I am unsure of why this is happening? Is there a problem with my export code? I thought that the img_number += 1 part would fix this issue. Or is it that the ds_merged variable gets overwritten in the for loop?

Any help is greatly appreciated!

Batch processing geoTIFFs: for-loop only exports last dataset

0 Answers