Цикл через год, месяц, день для netCDF имя файла

Question

Я использую цикл for в R, чтобы читать файл netCDF из папки и извлекать значения для данного списка долготы, широты. Это похоже на работу, кроме ОДНОЙ ПРОБЛЕМЫ. Когда цикл возвращает значения по дате, он создает 29-29 января после 28 февраля. Я хочу, как обычно, 01 марта после 28 или 29 февраля (для високосного года). Вот мой R код:

# given latitude, longitude list 
sb1 <- data.frame(longitude=1:10,latitude =1:10) 

# Extracting zonal or sub-basin average rainfall from netCDF file 

sb1_r <- c() 
date <- c() 
rain_month <- c() 
rain_year <- c() 


for (year in 1998:1998){ 
    for (month in 1:3){ 
    for (day in seq_along(1:31)){ 
     FileName <- paste('3B42_daily',year,sprintf("%02d",month),sprintf("%02d", day),'7.SUB.nc', sep='.') 
    if (!file.exists(FileName)){ 
    next 
    } else { 

     File <- nc_open(FileName) 
     rain <- ncvar_get(File, 'r') 

     sb1_r[day] <- mean(apply(sb1,1,function(x)rain[x[1],x[2]]),na.rm = TRUE) 
     date[day] <- paste(year,sprintf("%02d", month),sprintf("%02d", day),sep='-') 
     rain_month <- data.frame(date,sb1_r) 
     nc_close(File) 
     } 
    } 

    rain_year <- rbind(rain_year,rain_month) 
    } 

} 

Вы можете найти ежедневные данные NetCDF в течение трех месяцев по этой ссылке: https://drive.google.com/open?id=0B8rqKaYt0VEaMWVGc1gzdXI1U28

Answer 1

Вместо того, чтобы пытаться создать имена файлов делать противоположное. Извлеките имена файлов и для каждого файла получите дату из имени файла и sb1_r из файла. Вы можете сделать это быстрее, используя rbindlist из пакета data.table (но не обязательно).

# дано широта, долгота список SB1 < - data.frame (долгота = 1: 10, широта = 1: 10)

# Extracting zonal or sub-basin average rainfall from netCDF file 
filenames = list.files(path = ".", pattern = ".nc") 
rain_year = data.frame() 

require(ncdf4) 
for(FileName in filenames){ 
    File <- nc_open(FileName) 
    # Create Date 
    year <- strsplit(FileName, split = '[.]')[[1]][2] 
    month <- strsplit(FileName, split = '[.]')[[1]][3] 
    day <- strsplit(FileName, split = '[.]')[[1]][4] 
    date = paste(year, month, day, sep = "-") 
    # get value 
    rain <- ncvar_get(File, 'r') 
    sb1_r <- mean(apply(sb1,1,function(x)rain[x[1],x[2]]),na.rm = TRUE) 
    # update data.frame 
    rain_year = rbind(rain_year, data.frame(date = date, sb1_r = sb1_r)) 
    nc_close(File) 
} 

# Faster using data.table package 
require(data.table) 
temp = rbindlist(
    lapply(X = filenames, function(FileName){ 
    year <- as.integer(strsplit(FileName, split = '[.]')[[1]][2]) 
    month <- as.integer(strsplit(FileName, split = '[.]')[[1]][3]) 
    day <- as.integer(strsplit(FileName, split = '[.]')[[1]][4]) 
    date = paste(year, month, day, sep = "-") 
    File <- nc_open(FileName) 
    rain <- ncvar_get(File, 'r') 
    sb1_r <- mean(apply(sb1,1,function(x)rain[x[1],x[2]]),na.rm = TRUE) 
    return(data.frame(date = date, sb1_r = sb1_r)) 
    }) 
) 

Answer 2

Обратите внимание, что вышеуказанные коды в R являются НЕ правильно, если файл netcdf rainfall использует сетку с экви-зоной, что редко бывает! (И это не относится к файлам TRMM, используемым в этом примере). Это обычная ошибка при обработке данных с привязкой к сетке.

