Вычисление области перекрытия между двумя растровыми слоями различного происхождения и экстента в R

Question

Я совершенно не задаю вопросы о stackoverflow и более или менее новичком в R (и программировании в целом), так что несите меня.

У меня есть ASCII-файлы видов, которые показывают только присутствие. После чистки дальних интернета мне удалось загрузить, конвертировать в растре, маску вдоль желаемых границ (в моем случае, береговую линию Австралии) и нарисовать их, чтобы я мог визуализировать диапазоны на непроектированной карте.

Достигнув качественного аспекта этого, мне нужно перейти к количественному аспекту; то есть мне нужно рассчитать степень симпатии между видами. Для этого мне нужно сначала вычислить область перекрытия, где я столкнулся с проблемами. Вот что мне удалось сделать до сих пор:

> d 
class  : RasterLayer 
dimensions : 85, 270, 22950 (nrow, ncol, ncell) 
resolution : 0.08, 0.08 (x, y) 
extent  : 119.4993, 141.0993, -36.65831, -29.85831 (xmin, xmax, ymin, ymax) 
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory 
names  : layer 
values  : 2, 2 (min, max) 

> b 
class  : RasterLayer 
dimensions : 140, 222, 31080 (nrow, ncol, ncell) 
resolution : 0.08, 0.08 (x, y) 
extent  : 134.2456, 152.0056, -40.44268, -29.24268 (xmin, xmax, ymin, ymax) 
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory 
names  : layer 
values  : 2, 2 (min, max) 

x<-resample(b,d,method="ngb") 
y<-mask(x,d) 

>y 
class  : RasterLayer 
dimensions : 85, 270, 22950 (nrow, ncol, ncell) 
resolution : 0.08, 0.08 (x, y) 
extent  : 119.4993, 141.0993, -36.65831, -29.85831 (xmin, xmax, ymin, ymax) 
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory 
names  : layer 
values  : 2, 2 (min, max) 

у растровый перекрытия между d и Ь маскируется над д (при попытке замаскировать над б я получаю ошибку о том, что экстентов разные) , Как рассчитать его площадь? функция области() из пакета растрового плюет это:

area(y) 
class  : RasterLayer 
dimensions : 85, 270, 22950 (nrow, ncol, ncell) 
resolution : 0.08, 0.08 (x, y) 
extent  : 119.4993, 141.0993, -36.65831, -29.85831 (xmin, xmax, ymin, ymax) 
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory 
names  : layer 
values  : 63.65553, 68.75387 (min, max) 

Я не совсем уверен, что делать с этим. Это даже хороший/правильный/правильный способ получения областей? Почему экстенты различаются между y и b, но то же самое между d и y? Кроме того, каковы единицы значений из области (y)? Я полагаю, что единицы не имеют большого значения, потому что я буду принимать соотношение в конечном счете (разделив перекрытие на диапазон более ограниченных видов), но мне любопытно узнать для дальнейшего использования.

Прошу прощения, если это глупый вопрос. Я ценю любой вклад, который может быть у кого-то.

Answer 1

Лучший способ получить перекрытие - intersect. Вы можете создать кирпич перекрывающихся значений и использовать команду типа any, чтобы получить перекрывающиеся диапазоны, предполагая, что значение равно 1 или TRUE, в каждом диапазоне и 0, FALSE, или NA вне диапазона:

library(raster) 

wgs84 <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" 
d <- raster() 
extent(d) <- extent(119.4993, 141.0993, -36.65831, -29.85831) 
res(d) <- c(0.08, 0.08) 
projection(d) <- CRS(wgs84) 
values(d) <- sample(c(NA, 1), ncell(d), replace=TRUE) 

b <- raster() 
extent(b) <- c(134.2456, 152.0056, -40.44268, -29.24268) 
res(b) <- c(0.08, 0.08) 
projection(b) <- CRS(wgs84) 
values(b) <- sample(c(NA, 1), ncell(b), replace=TRUE) 

y <- intersect(b, d) 

x <- brick(resample(b, y, method = "ngb"),resample(d, y, method = "ngb")) 
x2 <- any(x, na.rm = TRUE) 

library(maps) 
map(regions = "australia") 
image(d, add = TRUE, col = "blue") 
image(b, add = TRUE, col = "green") 
plot(extent(y), add = TRUE) 
image(x2, add = TRUE, col = "red") 

area функция получает вас область каждый клеток приблизительная (для того, чтобы получить истинную площадь, вы должны перепроецировать его в области системы координат). Для того, чтобы получить общую приблизительную площадь перекрытия, добавить все значения ячеек, индекс область суммирования по значениям объединенного растра:

sum(values(area(x2))[which(values(x2))]) 
# [1] 361407.1