Boucle passant par l'année, le mois, le jour pour le nom de fichier netCDF

Question

J'utilise une boucle for dans R pour lire le fichier netCDF d'un dossier et extraire des valeurs pour une liste donnée de longitude, latitude. Cela ressemble à travailler, sauf UN PROBLÈME. Lorsque la boucle renvoie des valeurs par rapport à la date, elle crée du 29 au 31 janvier après le 28 février. Je veux, comme d'habitude, le 01 mars après le 28 ou 29 février (pour les années bissextiles). Voici mon code 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) 
    } 

} 

Vous pouvez trouver des données netCDF par jour pendant trois mois à ce lien: https://drive.google.com/open?id=0B8rqKaYt0VEaMWVGc1gzdXI1U28

Answer 1

Au lieu d'essayer de créer les noms de fichiers font le contraire. Extraire les noms de fichiers et obtenir pour chaque fichier la date du fichier et sb1_r du fichier. Vous pouvez le faire plus rapidement en utilisant rbindlist à partir du paquet data.table (mais ce n'est pas nécessaire).

# donné la latitude, la longitude de la liste SB1 < - data.frame (longitude = 1: 10, latitude = 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

Notez que les codes ci-dessus dans R sont pas correcte à moins le fichier netcdf des précipitations utilise une grille équi-zone, ce qui est rarement le cas! (Et ce n'est pas le cas pour les fichiers TRMM utilisés dans cet exemple). Ceci est une erreur courante lors de la manipulation de données maillées. Par exemple, si vous avez une grille latitude/longitude régulière, vous devez prendre en compte la réduction en cosinus de la cote de longitude lorsque vous vous déplacez vers les pôles. L'erreur est faible si votre sous-bassin est petit, mais peut devenir significatif si la zone est grande. Pour certains types de grilles, par ex. Si les grilles gaussiennes sont réduites, l'erreur peut être très significative si votre sous-domaine coïncide avec le changement discontinu des nombres de points de longitude, en particulier pour les champs "non lisses" tels que les précipitations.

J'effectuerais toujours un traitement zonal et sous-bassin en utilisant CDO pour m'assurer que le calcul est effectué correctement. Si vous utilisez CDO, les moyennes de zone et les fonctions de moyennage par zone représentent la grille native.

Ainsi mon code ressemblerait à quelque chose comme ceci:

#!/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) 
    } 
}