Запрос слоя растрового кирпича на основе другого растра в R
У меня есть файл NetCDF глобальных океанографических данных (OmegaA) с относительно грубым пространственным разрешением с 33 уровнями глубины. У меня также есть глобальный растр батиметрии с гораздо более высоким разрешением. Моя цель - использовать данные OmegaA на морском дне из файла NetCDF, используя данные батиметрии для определения желаемой глубины. Мой код до сих пор;
library(raster)
library(rgdal)
library(ncdf4)
# Aragonite data. Defaults to CRS WGS84
ncin <- nc_open("C:/..../GLODAPv2.2016b.OmegaA.nc")
ncin.depth <- ncvar_get(ncin, "Depth")# 33 depth levels
omegaA.brk <- brick("C:/.../GLODAPv2.2016b.OmegaA.nc")
omegaA.brk <-rotate(omegaA.bkr)# because netCDF is in Lon 0-360.
# depth raster. CRS WGS84
r<-raster("C:/....GEBCO.tif")
# resample the raster brick to the resolution that matches the bathymetry raster
omegaA.brk <-resample(omegaA.brk, r, method="bilinear")
# create blank final raster
omegaA.rast <- raster(ncol = r@ncols, nrow = r@nrows)
extent(omegaA.rast) <- extent(r)
omegaA.rast[] <- NA_real_
# create vector of indices of desired depth values
depth.values<-getValues(r)
depth.values.index<-which(!is.na(depth.values))
# loop to find appropriate raster brick layer, and extract the value at the desired index, and insert into blank raster
for (p in depth.values.index) {
dep.index <-which(abs(ncin.depth+depth.values[p]) == min(abs(ncin.depth+depth.values[p]))) ## this sometimes results in multiple levels being selected
brk.level <-omegaA.brk[[dep.index]] # can be more than on level if multiple layers selected above.
omegaA.rast[p] <-omegaA.brk[[1]][p] ## here I choose the first level if multiple levels have been selected above
print(paste(p, "of", length(depth.values.index))) # counter to look at progress.
}
Проблема: в результате получается растр с массивными пробелами (NA), где должны быть данные. Разрывы часто принимают отличительную форму - например, следуют по контуру или по длинной прямой линии. Я вставил обрезанный пример.
введите описание изображения здесь
Я думаю, что это может быть связано с тем, что либо 1) по какой-то причине утверждение 'which' в цикле не находит совпадения, либо 2) создается несоосность проекций, которые, как я прочитал, могут возникнуть при использовании 'Rotate'.
Я пытался убедиться, что все экстенты, разрешения, количество ячеек и CRS одинаковы, как они кажутся.
Чтобы ускорить процесс, я обрезал глобальный кирпичик и растровый растр в своей области интересов, снова проверив, что все пространственные разрешения и т. Д. И т. Д. Совпадают - я не включил эти шаги для простоты.
В убыток. Любая помощь приветствуется!
1 ответ
Без воспроизводимого примера такие проблемы трудно решить. Я не могу сказать, где ваша проблема, но я представлю вам подход, который я бы попробовал. Может быть, это хорошо, может быть, плохо, я не знаю, но это может вдохновить вас найти способ обойти вашу проблему.
Насколько я понимаю, у вас есть кирпич OmegaA (33 слоя / глубина) и растр батиметрии. Вы хотите получить значение OmegaA на дне моря. Вот как бы я это сделал:
- Сделать растр OmegaA в том же разрешении и степени, что и батиметрический
- Преобразуйте растр батиметрии в растровый кирпич из 33 трех слоев 0-1. Например, если морское дно находится на 200 м для одного конкретного пикселя, то этот пиксель на всем глубинном слое, кроме 200, равен 0 и 1 для 200. Чтобы запрограммировать это, я бы пошел долгий путь,
:
r_1 <- r
values(r_1) <- values(r)==10 # where 10 is the depth (it could be a range with < or >)
r_2 <- r
values(r_2) <- values(r)==20
...
r_33 <- r
values(r_33) <- values(r)==250
r_brick <- brick(r_1, r_2, ..., r_33)
- затем вы умножаете оба своих растровых кирпича. Они имеют одинаковое измерение, это должно быть легко. На выходе должен быть растровый кирпич из 33 слоев с 0 везде, где это не дно моря и значение OmegaA где-либо еще.
- Объедините весь полученный ранее слой кирпича в простой растр с суммой.
Это должно работать. Если у вас есть проблемы с обработкой растрового кирпича, вы можете превратить данные в базовые массивы R, это может быть проще.
Удачи.