Spatial polygon plots in wrong location on map background

Question

I want to plot the utilization distributions of two animals on a spatial map background and calculate area overlap. However, for some reason, although both layers have the same projection, one is off (plotted too far to the east). Consequently, any overlap I calculate then is incorrect.

So my question is what causes this issue and how can I fix this?

Below are my two rasters, transformed to SpatialPolygons so that I could use gIntersection() to calculate overlap. Note that while p1 plots just fine, p2 does not.

p1 <- class       : SpatialPolygons 
features    : 549 
extent      : 667950.6, 672950.6, 2840181, 2853981  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=17 +datum=WGS84 +units=m +no_defs 

p2 <- class       : SpatialPolygons 
features    : 257 
extent      : 670158.7, 673958.7, 2839623, 2851623  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=17 +datum=WGS84 +units=m +no_defs 

Below pic shows how the polygons are currently plotted. p1 is depicted in green, and p2 is in red.

Desired plotting location for p2 should look like this (this plot is created in ggplot2 from the same raster file):

I thought maybe it was because the extents are so different. So I extended the areas of both rasters to the same extent, after which I converted them to polygons and plotted. But this does not solve anything.

EDIT: rasters r1 and r2 that can be converted to SpatialPolygons using as(x, 'SpatialPolygons') can be downloaded here and here.

EDIT 2: The crs for both rasters should be:

"+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs"

EDIT 3: If it helps, the two rasters in question were generated by combining 2 sub-region rasters (central and north) after first extending the spatial extents of each raster to the spatial extent shared by both rasters and then rescaling to range from 0 to 1. However, some other rasters (not included here) are made up from 3 sub-region rasters. I don't know why this could matter, as the projections of the rasters are the same, so therefore should lead to correct spatial plots?

EDIT 4: Region rasters that were used to create r1 can be downloaded here and here. Similarly, region rasters for r2 are here and here.

EDIT 5: Below the code that i used to generate the output. Sorry for not introducing this earlier. r1 and r2 provided in the links are identical to r1.recl and r2.recl where the 0's have been replaced by NA.

     library(magrittr)
aeqd.crs <- "+proj=aeqd +lat_0=25.6871 +lon_0=-79.29336 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" #aeqd is original CRS
    UTMstring <- "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs" # crs of study site shapefile
        r1 = raster("r1")
        proj4string(r1) <- CRS(aeqd.crs)
        r1 %<>% raster::setMinMax()   
        r1 <- projectRaster(r1, crs = UTMstring) 
        r1[is.na(r1)] <- 0
        r1 <- r1 / max(r1@data@values)
        recl <- matrix(data = c(0, 0.05, 0,
                                0.05, 1, 1), ncol = 3, byrow = T)
        r1.recl <- reclassify(r1, recl)
        r1.recl[r1.recl == 0] <- NA
        r2 = raster("r2")
        proj4string(r2) <- CRS(aeqd.crs) 
        r2 %<>% raster::setMinMax()
        r2 <- projectRaster(r2, crs = UTMstring)
        r2[is.na(r2)] <- 0
        r2 <- r2 / max(r2@data@values)
        r2.recl <- reclassify(r2, recl)
        r2.recl[r2.recl == 0] <- NA
        
        # convert rasters to spatial polygons; vector shapes give access to spatial tools in the {rgeos} package
        p1 <- as(r1.recl, 'SpatialPolygons')
        p2 <- as(r2.recl, 'SpatialPolygons')
        
    overlap <- gIntersection(p1, p2) # Creates raster
    
        # Calculate common area
        plot(p1, col = "green")
        plot(p2, col = "red", add = TRUE)
    plot(overlap, col = "orange", add = TRUE)
        
        # calculate area of overlap (in m2)
        common.area <- gArea(overlap) / (1000^2)
        common.area
[1] 92716.54

Answer 1

Here are two approaches to compute overlap with your data.

Your data

library(terra)
prj <- "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs"
r1 <- rast("r1.asc")
crs(r1) <- prj
r2 <- rast("r2.asc")
crs(r2) <- prj

The "raster" approach

e <- union(ext(r1), ext(r2))
x1 <- extend(r1, e)
x2 <- resample(r2, x1) + 1
y <- sum(c(x1, x2), na.rm=TRUE)    
a <- cellSize(y)
zonal(a, y, "sum")
#  sum        area
#1   1 21881870.89
#2   2 10200711.13
#3   3    80006.73

Where "1" is the area covered by species r1, "2" by r2 and "3" by both.

You can also assume a fixed nominal resolution. In this case that is only a little less precise (but that is not always the case)

f <- freq(y)
f$area <- f$count * prod(res(y))
f[,-1]
#  value count     area
#1     1   547 21880000
#2     2   255 10200000
#3     3     2    80000

The "vector" approach

p1 <- as.polygons(r1)
p2 <- as.polygons(r2)
p <- union(p1, p2)
p$area <- expanse(p)
data.frame(p)
#  r1 r2        area
#1  1  1    92724.32
#2  1 NA 21869156.44
#3 NA  1 10187993.21
 

You could also do

intersect(p1, p2) |> expanse()
#[1] 92724.32

The results are similar for the raster and vector approach, but not the same. This is because the raster data are not aligned, which suggests that you could improve how these were created.

Plotting the two species shows the same pattern you have

p1$sp <- "1"
p2$sp <- "2"
p <- rbind(p1, p2)
plot(p, "sp")

And I can more-or-less reproduce r1 and r2 (but I need to guess things, in the future please describe your procedures with code, not with natural language) :

bn = rast("bull_North.asc")
bs = rast("bull_South.asc")
rn = rast("reef_North.asc")
rs = rast("reef_South.asc")
b2 = bn + bs
r2 = rn + rs
e <- intersect(ext(b2), ext(r2))
b <- crop(b2, e)
r <- crop(r2, e)
b <- ifel(b>0.05, 1, NA)
r <- ifel(r>0.05, 1, NA)
plot(b, col="blue")
plot(r, col="red", add=TRUE)

So, clearly the problem is already present in the original data.