Вложенные циклы for с NetCDF в R

Я очень-очень благодарен всем, кто хочет помочь мне с проблемой, с которой я действительно застрял. Но заранее: это сложная тема, и я изо всех сил стараюсь объяснить, что я собираюсь делать со своим кодом. Речь идет о климатических данных в файлах NetCDF, которые содержат месячные данные о температуре (tas) и осадках (pr) за периоды времени с 1971 по 2000 и с 2071 по 2100 годы. nc-файлы исторического периода содержат ок. 440x400 точек сетки (Карта Европы). nc-файлы будущего периода содержат точку сетки 1x1 (для интересующего города). Каждая точка сетки имеет 360 значений температуры или осадков (в зависимости от модели), по одному значению на каждый месяц из 30-летних периодов. Другими словами: каждая точка сетки имеет распределение 360 точек. Теперь я хочу итеративно вычислить статистическую разницу между распределением точек сетки одного города (2071–2100 гг.) и распределением точек сетки каждой Европы (1971–2000 гг.). Я получу одно усредненное абсолютное расстояние на точку сетки Европы. Идея состоит в том, чтобы найти в растре европейской сетки точку сетки, в которой распределение температуры или осадков наиболее похоже на распределение интересующего города в будущем. Я должен провести этот расчет для 30 различных моделей климата.

# List filenames of the directory

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

#Create array for desired ‘similarity indices’. One matrix per climate model run.

sim.array <- array(NA, dim = c(440,400,30))

#Looping through the models of the period 1971-2000. Some containing precipitation data others temperature (see if…else) 

for(k in 1:length(hist.files))   {
        hist.data <- nc_open(hist.files[k])   

   if(grepl("pr", hist.data$filename)){
    hist.tas <- ncvar_get(hist.data, "pr")
        }else{
    hist.tas <- ncvar_get(hist.data, "tas") 
    hist.tas <- kelvin.to.celsius(hist.tas, round=2)
   }

#Looping through the models of the 2071 to 2100 period (city). Some containing precipitation data others temperature (see if…else)

for(r in 1:length(rcp.files)) {
    rcp.data <- nc_open(rcp.files[r])
    if(grepl("pr", rcp.data$filename)){
    rcp.tas <- ncvar_get(rcp.data, "pr") 
        }else{
    rcp.tas <- ncvar_get(rcp.data, "tas")
    rcp.tas <- kelvin.to.celsius(rcp.tas, round=2)
        }

#This if statement because hist contains more models than rcp and I want to exclusively use the models contained in both of them.  

if(hist.data %in% rcp.data) {  

#Looping through the grid points of ‘hist’ model k. Lastly the function that calculates for each grid point of the model a difference value (always to the one grid point of ‘rcp’). My idea of the break statement was to loop nrow and ncol the same times, but I’m not sure if break does what I intended to.       

for(i in 1:nrow(hist.tas)) { 
       for(j in 1:ncol(hist.tas)) {
    sim.array[i,j,k] <- abs(sum(rcp.tas - hist.tas[i,j,])/360)
break
    }
  print(sim.array[i,j,k])
  }
 }
}   
}
sim.array[1,1,1]

Ну, я получаю массив, полный NA. Сообщение об ошибке отсутствует, но что-то идет не так! Кто может найти ошибку? Я ценю любую помощь. Заранее большое спасибо!

Обновление: ваши предложения кажутся разумным решением! До сих пор у меня не было времени их применить, но я сделаю это позже! Я думал о векторизации, но мне не удалось сделать векторы из трехмерных массивов без грязного кода, полного разных векторов в конце ... Я также не знал, как удалить модели, которые не соответствуют hist и rcp. С пересечением и %in% я знал индекс несоответствующих файлов ... но должен быть лучший способ, чем отмечать вручную все эти индексы для удаления, не так ли? Пожалуйста, взгляните на некоторые имена файлов истории:

