基于播种日期提取格点降水特定时段数据的技术求助

I’ll help you modify your code to pull fixed-length precipitation periods after sowing, handling both same-season and cross-year cases (like Rabi crops sown in late November/December). Here’s a practical solution aligned with your original code’s logic:

Step 1: Setup Dependencies and Parameters #

First, install and load the required packages:

install.packages(c("ncdf4", "abind"))
library(ncdf4)
library(abind)

Define your core parameters (adjust these to match your data and desired windows):

# Target windows: start days after sowing, fixed window length
window_starts_after_sowing <- c(30, 80)  # Example: 30th day post-sowing, 80th day post-sowing
window_length <- 30                      # Each window is exactly 30 days long
nyear <- 30                              # Number of years (matches your original 2:30 loop)
base_year <- 1980                        # Base year for your dataset

# Load sowing/harvest date arrays (replace with your actual file loading code)
# nc_sowing <- nc_open("path/to/sowing_dates.nc")
# sdate <- ncvar_get(nc_sowing, "sowing_date")  # [lon, lat, nyear]
# nc_close(nc_sowing)
# nc_harvest <- nc_open("path/to/harvest_dates.nc")
# hdate <- ncvar_get(nc_harvest, "harvest_date")  # [lon, lat, nyear]
# nc_close(nc_harvest)

# Get spatial dimensions from first precipitation file
first_prec_file <- paste0("E:/Paper_2018/Prec/prec_", base_year + 2, ".nc")
nc <- nc_open(first_prec_file)
lons <- ncvar_get(nc, "longitude")
lats <- ncvar_get(nc, "latitude")
nc_close(nc)
nlon <- length(lons)
nlat <- length(lats)
n_windows <- length(window_starts_after_sowing)

# Initialize output array: [lon, lat, number_of_windows, window_length, number_of_years]
# Split into separate arrays if you don't need all windows/years in one structure
precipitation_windows <- array(0, dim = c(nlon, nlat, n_windows, window_length, nyear))

Step 2: Modified Loop to Extract Windows #

This loop handles both Kharif (same-year sowing/harvest) and Rabi (cross-year) cases, ensuring every window stays exactly 30 days long:

for (year_idx in 2:nyear) {
  current_harvest_year <- base_year + year_idx
  prev_year <- current_harvest_year - 1
  cat("Processing harvest year:", current_harvest_year, "\n")
  
  # Load precipitation for current and previous years
  current_prec_file <- paste0("E:/Paper_2018/Prec/prec_", current_harvest_year, ".nc")
  nc_current <- nc_open(current_prec_file)
  prec_current <- ncvar_get(nc_current, "pr")  # [lon, lat, 365]
  nc_close(nc_current)
  
  prev_prec_file <- paste0("E:/Paper_2018/Prec/prec_", prev_year, ".nc")
  nc_prev <- nc_open(prev_prec_file)
  prec_prev <- ncvar_get(nc_prev, "pr")        # [lon, lat, 365]
  nc_close(nc_prev)
  
  # Get sowing/harvest dates for current harvest year
  sdate_curr_harvest <- sdate[, , year_idx]
  hdate_curr_harvest <- hdate[, , year_idx]
  
  # Determine sowing year/day for each grid point
  kharif_mask <- hdate_curr_harvest > sdate_curr_harvest  # Same-year sowing
  rabi_mask <- !kharif_mask                               # Previous-year sowing
  
  sowing_year <- array(current_harvest_year, dim = c(nlon, nlat))
  sowing_day <- sdate_curr_harvest
  # Update Rabi points to use previous year's sowing date
  sowing_year[rabi_mask] <- prev_year
  sowing_day[rabi_mask] <- sdate[, , year_idx - 1]
  
  # Process each target window
  for (win_idx in 1:n_windows) {
    win_start_offset <- window_starts_after_sowing[win_idx]
    win_start_abs <- sowing_day + win_start_offset
    win_end_abs <- win_start_abs + window_length - 1
    
    # Case 1: Entire window is within sowing year
    same_year_mask <- win_end_abs <= 365
    if (any(same_year_mask)) {
      precipitation_windows[same_year_mask, win_idx, , year_idx] <- 
        prec_current[same_year_mask, win_start_abs[same_year_mask]:win_end_abs[same_year_mask]]
    }
    
    # Case 2: Window crosses into the next year
    cross_year_mask <- win_end_abs > 365
    if (any(cross_year_mask)) {
      days_in_sowing_year <- 365 - win_start_abs[cross_year_mask] + 1
      days_in_next_year <- window_length - days_in_sowing_year
      
      # Split cross-year points into Kharif/Rabi
      cross_kharif <- cross_year_mask & kharif_mask
      cross_rabi <- cross_year_mask & rabi_mask
      
      # Handle Kharif windows spilling into next harvest year
      if (any(cross_kharif)) {
        next_year <- current_harvest_year + 1
        next_prec_file <- paste0("E:/Paper_2018/Prec/prec_", next_year, ".nc")
        nc_next <- nc_open(next_prec_file)
        prec_next <- ncvar_get(nc_next, "pr")
        nc_close(nc_next)
        
        kharif_prev_part <- prec_current[cross_kharif, win_start_abs[cross_kharif]:365]
        kharif_next_part <- prec_next[cross_kharif, 1:days_in_next_year[cross_kharif]]
        precipitation_windows[cross_kharif, win_idx, , year_idx] <- abind(kharif_prev_part, kharif_next_part, along = 3)
      }
      
      # Handle Rabi windows spanning sowing year to harvest year
      if (any(cross_rabi)) {
        rabi_prev_part <- prec_prev[cross_rabi, win_start_abs[cross_rabi]:365]
        rabi_next_part <- prec_current[cross_rabi, 1:days_in_next_year[cross_rabi]]
        precipitation_windows[cross_rabi, win_idx, , year_idx] <- abind(rabi_prev_part, rabi_next_part, along = 3)
      }
    }
  }
}

Key Improvements Explained #

1. Sowing Year Logic: Uses your original hdate check to distinguish Kharif/Rabi crops, ensuring we pull sowing dates from the correct year.
2. Fixed-Length Guarantees: Splits cross-year windows into two parts and combines them to maintain exactly 30 days, no matter the sowing date.
3. Grid-Wide Processing: Handles all grid points at once instead of a single (x,y) point, making the code scalable for large datasets.
4. Edge Case Handling: Explicitly manages Kharif windows that spill into the next harvest year and Rabi windows that span two calendar years.

Notes #

• Replace the placeholder sowing/harvest date loading code with your actual file-reading logic.
• If you don’t need all years in the output, adjust the array dimensions (e.g., average over years or save each year separately).
• For very large datasets, use parallel processing (like the foreach package) to speed up the loop.

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