如何排序TRMM每日NetCDF文件、创建栅格堆栈并验证时序与点位准确性

Hey there! Let's work through how to verify both the temporal order of your TRMM raster stack and the accuracy of your point extractions. Here's a practical, step-by-step breakdown:

1. Verify the Temporal Order of Your Raster Stack #

TRMM NetCDF files store timestamp information in the time variable, which we can extract to cross-check against your stack's layer order.

Step 1: Extract Dates Directly from NetCDF Files #

Use the ncdf4 package to pull the actual timestamp from each file, then map it to your stack layers:

library(ncdf4)
library(raster)

# Extract dates from each TRMM nc file
trmm_dates <- lapply(trmm_files, function(file_path) {
  nc_conn <- nc_open(file_path)
  # Pull the time variable and its units (usually "days since YYYY-MM-DD")
  time_vals <- ncvar_get(nc_conn, "time")
  time_units <- ncatt_get(nc_conn, "time", "units")$value
  # Convert to human-readable date
  origin_date <- sub("days since ", "", time_units)
  file_date <- as.Date(time_vals, origin = origin_date)
  nc_close(nc_conn)
  file_date
})

# Flatten the list to a vector and pair with layer indices
layer_date_map <- data.frame(
  Layer_Index = 1:nlayers(trmm_stack_flip),
  Date = do.call(c, trmm_dates),
  Source_File = basename(trmm_files)
)

# Print the first 10 rows to check order
head(layer_date_map, 10)

Step 2: Cross-Check with Filenames (If Applicable) #

If your TRMM filenames include dates (e.g., trmm_20010101.nc), extract those dates too to confirm consistency:

# Adjust the format string to match your filename pattern
file_name_dates <- as.Date(basename(trmm_files), format = "trmm_%Y%m%d.nc")

# Check if dates from files match dates from NetCDF metadata
all.equal(do.call(c, trmm_dates), file_name_dates)

If this returns TRUE, your stack is correctly ordered by date.

2. Validate Point Extraction Accuracy #

There are a couple of quick checks to ensure your points are being mapped correctly:

Fix a Common Coordinate Order Issue #

First, a critical note: the extract() function expects coordinates in (longitude, latitude) order, but your xym matrix is cbind(lat, lon). This will shift your points to the wrong location! Fix it with:

# Correct coordinate order: lon first, lat second
xym <- cbind(lon, lat)

Visualize Points on the Raster #

Plot the first layer of your flipped stack and overlay your points to confirm they fall within your study area:

# Plot the first precipitation layer
plot(trmm_stack_flip[[1]], main = "TRMM HQprecipitation (First Layer)")
# Add your points in red for visibility
points(xym, pch = 19, col = "red", cex = 1.5)

If points are outside the raster extent, extract() will return NA for those locations—adjust your coordinates if needed.

Compare Direct Raster Values to Extracted Values #

Pull the value of a single point directly from the raster and compare it to your extracted data to confirm accuracy:

# Pick the first point (lon=91.85, lat=26.0166666)
target_cell <- cellFromXY(trmm_stack_flip[[1]], xym[1, ])
direct_raster_value <- trmm_stack_flip[[1]][target_cell]
extracted_value <- data2[1, 1]

# Print comparison
cat("Direct raster value:", direct_raster_value, "\n")
cat("Extracted value:", extracted_value, "\n")
cat("Values match?", direct_raster_value == extracted_value, "\n")

3. Validate Temporal Consistency of Extracted Data #

Once you confirm the stack order and extraction accuracy, check if the time series makes sense for your region:

Convert Extracted Data to a Time Series Data Frame #

# Convert the extracted matrix to a data frame with date column names
extracted_df <- as.data.frame(data2)
colnames(extracted_df) <- do.call(c, trmm_dates)

# Transpose to have dates as rows and points as columns
ts_df <- t(extracted_df)
ts_df <- cbind(Date = do.call(c, trmm_dates), as.data.frame(ts_df))
colnames(ts_df)[2:ncol(ts_df)] <- paste0("Point_", 1:nrow(xym))

Plot the Time Series #

Use ggplot2 to visualize precipitation trends—you should see logical seasonal patterns (e.g., monsoon peaks for your Southeast Asian coordinates):

library(ggplot2)

ggplot(ts_df, aes(x = Date, y = Point_1)) +
  geom_line(color = "darkblue") +
  labs(title = "TRMM HQprecipitation Time Series (Point 1)",
       y = "Precipitation (mm)", x = "Date") +
  theme_minimal()

If the plot shows unexpected gaps or erratic jumps, double-check your stack order or coordinate corrections.

Bonus: Confirm Your Raster Flip/Transpose Was Correct #

TRMM data often comes with reversed latitude or rotated grids. Verify your flip/transpose matches the original NetCDF metadata:

# Pull original lat/lon from the first file
nc_conn <- nc_open(trmm_files[1])
original_lat <- ncvar_get(nc_conn, "lat")
original_lon <- ncvar_get(nc_conn, "lon")
nc_close(nc_conn)

# Check if original lat is descending (common in TRMM)
if (original_lat[1] > original_lat[length(original_lat)]) {
  cat("Original latitudes are descending—your flip(y) was correct!\n")
}

# Check if flipped/transposed raster matches original coordinates
raster_lats <- yFromRow(trmm_stack_flip, 1:nrow(trmm_stack_flip))
raster_lons <- xFromCol(trmm_stack_flip, 1:ncol(trmm_stack_flip))

cat("Latitudes match after flip?", all.equal(sort(original_lat), raster_lats), "\n")
cat("Longitudes match after transpose?", all.equal(original_lon, raster_lons), "\n")

内容的提问来源于stack exchange，提问作者Jedd