Extract values from raster, R and ArcGIS different results

Question

My goal is to extract values from the raster called pftnc by using the other raster sp after converting it to a point shapefile. Files available here pftnc and sp

I have two raster files

pftnc=raster("pftnc.tif")
sp=raster("sp.tif")
sp
    class       : RasterLayer 
    dimensions  : 142, 360, 51120  (nrow, ncol, ncell)
    resolution  : 96486, 96515  (x, y)
    extent      : -17367529, 17367529, -6356742, 7348382  (xmin, xmax, ymin, ymax)
    coord. ref. : +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names       : Puma_concolor 
values      : 0, 1  (min, max)

pftnc
    class       : RasterLayer 
    dimensions  : 720, 1440, 1036800  (nrow, ncol, ncell)
    resolution  : 0.25, 0.25  (x, y)
    extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
    coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
    value     : 0,1

In R, first I use aggregate to upscale one to have the same dimensions, I crop and then I project.

pftagg <- aggregate(pftnc, fact=4)
pftext <- c(xmin= -180, xmax= 180, ymin= -60, ymax= 90)
pftagg=crop(pftagg,pftext)
pft1 <- projectRaster(pftagg, sp)

Then I convert raster to points and extract the values

sppoints=rasterToPoints(sp,fun=function(x){x==1})
result=extract(pft1,sppoints[,1:2])

I do the equivalent operations in ArcGIS 10.5, but I don't reproject or crop because when I import the layers everything lines up perfectly.

So the problem now is that I have inconsistencies between the two methods. Aside from a few missing points from the ArcGIS operation, I get slightly different means and standard deviation. I also tried to import in ArcGIS the raster after cropping, aggregating, etc, and I still get these differences. So things in ArcGIS are at least consistent with itself. I am wondering if these small differences are due to how the two software are dealing with projecting the latlong raster to CEA projection.

Answer 1

When you say "everything lines up perfectly" please be mindful that ArcGIS is doing quite a bit behind the scenes that you may not be aware of. Besides on-the-fly projection, if the data are of a difference resolution a resampling is performed, which defaults to nearest. If your data is continuous than this default is quite incorrect.

I would highly recommend setting an analysis environment in ArcGIS and then explicitly resampling your data using binlinear to create a new raster that aligns, or do the equivalent in R. These are the rasters that you would use in both analysis.

I would also ask, why are you creating points and extracting values? If you are matching resolution and extent, than all you have to do is stack the rasters and then use raster::rasterToPoints. This will write out an sp SpatialPointsDataFrame with values for all the rasters in the stack. You can also coerce a raster or stack into a SpatialPixelsDataFrame. Your entire raster::extract step is unnecessary. Besides, it is really not good practice to coerce rasters into vectors. Raster data structures are very optimal due to storing arrays and not needing [x,y] coordinates for each pixel. The raster::overlay, raster::calc and raster::focal functions can cover much of your analytical needs as well as remaining memory safe.

Here is an illustration where we create a binary and a gaussian raster. We then apply a function that turns all data in the gaussian raster to NA if the corresponding cell in the binary raster is not equal to 1.

First, add libraries and create some data

library(raster)
library(sp)

r <- raster(nrows=180, ncols=360, xmn=571823.6, xmx=616763.6, ymn=4423540, 
            ymx=4453690, resolution=270, crs = CRS("+proj=utm +zone=12 +datum=NAD83 
            +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
    r[] <- rpois(ncell(r), lambda=1)
    r <- calc(r, fun=function(x) { x[x >= 1] <- 1; return(x) } )

r2 <- r
r2[] <- runif(ncell(r))

r <- stack(r,r2)
  names(r) <- c("binary","gauss")
    plot(r)

Now, apply a function to the stack to create NA if the binary raster is not equal to 1.

r[[2]] <- calc(r, fun=function(x) { ifelse( x[1] != 1, NA, x[2]) } )    
  plot(r)

Now, you can write out the raster (r[[2]]) or coerce to a SpatialPointsDataFrame and write to a shapefile.