I am using MODIS11A2 LST 8-day average raster .tif data to calculate daily average for the year. For year 2001, I have 44 files (2 files are missing). I want that a pixel should be excluded from the average only when it has NA (nodata value) in all the rasters. Suppose if a pixel at x,y location has NA value only in one raster and has valid data in other rasters, then the average of that pixel should be calculated as sum of value of that pixel in 43 rasters divided by 43 (not 44). I have created following code to do this but my code is excluding the pixel in question from the average even if other rasters have valid data for that pixel.

# Load necessary libraries
library(raster)
library(rgdal)

# Set the path to your folder containing the .tif files
folder_path <- "D:/A/LST/mosaic_clipped_reprojected/test"

# List all .tif files in the folder
tif_files <- list.files(folder_path, pattern = "*.tif$", full.names = TRUE)

# Initialize variables to store the sum of rasters and valid count
total_sum <- NULL
valid_count <- NULL

# Set NoData value for MODIS LST data
nodata_value <- -9999

# Loop through all the tif files and process them
for (tif_file in tif_files) {
  # Read the raster file
  raster_data <- raster(tif_file)

  # Set the NoData value for the raster (this step ensures NoData values are handled 
correctly)
  raster_data[raster_data == nodata_value] <- NA  # Convert NoData to NA

  # Initialize total_sum and valid_count on the first iteration
  if (is.null(total_sum)) {
    total_sum <- raster_data
    valid_count <- !is.na(raster_data)  # Create a mask for valid data (TRUE = valid, 
FALSE = NA)
  } else {
    # Add current raster data to the total_sum, and count valid pixels
    total_sum <- total_sum + raster_data
    valid_count <- valid_count + !is.na(raster_data)  # Count valid pixels (add 1 for 
valid, 0 for NA)
  }
}

# After all rasters are processed, replace invalid count (0 values) with NA to avoid 
dividing by 0
valid_count[valid_count == 0] <- NA

# Calculate the daily average by dividing the total sum by the valid count for each 
pixel
# This ensures that only valid pixels contribute to the sum and the division.
daily_average <- total_sum / valid_count

# Define the reference raster for resampling (choose one of your original rasters)
reference_raster <- raster(tif_files[1])  # Example: using the first raster file as 
reference

# Resample the daily average to the reference raster using bilinear interpolation
resampled_daily_average <- resample(daily_average, reference_raster, method = 
"bilinear")

# Define the output file path for the daily average raster
output_file <- "D:/A/LST/AnnualAvgLST/daily_average_resampled_test1.tif"

# Write the resampled daily average raster to a new file
writeRaster(resampled_daily_average, output_file, format = "GTiff", overwrite = TRUE)

# Print message to indicate the result
cat("Resampled daily average raster has been saved to:", output_file, "\n")

I am using MODIS11A2 LST 8-day average raster .tif data to calculate daily average for the year. For year 2001, I have 44 files (2 files are missing). I want that a pixel should be excluded from the average only when it has NA (nodata value) in all the rasters. Suppose if a pixel at x,y location has NA value only in one raster and has valid data in other rasters, then the average of that pixel should be calculated as sum of value of that pixel in 43 rasters divided by 43 (not 44). I have created following code to do this but my code is excluding the pixel in question from the average even if other rasters have valid data for that pixel.

# Load necessary libraries
library(raster)
library(rgdal)

# Set the path to your folder containing the .tif files
folder_path <- "D:/A/LST/mosaic_clipped_reprojected/test"

# List all .tif files in the folder
tif_files <- list.files(folder_path, pattern = "*.tif$", full.names = TRUE)

# Initialize variables to store the sum of rasters and valid count
total_sum <- NULL
valid_count <- NULL

# Set NoData value for MODIS LST data
nodata_value <- -9999

# Loop through all the tif files and process them
for (tif_file in tif_files) {
  # Read the raster file
  raster_data <- raster(tif_file)

  # Set the NoData value for the raster (this step ensures NoData values are handled 
correctly)
  raster_data[raster_data == nodata_value] <- NA  # Convert NoData to NA

  # Initialize total_sum and valid_count on the first iteration
  if (is.null(total_sum)) {
    total_sum <- raster_data
    valid_count <- !is.na(raster_data)  # Create a mask for valid data (TRUE = valid, 
FALSE = NA)
  } else {
    # Add current raster data to the total_sum, and count valid pixels
    total_sum <- total_sum + raster_data
    valid_count <- valid_count + !is.na(raster_data)  # Count valid pixels (add 1 for 
valid, 0 for NA)
  }
}

# After all rasters are processed, replace invalid count (0 values) with NA to avoid 
dividing by 0
valid_count[valid_count == 0] <- NA

# Calculate the daily average by dividing the total sum by the valid count for each 
pixel
# This ensures that only valid pixels contribute to the sum and the division.
daily_average <- total_sum / valid_count

# Define the reference raster for resampling (choose one of your original rasters)
reference_raster <- raster(tif_files[1])  # Example: using the first raster file as 
reference

# Resample the daily average to the reference raster using bilinear interpolation
resampled_daily_average <- resample(daily_average, reference_raster, method = 
"bilinear")

# Define the output file path for the daily average raster
output_file <- "D:/A/LST/AnnualAvgLST/daily_average_resampled_test1.tif"

# Write the resampled daily average raster to a new file
writeRaster(resampled_daily_average, output_file, format = "GTiff", overwrite = TRUE)

# Print message to indicate the result
cat("Resampled daily average raster has been saved to:", output_file, "\n")

• r
• tiff
• r-raster
• no-data
• r-modis

2 Answers 2

Patrick answered your specific question, but a bigger issue is that your script could be much simpler. There is no need for a loop, and to compute the average, you should use the mean function. I am showing this with "terra" because "raster" is obsolete, but the same would apply to a "raster" script.

library(terra)

folder_path <- "D:/A/LST/mosaic_clipped_reprojected/test"
tif_files <- list.files(folder_path, pattern = "\\.tif$", full.names = TRUE)

rd <- rast(tif_files)
# NAflag(rd) <- -99999  # is this necessary?

daily_average <- mean(rd, na.rm=TRUE)

output_file <- "D:/A/LST/AnnualAvgLST/daily_average_resampled_test1.tif"
writeRaster(daily_average, output_file, overwrite = TRUE)

Setting the NAflag is probably not necessary. There is no point in using resample the way you do as the input and output raster geometries are identical.

The problem in your code is here:

total_sum <- total_sum + raster_data

If either term in the summation is NA then the sum will also be NA. You can solve this easily:

total_sum <- sum(total_sum, raster_data, na.rm = TRUE)