Я очень-очень благодарен всем, кто хочет помочь мне с проблемой, с которой я действительно застрял. Но заранее: это сложная тема, и я изо всех сил стараюсь объяснить, что я собираюсь делать со своим кодом. Речь идет о климатических данных в файлах 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? Спасибо большое!
summary
,plot
)? Есть ли NA (-32767 или -9999)?sum(..., na.rm = TRUE)
поможет? - person Tung   schedule 20.02.2018