2

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.

  • @JeffreyEvans I should have mentioned that in ArcGIS I also used the aggregate tool, so the resampling is not the issue. I need to extract values from raster pftnc in those cells where the raster sp has value=1. Maybe I should try to avoid the on-the-fly projection and reproject the raster myself. – Herman Toothrot Oct 23 '18 at 21:01
  • Keep in mind that reprojection requires resampling as well, and the global raster resampling method defaults to nearest. So, unless you run the tool the results will be incorrect. These data preparation steps should be done explicitly so you do not end up with unexpected results. – Jeffrey Evans Oct 23 '18 at 21:07
  • Try either projectRaster(pftnc, sp) or resample(pftnc, sp) - you can reverse the order if you prefer – mdsumner Oct 23 '18 at 22:34
2

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.

  • I appreciate the advice and I tried your suggestions. The aggregate function in R and ArcGIS produces the same results, so that's good. However, after I use "extract values to points" in ArcGIS or stack+rasterToPoints in R, I get different results. In addition, if I also reproject the raster in ArcGIS and then extrat values, I get a 3rd result. So I think there are some minor inconsistencies between the R reprojection and the ArcGIS reprojection (onthefly or via reprojection tool). I am not sure which result to 'believe'. – Herman Toothrot Oct 24 '18 at 10:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.