Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Cryptic polygon issue leading to incorrect rasterization

Tags:

r

raster

r-sf

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')

enter image description here 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)

enter image description here

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.

like image 206
Pascal Avatar asked Feb 18 '26 02:02

Pascal


1 Answers

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)

enter image description here

like image 96
sebdalgarno Avatar answered Feb 20 '26 16:02

sebdalgarno



Donate For Us

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