Outputs are not identical to raster::rasterize()

[noamross]

Due to some rounding inconsistency, outputs from fasterize() are not exactly the same as rasterize(), as in the example below. This could be in the implementation of std::ceil as opposed to whatever rounding or edge rule is used in raster. I'm not sure whether this is a problem or not, but it would be good to know why it is.

library(fasterize)
#> 
#> Attaching package: 'fasterize'
#> The following object is masked from 'package:graphics':
#> 
#>     plot
library(raster)
#> Loading required package: sp
library(sf)
#> Linking to GEOS 3.4.2, GDAL 2.1.2, proj.4 4.9.1
p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p1 <- list(p1, hole)
p2 <- list(rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0)))
p3 <- list(rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0)))
pols <- st_sf(value = rep(1,3),
             geometry = st_sfc(lapply(list(p1, p2, p3), st_polygon)))
pols_sp <- as(pols, "Spatial")
r <- raster(pols, res = 1)
f <- fasterize(pols, r, field = "value", fun="sum")
r <- rasterize(pols_sp, r, field = "value", fun="sum")
identical(r, f)
#> [1] FALSE
plot(r-f)

Session info

devtools::session_info()
#> Session info -------------------------------------------------------------
#>  setting  value                       
#>  version  R version 3.3.2 (2016-10-31)
#>  system   x86_64, darwin13.4.0        
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  tz       America/New_York            
#>  date     2017-03-19
#> Packages -----------------------------------------------------------------
#>  package   * version     date       source                            
#>  backports   1.0.5       2017-01-18 cran (@1.0.5)                     
#>  bitops      1.0-6       2013-08-17 CRAN (R 3.3.0)                    
#>  DBI         0.6         2017-03-10 Github (rstats-db/DBI@4d9571a)    
#>  devtools    1.12.0.9000 2017-02-13 Github (hadley/devtools@d8ab190)  
#>  digest      0.6.12      2017-01-27 cran (@0.6.12)                    
#>  evaluate    0.10        2016-10-11 cran (@0.10)                      
#>  fasterize * 0.1.0       2017-03-20 local                             
#>  htmltools   0.3.5       2016-03-21 CRAN (R 3.3.0)                    
#>  knitr       1.15.1      2016-11-22 cran (@1.15.1)                    
#>  lattice     0.20-34     2016-09-06 CRAN (R 3.3.2)                    
#>  magrittr    1.5         2014-11-22 CRAN (R 3.3.0)                    
#>  memoise     1.0.0.9001  2016-12-16 Github (hadley/memoise@884d565)   
#>  pkgbuild    0.0.0.9000  2016-11-28 Github (r-pkgs/pkgbuild@65eace0)  
#>  pkgload     0.0.0.9000  2017-02-13 Github (r-pkgs/pkgload@fc907a1)   
#>  raster    * 2.5-8       2016-06-02 CRAN (R 3.3.0)                    
#>  Rcpp        0.12.9.4    2017-03-09 Github (RcppCore/Rcpp@0566d7c)    
#>  RCurl       1.95-4.8    2016-03-01 CRAN (R 3.3.0)                    
#>  rmarkdown   1.3.9004    2017-03-16 Github (rstudio/rmarkdown@d5d22b9)
#>  rprojroot   1.2         2017-01-16 cran (@1.2)                       
#>  sf        * 0.3-4       2017-02-06 CRAN (R 3.3.2)                    
#>  sp        * 1.2-4       2016-12-22 cran (@1.2-4)                     
#>  stringi     1.1.2       2016-10-01 cran (@1.1.2)                     
#>  stringr     1.2.0       2017-02-18 cran (@1.2.0)                     
#>  udunits2    0.13        2016-11-17 cran (@0.13)                      
#>  units       0.4-2       2017-01-13 cran (@0.4-2)                     
#>  withr       1.0.2       2016-06-20 CRAN (R 3.3.0)                    
#>  XML         3.98-1.5    2016-11-10 cran (@3.98-1.)                   
#>  yaml        2.1.14      2016-11-12 cran (@2.1.14)

[mdsumner]