> hist.files.tas <- list.files("/historical", full.names = TRUE, pattern = "tas")
> hist.files.tas
 [1] "/historical/tas_CNRM-CERFACS-CNRM-CM5_CLMcom-CCLM4-8-17_r1i1p1.nc"   
 [2] "/historical/tas_CNRM-CERFACS-CNRM-CM5_CNRM-ALADIN53_r1i1p1.nc"       
 [3] "/historical/tas_CNRM-CERFACS-CNRM-CM5_RMIB-UGent-ALARO-0_r1i1p1.nc"  
 [4] "/historical/tas_CNRM-CERFACS-CNRM-CM5_SMHI-RCA4_r1i1p1.nc"           
 [5] "/historical/tas_ICHEC-EC-EARTH_CLMcom-CCLM4-8-17_r12i1p1.nc"         
 [6] "/historical/tas_ICHEC-EC-EARTH_DMI-HIRHAM5_r3i1p1.nc"                
 [7] "/historical/tas_ICHEC-EC-EARTH_KNMI-RACMO22E_r12i1p1.nc"             
 [8] "/historical/tas_ICHEC-EC-EARTH_KNMI-RACMO22E_r1i1p1.nc"              
 [9] "/historical/tas_ICHEC-EC-EARTH_SMHI-RCA4_r12i1p1.nc"                 
[10] "/historical/tas_IPSL-IPSL-CM5A-MR_INERIS-WRF331F_r1i1p1.nc"          
[11] "/historical/tas_IPSL-IPSL-CM5A-MR_SMHI-RCA4_r1i1p1.nc"               
[12] "/historical/tas_MOHC-HadGEM2-ES_CLMcom-CCLM4-8-17_r1i1p1.nc"         
[13] "/historical/tas_MOHC-HadGEM2-ES_KNMI-RACMO22E_r1i1p1.nc"             
[14] "/historical/tas_MOHC-HadGEM2-ES_SMHI-RCA4_r1i1p1.nc"   

Есть еще модели с переменными tasmax и tasmin. Всего в hist 71 файл, а в rcp только 30. Не могли бы вы привести пример того, как написать автоматический код для удаления несовпадающих файлов hist? Спасибо большое!


person Mike    schedule 19.02.2018    source источник
comment
Вы просматривали 1 отдельный файл (summary, plot)? Есть ли NA (-32767 или -9999)? sum(..., na.rm = TRUE) поможет?   -  person Tung    schedule 20.02.2018
comment
Да, NA, к сожалению, не являются причиной ошибки... Спасибо!   -  person Mike    schedule 20.02.2018


Ответы (2)


Мне кажется, что приведенное ниже не имеет смысла и всегда ЛОЖЬ:

if (hist.data %in% rcp.data)

Так что ничего не происходит с sim_array

Я бы начал с чего-то вроде этого:

hist.files.pr <- list.files("/historical", full.names = TRUE, pattern="pr")
hist.files.tas <- list.files("/historical", full.names = TRUE, pattern="tas")
rcp.files.pr <- list.files("/rcp", full.names = TRUE, pattern="pr")
rcp.files.tas <- list.files("/rcp", full.names = TRUE, pattern="tas")

На этом этапе вы можете удалить файлы из «hist» для моделей, которых нет в «rcp».

hist.files.tas <- c( "/historical/tas_CNRM-CERFACS-CNRM-CM5_CLMcom-CCLM4-8-17_r1i1p1.nc", "/historical/tas_CNRM-CERFACS-CNRM-CM5_CNRM-ALADIN53_r1i1p1.nc", "/historical/tas_CNRM-CERFACS-CNRM-CM5_RMIB-UGent-ALARO-0_r1i1p1.nc", "/historical/tas_CNRM-CERFACS-CNRM-CM5_SMHI-RCA4_r1i1p1.nc", "/historical/tas_ICHEC-EC-EARTH_CLMcom-CCLM4-8-17_r12i1p1.nc", "/historical/tas_ICHEC-EC-EARTH_DMI-HIRHAM5_r3i1p1.nc", "/historical/tas_ICHEC-EC-EARTH_KNMI-RACMO22E_r12i1p1.nc", "/historical/tas_ICHEC-EC-EARTH_KNMI-RACMO22E_r1i1p1.nc", "/historical/tas_ICHEC-EC-EARTH_SMHI-RCA4_r12i1p1.nc", "/historical/tas_IPSL-IPSL-CM5A-MR_INERIS-WRF331F_r1i1p1.nc", "/historical/tas_IPSL-IPSL-CM5A-MR_SMHI-RCA4_r1i1p1.nc", "/historical/tas_MOHC-HadGEM2-ES_CLMcom-CCLM4-8-17_r1i1p1.nc", "/historical/tas_MOHC-HadGEM2-ES_KNMI-RACMO22E_r1i1p1.nc", "/historical/tas_MOHC-HadGEM2-ES_SMHI-RCA4_r1i1p1.nc")

