Запрос слоя растрового кирпича на основе другого растра в R

Question

У меня есть файл 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 одинаковы, как они кажутся.

Чтобы ускорить процесс, я обрезал глобальный кирпичик и растровый растр в своей области интересов, снова проверив, что все пространственные разрешения и т. Д. И т. Д. Совпадают - я не включил эти шаги для простоты.

В убыток. Любая помощь приветствуется!

Answer 1

Без воспроизводимого примера такие проблемы трудно решить. Я не могу сказать, где ваша проблема, но я представлю вам подход, который я бы попробовал. Может быть, это хорошо, может быть, плохо, я не знаю, но это может вдохновить вас найти способ обойти вашу проблему.

Насколько я понимаю, у вас есть кирпич OmegaA (33 слоя / глубина) и растр батиметрии. Вы хотите получить значение OmegaA на дне моря. Вот как бы я это сделал:

1. Сделать растр OmegaA в том же разрешении и степени, что и батиметрический
2. Преобразуйте растр батиметрии в растровый кирпич из 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)

3. затем вы умножаете оба своих растровых кирпича. Они имеют одинаковое измерение, это должно быть легко. На выходе должен быть растровый кирпич из 33 слоев с 0 везде, где это не дно моря и значение OmegaA где-либо еще.
4. Объедините весь полученный ранее слой кирпича в простой растр с суммой.

Это должно работать. Если у вас есть проблемы с обработкой растрового кирпича, вы можете превратить данные в базовые массивы R, это может быть проще.

Удачи.