Автоматизировать ncvar_get-чтение различных имен переменных в файлах NetCDF

У меня есть набор данных NetCDF с двумя климатическими сценариями (rcp и hist), каждый из которых содержит 25 файлов. Каждый файл содержит данные для переменной «pr», «tas», «tasmax» или «tasmin». Я написал цикл for для многократного чтения файлов hist и rcp, читал их с помощью nc_open, извлекал переменную с помощью ncvar_get и, наконец, выполнял вычисления в форме mean(abs(hist - rcp)), чтобы получить среднее абсолютное расстояние между каждой парой of hist и rcp.Проблема: поскольку ncvar_get требует точного имени переменной текущего файла, я написал блок if else (см. ниже), который должен найти имя переменной текущего файла и применить его для ncvar_get.Запуск кода, который я получаю следующая ошибка:

[1] "vobjtovarid4: error #F: I could not find the requsted var (or dimvar) in the file!"

[1] "var (or dimvar) name: tas"

[1] "file name: /data/historical/tasmax_ICHEC-EC-EARTH_DMI-HIRHAM5_r3i1p1.nc" Error in vobjtovarid4(nc, varid, verbose = verbose, allowdimvar = TRUE) : Variable not found



 #Extract of the files in the hist list. Same file names in the rcp list, but different directory

    > hist.files.cl <- list.files("/historical", full.names = TRUE)

    > hist.files.cl

 [1] "/historical/pr_CNRM-CERFACS-CNRM-CM5_ALADIN53_r1i1p1.nc"           
 [2] "/historical/pr_CNRM-CERFACS-CNRM-CM5_ALARO-0_r1i1p1.nc"            
 [3] "/historical/pr_ICHEC-EC-EARTH_HIRHAM5_r3i1p1.nc"                   
 [4] "/historical/pr_ICHEC-EC-EARTH_RACMO22E_r12i1p1.nc"                 
 [5] "/historical/pr_ICHEC-EC-EARTH_RCA4_r12i1p1.nc"                     
 [6] "/historical/pr_MPI-M-MPI-ESM-LR_RCA4_r1i1p1.nc"                    
 [7] "/historical/pr_MPI-M-MPI-ESM-LR_REMO2009_r1i1p1.nc"                
 [8] "/historical/pr_MPI-M-MPI-ESM-LR_REMO2009_r2i1p1.nc"                
 [9] "/historical/tas_CNRM-CERFACS-CNRM-CM5_CNRM-ALADIN53_r1i1p1.nc"     
[10] "/historical/tas_CNRM-CERFACS-CNRM-CM5_RMIB-UGent-ALARO-0_r1i1p1.nc"
[11] "/historical/tas_ICHEC-EC-EARTH_DMI-HIRHAM5_r3i1p1.nc"              
[12] "/historical/tas_ICHEC-EC-EARTH_KNMI-RACMO22E_r12i1p1.nc"           
[13] "/historical/tas_ICHEC-EC-EARTH_SMHI-RCA4_r12i1p1.nc"               
[14] "/historical/tas_MPI-M-MPI-ESM-LR_MPI-CSC-REMO2009_r1i1p1.nc"       
[15] "/historical/tas_MPI-M-MPI-ESM-LR_MPI-CSC-REMO2009_r2i1p1.nc"       
[16] "/historical/tasmax_ICHEC-EC-EARTH_DMI-HIRHAM5_r3i1p1.nc"           
[17] "/historical/tasmax_ICHEC-EC-EARTH_KNMI-RACMO22E_r12i1p1.nc"        
[18] "/historical/tasmax_ICHEC-EC-EARTH_SMHI-RCA4_r12i1p1.nc"        


euc.distance <- list()


