Как добавить средние растры в пределах цикла, который создает растры? R

Question

У меня есть несколько каталогов с 700 + бинарными закодированными растрами, которые я беру в среднем на выходные растры в каталоге. однако в настоящее время я создаю растры 1 на 1 в цикле for, а затем загружаю вновь созданные растры обратно в R, чтобы взять сумму, чтобы получить общее количество осадков в месяц.

Однако, поскольку мне не нужны отдельные растры, только средний растровый, у меня есть догадка, что я мог бы сделать это все в 1 петле, а не сохранить растры, а только выходного среднего растра, но я иду короткий в том, как программировать это в R.

setwd("~/Desktop/CMORPH/Levant-Clip/200001") 

dir.output <- '~/Desktop/CMORPH/Levant-Clip/200001' ### change as needed to give output location 
path <- list.files("~/Desktop/CMORPH/MonthlyCMORPH/200001",pattern="*.bz2", full.names=T, recursive=T) 

for (i in 1:length(path)) { 
    files = bzfile(path[i], "rb") 
    data <- readBin(files,what="double",endian = "little", n = 4948*1649, size=4) #Mode of the vector to be read 
    data[data == -999] <- NA #covert missing data from -999(CMORPH notation) to NAs 
    y<-matrix((data=data), ncol=1649, nrow=4948) 
    r <- raster(y) 
    e <- extent(-180, 180, -90, 83.6236) ### choose the extent based on the netcdf file info 
    tr <- t(r) #transpose 
    re <- setExtent(tr,extent(e)) ### set the extent to the raster 
    ry <- flip(re, direction = 'y') 
    projection(ry) <- "+proj=longlat +datum=WGS84 +ellps=WGS84" 
    C_Lev <- crop(ry, Levant) ### Clip to Levant 
    M_C_Lev<-mask(C_Lev, Levant) 
    writeRaster(M_C_Lev, paste(dir.output, basename(path[i]), sep = ''), format = 'GTiff', overwrite = T) ###the basename allows the file to be named the same as the original 
} 
# 
raspath <- list.files ('~/Desktop/CMORPH/Levant-Clip/200001',pattern="*.tif",  full.names=T, recursive=T) 
rasstk <- stack(raspath) 
sum200001<-sum(rasstk) 
writeRaster(avg200001, paste(dir.output, basename(path[i]), sep = ''), format = 'GTiff', overwrite = T) ###the basename allows the file to be named the same as the original 

в настоящее время этот код занимает около 75 минут, чтобы выполнить, и у меня есть около 120 больше каталогов, чтобы идти, и я ищу для более быстрых решений.

спасибо за все и любые комментарии и ввод. лучше, Evan

Answer 1

Конкретизируя мой предыдущий комментарий, вы можете попробовать:

setwd("~/Desktop/CMORPH/Levant-Clip/200001") 

dir.output <- '~/Desktop/CMORPH/Levant-Clip/200001' ### change as needed to give output location 
path <- list.files("~/Desktop/CMORPH/MonthlyCMORPH/200001",pattern="*.bz2", full.names=T, recursive=T) 
raster_list = list() 
for (i in 1:length(path)) { 
    files = bzfile(path[i], "rb") 
    data <- readBin(files,what="double",endian = "little", n = 4948*1649, size=4) #Mode of the vector to be read 
    data[data == -999] <- NA #covert missing data from -999(CMORPH notation) to NAs 
    y<-matrix((data=data), ncol=1649, nrow=4948) 
    r <- raster(y) 
    if (i == 1) { 
    e <- extent(-180, 180, -90, 83.6236) ### choose the extent based on the netcdf file info 

    } 
    tr <- t(r) #transpose 
    re <- setExtent(tr,extent(e)) ### set the extent to the raster 
    ry <- flip(re, direction = 'y') 
    projection(ry) <- "+proj=longlat +datum=WGS84 +ellps=WGS84" 
    C_Lev <- crop(ry, Levant) ### Clip to Levant 
    M_C_Lev<-mask(C_Lev, Levant) 
    raster_list[[i]] = M_C_Lev 
} 
# 

rasstk <- stack(raster_list, quick = TRUE) # OR rasstk <- brick(raster_list, quick = TRUE) 
avg200001<-mean(rasstk) 
writeRaster(avg200001, paste(dir.output, basename(path[i]), sep = ''), format = 'GTiff', overwrite = T) ###the basename allows the file to be named the same as the original 

Использование «быстрых» вариантов в stack определенно должны ускорить вещи, в частности, если у вас есть много растров.

Другая возможность - сначала вычислить среднее значение, а затем выполнить «пространственное прохождение». Например:

