Я бы возразил против вашего предположения, что это "очень большой" RasterStack, но, кроме того, я думаю, что то, что вы хотите сделать, должно быть прямо вперед.
Итак, сначала я загружаю и обрезаю данные в пределах Германии:
library(raster)
library(ncdf4)
library(R.utils)
library(maptools)
r <- stack('tg_0.25deg_reg_1995-2010_v17.0.nc')
data(wrld_simpl)
Germany <- wrld_simpl[wrld_simpl@data$NAME == "Germany",]
r_crop <- crop(r,Germany)
#Let's take a look:
plot(r_crop[[1]])
plot(Germany,add=T)
Форма границы не очень красивая, но она делает свою работу. Кроме того, вы можете видеть, что на севере значения с NoData по-прежнему правильно обозначаются так:
r_crop[[1]][1,1]
# NA
На следующих шагах я просто использую имена слоев, чтобы извлечь годы, а затем использую lapply, чтобы вычислить средние значения для каждого года:
nms <- names(r_crop)
yrs <- unique(sub('X(\\d+).+','\\1',nms))
yrs[1]
# [1] "1995"
annual_means <- lapply(yrs,function(x) mean(r_crop[[grep(x,nms)]],na.rm=TRUE))
Это даст вам список с именем annual_means с растром в год, представляющим среднегодовое значение за этот год. Теперь вы можете либо объединить их вместе (с do.call(stack,annual_means)), обработать их по отдельности, либо, как вы, вероятно, захотите, записать их на диск как csv:
# first take a look
plot(annual_means[[1]])
# write to disk
write.table(as.matrix(annual_means[[1]]),'ANNUAL_MEAN_TEMP_1995.csv',quote = F,row.names = F,col.names = F,sep = ';')
Редактировать
annual_means - список с растром на элемент, представляющий среднюю температуру, рассчитанную на основе ежедневных наблюдений исходного набора данных. Поэтому в списке будет столько элементов, сколько и лет.
Приведенный выше пример write.table был показан только в течение одного из этих лет. Это означает, что если вы хотите получить такой результат, вам нужно будет повторить шаг для всех элементов списка.
То, что делает шаг write.table, - это просто преобразование растра в матрицу и запись его на диск. Результатом будет матрица с таким количеством строк и столбцов, что и сам растр, с каждой ячейкой, разделенной точкой с запятой (мое личное предпочтение).
Edit2:
Просто чтобы проиллюстрировать мои комментарии сверху:
У вас есть данные за 16 лет, как видно из вектора yrs:
yrs
#[1] "1995" "1996" "1997" "1998" "1999" "2000" "2001" "2002" "2003" "2004"
#[11] "2005" "2006" "2007" "2008" "2009" "2010"
Теперь annual_means - это список длиной 16 с растром в один слой в год, который представляет собой среднее значение за весь год, рассчитанное для всей Германии из ежедневных данных.
Вот пример вывода:
annual_means[[1]]
# class : RasterLayer
# dimensions : 31, 37, 1147 (nrow, ncol, ncell)
# resolution : 0.25, 0.25 (x, y)
# extent : 5.75, 15, 47.25, 55 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# data source : in memory
# names : layer
# values : 3.329288, 11.32734 (min, max)
Как вы можете видеть, растр имеет разрешение 0,25 градуса (что является исходным разрешением ваших данных), что приводит к растру с 31 строкой и 37 столбцами, покрывающими Германию.
Чтобы получить желаемый результат:
Сначала я назову записи списка с соответствующими годами, чтобы сделать его немного более заметным (вы можете пропустить это):
names(annual_means) <- yrs
Теперь я извлеку координаты для каждого растра и создам кадр данных со значениями (используя lapply для перебора списка):
result <- lapply(annual_means, function(x) data.frame(long = coordinates(x)[,1],lat = coordinates(x)[,2],temp_mean =x[]))
Теперь мы можем проверить верхнюю часть кадра данных, например, 2000 год:
head(result$`2000`)
# long lat temp_mean
# 1 5.875 54.875 NaN
# 2 6.125 54.875 NaN
# 3 6.375 54.875 NaN
# 4 6.625 54.875 NaN
# 5 6.875 54.875 NaN
# 6 7.125 54.875 NaN
Как видите, все первые пиксели - это NoData (как показано на графике), что вам и нужно.
Таким образом, в итоге result - это список, каждый элемент которого представляет собой фрейм данных определенного года, содержащий столбцы long, lat и temp_mean.
Чтобы на 100% повторить желаемый результат, теперь можно повторить цикл по списку result, чтобы изменить имя столбца temp_mean на конкретный год (это совершенно необязательно):
for (i in seq_along(result)){
colnames(result[[i]])[3] <- paste0('Net_',names(result)[i],'_Temp')
}
Даю вам:
head(result$`2000`)
# long lat Net_2000_Temp
# 1 5.875 54.875 NaN
# 2 6.125 54.875 NaN
# 3 6.375 54.875 NaN
# 4 6.625 54.875 NaN
# 5 6.875 54.875 NaN
# 6 7.125 54.875 NaN
Edit3:
Чтобы получить один фрейм данных всеми средствами, вы можете сделать это:
ann_mean_temp_1996_1999 <- cbind(result[[1]][,1:2],do.call(cbind,lapply(result,function(x) x[,3])))
colnames(ann_mean_temp_1996_1999)[3:ncol(ann_mean_temp_1996_1999)]<- unlist(lapply(result,function(x) colnames(x)[3]))
Первый lapply связывает long / lat (который не изменяется за все годы) вместе с 3-м столбцом каждого элемента списка (который является T-MEAN).
Второй lapply извлекает и снова присваивает имена столбцов для температур, которые, похоже, теряются в процессе. Вероятно, есть более элегантное решение для этого, чем два раза использовать lapply, но оно делает свою работу.