I've seen a difference in point in cell tests we used tabulate for in much older code, it's something I wanted to figure out. It would be nice to find it was the same reason so I will try to have a look.

[mdsumner]

Also pretty sure raster's tests leave the right (and top?) side as open, i.e. points on the boundary are counted as outside. It might be worth considering

[mdsumner]

I don't have the answer but you can see differences at different resolutions, here using intersection to identify which tests are different and plotting the value given by fasterize. At resolution 1 or 10 fasterize over counts rasterize, but at even values between those it under counts.

library(fasterize)
library(raster)
library(sf)
p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p1 <- list(p1, hole)
p2 <- list(rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0)))
p3 <- list(rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0)))
pols <- st_sf(value = rep(1,3),
              geometry = st_sfc(lapply(list(p1, p2, p3), st_polygon)))
pols_sp <- as(pols, "Spatial")
## at res = 1 or 10 fasterize gets value 2, when rasterize gets value 1
## but at 2, 4, 6, and 8 fasterize under counts what rasterize gets by 1
r <- raster(pols, res = 10)  ## values of 2, 4, 6, 8
centrepoints <- st_multipoint(coordinates(r))
f <- fasterize(pols, r, field = "value", fun="sum")
r <- rasterize(pols_sp, r, field = "value", fun="sum")
plot(abs(r-f) > 0); plot(st_geometry(pols), add = TRUE)
plot(st_intersection(pols, centrepoints), add = TRUE, cex = 1:3)

asub <- abs(values(r) - values(f)) > 0
text(coordinates(r)[asub, ], lab = values(f)[asub])

Probably not that useful, but just in case it is.

(BTW, I really appreciate this package as example code, it's really great for learning how to go about this level of programming!)

[noamross]

Thanks! Geospatial stuff is fairly new to me, so glad to have the input.

[mikoontz]

A puzzle! I'd love to use this package for some work on gigantor multipolygons (CalVeg designations of California Wildlife Habitat Relationship types-- essentially type of plant cover). I'll dig around a bit, too!

Quick note that, even when all the raster values are the same, identical(r, f) will return FALSE because the @DaTa@names slots are different:

library(viridis)
library(fasterize)
library(raster)
library(sf)

p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p1 <- list(p1, hole)
p2 <- list(rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0)))
p3 <- list(rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0)))
pols <- st_sf(value = rep(1,3),
              geometry = st_sfc(lapply(list(p1, p2, p3), st_polygon)))
pols_sp <- as(pols, "Spatial")
r <- raster(pols, nrow = 1, ncol = 1)
f <- fasterize(pols, r, field = "value", fun="sum")
r <- rasterize(pols_sp, r, field = "value", fun="sum")
identical(r, f)
#> [1] FALSE
plot(r-f, col = viridis(2))

r@data@names
[1] ""
f@data@names
[1] "layer"

f@data@names <- ""

identical(r, f)
#> [1] TRUE

[mikoontz]

Here's some code I've been using to flag non-identical pixel values across raster layers. Just subtracting the rasters can be misleading because subtraction on an NA will always return NA, which then won't be plotted. So this test is whether the values are different or whether one of the two pixels is an NA.

library(fasterize)
library(raster)
library(sf)

# Function that determines, pixel-by-pixel, whether two raster layers differ
rasterCheck <- function(x, y) {
  setValues(x, 
            ifelse(
              test = (x[] != y[]) | (is.na(x[]) + is.na(y[]) == 1), 
              yes = 1, 
              no = NA)
            )
}

p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p1 <- list(p1, hole)
p2 <- list(rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0)))
p3 <- list(rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0)))
pols <- st_sf(value = rep(1,3),
              geometry = st_sfc(lapply(list(p1, p2, p3), st_polygon)))
pols_sp <- as(pols, "Spatial")
r <- raster(pols, nrow = 22, ncol = 22)
f <- fasterize(pols, r, field = "value", fun="sum")
r <- rasterize(pols_sp, r, field = "value", fun="sum")

plot(rasterCheck(r, f), col = "red", legend = FALSE, axes = FALSE, bty = "n")
plot(pols_sp, add = TRUE)
plot(as(r, "SpatialGrid"), add = TRUE)