-
Notifications
You must be signed in to change notification settings - Fork 15
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Outputs are not identical to raster::rasterize() #6
Comments
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. |
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 |
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!) |
Thanks! Geospatial stuff is fairly new to me, so glad to have the input. |
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,
|
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
|
Due to some rounding inconsistency, outputs from
fasterize()
are not exactly the same asrasterize()
, as in the example below. This could be in the implementation ofstd::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.Session info
The text was updated successfully, but these errors were encountered: