I have a number of spatial polygons that I am rasterizing, using the sf and raster packages. It has worked out great, except for this one case, where I have run out of ideas as to how to fix the problem:
I have a polygon as a simple feature polygon, and a raster template. Both have been transformed to an equal area Behrmann projection. There was a geometry validity issue with the polygon, which I fixed with st_make_valid from the lwgeom package.
library(sf)
library(lwgeom)
library(raster)
spfile <- 'sp.rds'
rasterfile <- 'rasterTemplate_150km.tif'
sp <- readRDS(spfile)
rasterTemplate <- raster(rasterfile)
# are there any geometry validity issues? yes!
st_is_valid(sp)
[1] FALSE
Warning message:
In evalq((function (..., call. = TRUE, immediate. = FALSE, noBreaks. = FALSE, :
Ring Self-intersection at or near point 9947244.3466675151 1450099.5373749144
sp2 <- st_make_valid(sp)
st_is_valid(sp2)
[1] TRUE
plot(sp2, col='blue')
Visually, the polygon looks good.
Now I rasterize it (I need to convert to SpatialPolygon for compatibility with raster):
sp3 <- as(sp2, 'Spatial')
cover <- rasterize(sp3, rasterTemplate, getCover=TRUE)
plot(cover, box=F, axes=F)

The problem is that there is a bar going through Australia.
Does this mean that there is still a geometry validity issue with that polygon that is not detected by st_is_valid? If so, that makes this type of issue hard to work with, as it is not easily detected, except visually...
Going back, I can confirm that the rasterization problem happens with as(sp, 'Spatial), so st_make_valid doesn't create this problem.
How can I fix this problem?
UPDATES
A suggestion below was made to use st_cast. I tried sp4 <- as(st_cast(sp, "POLYGON"), 'Spatial'). The same rasterization problem remains.
I also tried buffering by 0 with sp5 <- as(st_buffer(sp, dist = 0), 'Spatial'), but this also does not solve the problem.
The polygon can be downloaded as a .rds file here. The raster can be downloaded here.
The fasterize function from fasterize package apparently does not cause the same error. Fasterize will also be incorporated into the stars package (see this issue), which should be on CRAN soon.
There is also no need to st_cast or st_make_valid.
library(sf)
library(fasterize)
df <- readRDS('sp.rds') %>% st_sf(field = 1)
template <- raster('rasterTemplate_150km.tif')
r <- fasterize(df, template, field = 'field')
plot(r)

If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With