Calculating area of overlap between 2 raster layers with different origins and extents in R

Question

I'm completely new to asking questions on stackoverflow and more or less a novice at R (and programming in general) so bear with me.

I have ASCII files of species ranges which only show presence only. After scouring the far reaches of the internet, I've managed to upload, convert to raster, mask along desired borders (in my case, the coastline of Australia), and plot them so I can visualize the ranges on an unprojected map.

Having accomplished the qualitative aspect of this, I need to get to the quantitative aspect; that is, I need to calculate the degree of sympatry between species. In order to do that, I need to first calculate the area of overlap, which is where I run into problems. Here's what I've managed to do so far:

> 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)

y is a raster of the overlap between d and b masked over d (when I try to mask over b I get error saying that the extents are different). How do I calculate the area of it? the area() function from the raster package spits this out:

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)

I'm not entirely sure what to do with this. Is this even a good/accurate/correct way of getting areas? Why are the extents different between y and b but the same between d and y? Furthermore, what are the units of the values from area(y)? I suppose the units don't matter too much because I will be taking a ratio eventually (by dividing the overlap by the range of the more restricted species), but I am curious to know for future reference.

I'm sorry if this is a stupid question. I appreciate any input someone might have.

Answer 1

The best way to get the overlap is with intersect. You can create a brick of the overlapping values and use a command like any to get the overlapping ranges, assuming the value is 1 or TRUE in each range and 0, FALSE, or NA outside the range:

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")

The area function gets you the approximate area of each cell (in order to get the true area, you should reproject it to an area coordinate system). To get the total approximate area of the overlap, add all the cell values, index the area summation by the values of the combined raster:

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