for(i in 1:length(hist.files.cl)) {

#Open ith file in list of hist files as well as in list of rcp files

  hist.data <- nc_open(hist.files.cl[i])   
  rcp.data <- nc_open(rcp.files.cl[i])

  if(grepl("pr", hist.data$filename)){
    hist.var <- ncvar_get(hist.data, "pr")
    rcp.var <- ncvar_get(rcp.data, "pr")
    }else if (grepl("tas", hist.data$filename)){
    hist.var <- ncvar_get(hist.data, "tas")
    rcp.var <- ncvar_get(rcp.data, "tas")
    }else if (grepl("tasmax", hist.data$filename)){
    hist.var <- ncvar_get(hist.data, "tasmax")
    rcp.var <- ncvar_get(rcp.data, "tasmax")
    }else{
    hist.var <- ncvar_get(hist.data, "tasmin")
    rcp.var <- ncvar_get(rcp.data, "tasmin")
    }
 #Converting temperature variable from K to °C: 

  if(grepl("tas", hist.data$filename)){

    hist.var <- hist.var-273.15
    rcp.var <- rcp.var-273.15
  }

 #Find for the ith rcp file with dim=(1,1,360) in the ith hist file with dim=(385,373,360) the grid point with the best fitting distribution (each grid point consists of a distribution of 360 time steps).The calculation may contain errors...

  euc.distance[[i]] <- apply(hist.var, c(1,2), function(x) mean(abs(rcp.var - x)))
  min_values <-  which(rank(euc.distance[[i]], ties.method='min') <= 10)
}

Так как кат указал на возможную причину ошибки, но предложенный подход к извлечению интересующей части (=имя переменной) из имени файла не работает. Раньше я пытался автоматизировать извлечение имени переменной с помощью stringr("filename",startposition, endposition), пока не заметил, что в этом нет смысла, потому что каждое имя переменной (pr, tas, tasmax, tasmin) имеет другую строку длина. Какие возможности вы видите для меня? Спасибо большое!


r if-statement automation netcdf4
person Mike    schedule 02.03.2018    source источник
comment
на самом деле tas включен в tasmax, и именно поэтому вы сталкиваетесь с проблемой. Вам придется быть более точным в своем grepl звонке. Вероятно, вы можете извлечь именно ту часть, которая представляет интерес, а затем выполнить только один вызов (поскольку эта часть является именно тем именем переменной, которое вам нужно). Что-то вроде keypart <- sub("regex", "what I need", hist.data$filename) ; hist.var <- ncvar_get(hist.data, keypart) ; rcp.var <- ncvar_get(rcp.data, keypar): один блок вместо 4...   -  person Cath    schedule 02.03.2018
comment
Если вы предоставите нам дополнительную информацию (см. Как спросить и минимальный воспроизводимый пример) нам будет намного проще помочь   -  person Cath    schedule 02.03.2018
comment
Да, ваше объяснение имеет смысл! Я еще пытался написать общую функцию для автоматического извлечения интересующей части для каждого файла. Я использовал, например. stringr(filename, 55, 57), но я заметил, что в моем случае это не имеет смысла, так как из-за разной длины имен переменных pr имеет другие индексы в строке, чем tas, tasmax или tasmin. Я отредактирую свою вышеуказанную проблему и предоставлю имена файлов и т. Д. Спасибо!   -  person Mike    schedule 02.03.2018
comment
в вашем случае keypart <- sub("^([a-z]+)_.+", "\\1", basename(hist.files.cl)) (вы получите вектор всех ключевых частей, соответствующих именам файлов)   -  person Cath    schedule 02.03.2018
comment
Отлично, это работает! Большое спасибо Кэт!   -  person Mike    schedule 02.03.2018

Ответы (1)


arrow_upward
3
arrow_downward

Чтобы немного дополнить мой комментарий, если вам нужно работать с каждым файлом, вы можете сделать это сразу, сложив все в список.

Итак, сначала получите «ключевую часть» для каждого файла:

keyparts <- sub("^([a-z]+)_.+", "\\1", basename(hist.files.cl))
keyparts
# [1] "pr"     "pr"     "pr"     "pr"     "pr"     "pr"     "pr"     "pr"    
# [9] "tas"    "tas"    "tas"    "tas"    "tas"    "tas"    "tas"    "tasmax"
#[17] "tasmax" "tasmax"

Затем вы можете использовать lapply, чтобы делать то, что вам нужно сделать для всех файлов одновременно:

my_res <- lapply(seq(keyparts), 
                 function(i){

         hist.data <- nc_open(hist.files.cl[i])   
         rcp.data <- nc_open(rcp.files.cl[i])

         hist.var <- ncvar_get(hist.data, keyparts[i])
         rcp.var <- ncvar_get(rcp.data, keyparts[i])

         if(keyparts[i]=="tas"){
           hist.var <- hist.var-273.15
           rcp.var <- rcp.var-273.15
         }


        euc.distance <- apply(hist.var, c(1,2), function(x) mean(abs(rcp.var - x)))
        min_values <-  which(rank(euc.distance[[i]], ties.method='min') <= 10)

        return(list(euc.distance=euc.distance, min.values=min.values))

                   })
person Cath    schedule 03.03.2018