Например, если у вас есть обычная сетка lat/lon, вам необходимо учитывать уменьшение косинуса измерения долготы при движении к полюсам. Ошибка небольшая, если ваш суббассейн мал, но может стать значительным, если область большая. Для некоторых видов сеток, например. уменьшенные гауссовские решетки, ошибка может быть очень значительной, если ваш подобъект так совпаден с прерывистым изменением чисел долготы, особенно для «негладких» полей, таких как осадки.

Я бы всегда выполнял обработку зон и суббассейнов с помощью CDO, чтобы убедиться, что расчет выполнен правильно. Если вы используете CDO, то средние области и зональные функции усреднения учитывают собственную сетку.

Таким образом, мой код будет выглядеть примерно так:

#!/bin/bash 

# define the lat-lon bounds of your sub area 
lat1=20 
lat2=30 
lon1=0 
lon2=20 

# merge all the daily files into one file 
# do this one month at a time as some system limit number of open files 

year=1998 # can make this a loop if you want multiple years 
for month in `seq -f "%02g" 1 3` ; do 
    files=`ls 3B42_daily1998${month}*.nc` 
    cdo mergetime $files TRMM_${month}.nc 
done 
cdo mergetime $TRMM_*.nc TRMM_timeseries.nc 

# now extract the subdomain 
cdo sellonlatbox,$lon1,$lon2,$lat1,$lat2 TRMM_timeseries.nc TRMM_box.nc 

# CORRECT area average 
cdo fldmean TRMM_box.nc TRMM_box_av.nc 

# zonal average 
cdo zonmean TRMM_box.nc TRMM_box_zon.nc 

Answer 3

#!/usr/bin/env Rscript 
argv<-commandArgs(trailingOnly=TRUE) 
if(length(argv)==2 & argv[1] <= argv[2]) { 
    if (is.na(strptime(sprintf("%s",argv[1]),"%Y%m%d"))) { 
    cat("arguments valid check error: ", argv[1], "\n") 
    stop() 
    } 
    if (argv[2]!=format(strptime(sprintf("%s",argv[2]),"%Y%m%d"),"%Y%m%d")) { 
    cat("arguments valid check error: ", argv[2], "\n") 
    stop() 
    } 
} else if (length(argv)==2 & argv[1] > argv[2]) { 
    print(sprintf("error: %s is greater than %s",argv[1],argv[2])) 
    stop() 
} else if (length(argv)!=2) { 
    script.name<-basename(strsplit(commandArgs(trailingOnly=FALSE)[4],"=")[[1]][2]) 
    print(sprintf("Usage: %s startDate endDate",script.name)) 
    stop() 
} 

filelist<-c() 
for (Ymd in format(seq(
    as.POSIXct(sprintf("%s",argv[1]),format="%Y%m%d"), 
    as.POSIXct(sprintf("%s",argv[2]),format="%Y%m%d"), 
    by="24 hour"),"%Y%m%d")) { 
    filelist<-append(filelist, sprintf("%s.%s.%s.%s.%s","3B42_daily",substr(Ymd,1,4),substr(Ymd,5,6),substr(Ymd,7,8),"7.SUB.nc")) 
} 

sb1_r <- c() 
date <- c() 
rain_month <- c() 
rain_year <- c() 

for (i in 1:length(filelist)) { 
if (!file.exists(filelist[i])){ 
    next 
} else { 
    year <- as.numeric(substr(filelist[i],12,15)) 
    month <- as.numeric(substr(filelist[i],17,18)) 
    day <- as.numeric(substr(filelist[i],20,21)) 
    File <- nc_open(filelist[i]) 
    rain <- ncvar_get(File, 'r') 

    sb1_r[day] <- mean(apply(sb1,1,function(x)rain[x[1],x[2]]),na.rm = TRUE) 
    date[day] <- paste(year,sprintf("%02d", month),sprintf("%02d", day),sep='-') 
    rain_month <- data.frame(date,sb1_r) 
    nc_close(File) 
    } 
}