Schleife durch Jahr, Monat, Tag für NetCDF-Dateiname

Question

Ich verwende eine for-Schleife in R, um NetCDF-Datei aus einem Ordner zu lesen und Werte für gegebene Länge, Breite und Breite zu extrahieren. Es sieht so aus als würde es funktionieren, außer EIN PROBLEM. Wenn die Schleife Werte gegen das Datum zurückgibt, erstellt sie den 29. bis 31. Januar nach dem 28. Februar. Ich möchte wie üblich den 1. März nach dem 28. oder 29. Februar (für das Schaltjahr). Hier ist mein R-Code:

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

} 

Sie können täglich netCDF-Daten für drei Monate auf diesen Link: https://drive.google.com/open?id=0B8rqKaYt0VEaMWVGc1gzdXI1U28

Answer 1

Anstatt zu versuchen, die Dateinamen zu tun, das Gegenteil zu erstellen. Extrahiere die Dateinamen und erhalte für jede Datei das Datum aus dem Dateinamen und sb1_r aus der Datei. Sie können das schneller mit rbindlist aus dem data.table-Paket tun (wird aber nicht benötigt).

# gegebenen Breitengrad, Längengrad Liste sb1 < - data.frame (Länge = 1: 10, Breite = 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

Beachten Sie, dass die oben genannten Codes in R sind NICHT richtig es sei denn, Die Niederschlags-Netcdf-Datei verwendet ein Äqui-Area-Raster, was selten der Fall ist! (Dies ist bei den in diesem Beispiel verwendeten TRMM-Dateien nicht der Fall). Dies ist ein häufiger Fehler beim Umgang mit gerasterten Daten.

Wenn Sie beispielsweise ein reguläres Lat/Lon-Gitter haben, müssen Sie die Cosinus-Reduktion der Längengrad-Dimension berücksichtigen, während Sie sich in Richtung der Pole bewegen. Der Fehler ist klein, wenn Ihr Sub-Becken klein ist, kann aber signifikant werden, wenn die Fläche groß ist. Für einige Arten von Gittern, z.B. Reduzierte Gauß'sche Gitter, der Fehler kann sehr signifikant sein, wenn Ihre Subdomain gerade so zufällig mit der diskontinuierlichen Änderung der Längenpunktzahlen zusammenfällt, besonders für "nicht glatte" Felder wie Niederschlag.

Ich würde immer zonale und Unterbecken Verarbeitung mit CDO durchführen, um sicherzustellen, dass die Berechnung korrekt durchgeführt wird. Wenn Sie CDO verwenden, berücksichtigen die Mittelwerte für den Bereich und die zonalen Mittelwerte das native Raster.

So würde mein Code wie folgt aussehen:

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