# in this example, fut files is a subset of hist files; that should be OK if their filename structure is the same

rcp.files.tas <- hist.files.tas[1:7]

getModels <- function(ff) {
    base <- basename(ff)
    s <- strsplit(base, "_")
    sapply(s, function(i) i[[2]])
}

getHistModels <- function(hist, fut) {
    h <- getModels(hist)
    uh <- unique(h)
    uf <- unique(getModels(fut))
    uhf <- uh[uh %in% uf]
    hist[h %in% uhf]
}


hist.files.tas.selected <- getHistModels(hist.files.tas, rcp.files.tas)
# hist.files.pr.selected <- getHistModels(hist.files.pr, rcp.files.pr)

Двойного цикла (k, r), вероятно, можно было бы избежать, выполнив что-то вроде этого:

library(raster)
his.pr <- values(stack(hist.files.pr.selected, var="pr")))
his.tas <- values(stack(hist.files.tas.selected, var="tas"))
rcp.pr <- values(stack(hist.files.pr, var="pr"))
rcp.tas <- values(stack(hist.files.tas, var="tas"))

И цикла (i, j) по строкам и столбцам, вероятно, тоже можно избежать. R векторизован. То есть вы можете делать такие вещи, как (1:10) - 2.

В любом случае, ваш код очень сложно читать со всеми этими вложенными циклами. Если они вам действительно нужны, было бы лучше вызывать функции. Для получения дополнительной помощи предоставьте некоторые примеры данных вместо файлов, которых у нас нет, или сделайте несколько файлов доступными.

person Robert Hijmans    schedule 20.02.2018
comment
Привет, Роберт, спасибо за ответ!! Пожалуйста, смотрите мое выше обновление относительно имен файлов. - person Mike; 20.02.2018
comment
Я расширил свой ответ, используя эти имена файлов - person Robert Hijmans; 20.02.2018

Поскольку на самом деле в моем наборе данных есть еще две переменные «tasmax» и «tasmin», помимо «tas» и «pr», подход Роберта заключался бы в том, чтобы много писать для моего случая. Поэтому я попробовал другой способ, который в итоге сработал, хотя и не перечисляет файлы каждой переменной отдельно (недостаток, да!).

Список и сопоставление файлов исторических и rcp:

Чтобы сопоставить файлы, мне нужны чистые имена файлов без каталога, иначе which(!hist %in% rcp) всегда будет FALSE (как показано Робертом).

hist ‹- list.files("/исторический") rcp ‹- list.files("/rcp26")

no.match.h ‹- which(!hist %in% rcp) no.match.r ‹- which(!rcp %in% hist)

Поскольку мне нужно для nc_open имя файла, включая каталог, я должен создать соответствующий список файлов и вычесть несоответствующие файлы.

hist.files ‹- list.files("/data/scratch/lorchdav/cordex_eur/monmean/исторический", full.names = TRUE) rcp.files ‹- list.files("/data/scratch/lorchdav/cordex_ber_mean/rcp26 ", полные имена = ИСТИНА)

hist.files.cl ‹- hist.files[-no.match.h] hist.files.cl

rcp.files.cl ‹- rcp.files[-no.match.r] rcp.files.cl

person Mike    schedule 02.03.2018