for (i in 1:length(path)) { 
    files = bzfile(path[i], "rb") 
    data <- readBin(files,what="double",endian = "little", n = 4948*1649, size=4) #Mode of the vector to be read 
    data[data == -999] <- NA #covert missing data from -999(CMORPH notation) to NAs 

    if (i == 1) { 
    totdata <- data 
    num_nonNA <- as.numeric(!is.na(data)) 
    } else { 
totdata = rowSums(cbind(totdata,data), na.rm = TRUE) 
# We have to count the number of "valid" entries so that the average is correct ! 
num_nonNA = rowSums(cbind(num_nonNA,as.numeric(!is.na(data))),na.rm = TRUE) 
    } 
} 

avg_data = totdata/num_nonNA # Compute the average 

# Now do the "spatial" processing 

y<-matrix(avg_data, ncol=1649, nrow=4948) 
r <- raster(y) 
e <- extent(-180, 180, -90, 83.6236) ### choose the extent based on the netcdf file info 
tr <- t(r) #transpose 
re <- setExtent(tr,extent(e)) ### set the extent to the raster 
ry <- flip(re, direction = 'y') 
projection(ry) <- "+proj=longlat +datum=WGS84 +ellps=WGS84" 
C_Lev <- crop(avg_data, Levant) ### Clip to Levant 
M_C_Lev<-mask(C_Lev, Levant) 
writeRaster(M_C_Lev, paste(dir.output, basename(path[i]), sep = ''), format = 'GTiff', overwrite = T) ###the basename allows the file to be named the same as the original 

Это может быть быстрее или медленнее, в зависимости от «сколько» вы кадрирование исходные данные.

НТН,

Lorenzo

Answer 2

Я добавляю еще один ответ, чтобы прояснить и упростить вещи немного, и в связи с комментариями в чате. Приведенный ниже код должен делать то, что вы просите: то есть циклически перебирать файлы, читать «данные», вычислять сумму по всем файлам и преобразовывать их в растр с указанными размерами.

Обратите внимание, что для целей тестирования здесь я заменял свой цикл на именах файлов с простыми 1 720 циклом, а файл чтением с созданием массивов одинаковой длиной, как у вас заполнены значениями от 1 до 4, и некоторых NA!

totdata <- array(dim = 4948*1649) # Define Dummy array 
for (i in 1:720) { 
    message("Working on file: ", i) 
    data <- array(rep(c(1,2,3,4),4948*1649/4), dim = 4948*1649) # Create a "fake" 4948*1649 array each time to simulate data reading 
    data[1:1000] <- -999 # Set some values to NA 
    data[data == -999] <- NA #convert missing data from -999 

    totdata <- rowSums(cbind(totdata, data), na.rm = T) # Let's sum the current array with the cumulative sum so far 
} 

# Now reshape to matrix and convertt to raster, etc. 
y <- matrix(totdata, ncol=1649, nrow=4948) 
r <- raster(y) 
e <- extent(-180, 180, -90, 83.6236) ### choose the extent based on the netcdf file info 
tr <- t(r) #transpose 
re <- setExtent(tr,e) ### set the extent to the raster 
ry <- flip(re, direction = 'y') 
projection(ry) <- "+proj=longlat +datum=WGS84 +ellps=WGS84" 

Это создает "правильный" растр:

> ry 
class  : RasterLayer 
dimensions : 1649, 4948, 8159252 (nrow, ncol, ncell) 
resolution : 0.07275667, 0.1052902 (x, y) 
extent  : -180, 180, -90, 83.6236 (xmin, xmax, ymin, ymax) 
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory 
names  : layer 
values  : 0, 2880 (min, max) 

contatining сумму различных массивов: Вы можете заметить, что максимальное значение равно 720 * 4 = 2880 (только предостережение: если у вас есть ячейки, которые всегда находятся в NA, вы получите 0 вместо NA)

На моем ноутбуке это работает примерно через 5 минут!

На практике:

1. , чтобы избежать проблем с памятью, я не читаю в памяти все данные. Каждый из ваших массивов имеет более или менее 64 МБ, поэтому я не могу загрузить их все , а затем сделать сумму (если у меня нет 50 ГБ ОЗУ, чтобы выбросить - и даже в этот случай был бы медленным). Вместо этого я использую ассоциативную ассоциацию суммирования, вычисляя «кумулятивную» сумму при каждом цикле . Таким образом, вы работаете только с двумя 8-миллиметровыми массивами на времени: тот, который вы читаете из файла «i», и тот, который содержит текущую сумму.
2. , чтобы избежать ненужных вычислений здесь. Я суммирую непосредственно 0-мерные массивы, которые я получаю от чтения двоичного кода. Вам не нужно , чтобы изменить к матричным массивов в цикле, потому что вы можете сделать это на финале «суммированный» массив, который затем можно преобразовать в матричной форме

Я надеюсь, что это будет работать для вас и что я не пропущу что-то очевидное!

Насколько я понимаю, если использование этого подхода все еще медленное, у вас возникают проблемы в другом месте (например, при чтении данных: на 720 файлах, 3 секунды, потраченные на чтение для каждого файла, составляют примерно 35 минут обработки).

НТН,

Lorenzo