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

Я использую цикл for в R для чтения файла netCDF из папки и извлечения значений для заданного списка долготы и широты. Вроде работает, кроме ОДНОЙ ПРОБЛЕМЫ. Когда цикл возвращает значения по дате, он создает с 29 по 31 января после 28 февраля. Я хочу, как обычно, 1 марта после 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


r loops netcdf cdo-climate
person Wahid Palash    schedule 01.08.2016    source источник
comment
У вас for (day in seq_along(1:31)) для января, февраля и марта. Но в феврале всего 28 дней. Может в этом проблема? Если да, вам нужно настроить цикл.   -  person dataanalyst    schedule 01.08.2016
comment
@Gandalf Но у меня нет файлов NetCDF с именем 3B42_daily.1998.02.29.7.SUB и так далее. Чтобы избежать этого, я добавляю if (! File.exists (FileName)) {в свой код.   -  person Wahid Palash    schedule 01.08.2016
comment
Просто чтобы указать, что использование функции mean не даст вам правильного ответа при использовании, например, обычная сетка широты и долготы, так как ячейки сетки различаются по размеру. Таким образом, значение в каждой ячейке должно быть взвешено по площади ячейки. Намного лучше просто использовать CDO, который учитывает это автоматически - см. Ниже.   -  person Adrian Tompkins    schedule 28.01.2017

Ответы (3)


arrow_upward
1
arrow_downward

Вместо того, чтобы пытаться создать имена файлов, сделайте наоборот. Извлеките имена файлов и для каждого файла получите дату из имени файла и 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))
  })
)
person conighion    schedule 01.08.2016

arrow_upward
2
arrow_downward

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

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

Я бы всегда выполнял зональную и суб-бассейновую обработку с использованием 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
person Adrian Tompkins    schedule 18.01.2017

arrow_upward

arrow_downward

person    schedule
comment
Не могли бы вы пояснить, что вы изменили, по сравнению с другими вопросами? Сейчас сложно понять преимущества вашего решения. - person Malte Hartwig; 04.01.2018