terra::vect does not work for larger polygons

I am trying to convert watershed boundaries delineated using the R package streamstats to SpatVector objects using terra::vect. I uploaded an Rdata file with a list of a subset of 6 of the sf object watersheds I am working with to google drive to make a reproducible example, and below is my code that shows that the last watershed in the list (which is also much larger than the other 5) does not work in vect:

library(sf)
library(terra) # terra 1.6.17

## import the data:

load('C:/PhD/Research/Code/Manuscript/top6.Rdata')

## Here is an example of a watershed that does work with terra::vect:

# convert the 5th watershed to a sf polygon object:

v <- st_cast(st_boundary(l_DA_6[[5]]$geometry), "POLYGON")

# plot it: 

plot(v) # looks good

# now convert the sf object to SpatVector:

SpatV<-vect(v)

# and plot the SpatVector:

plot(SpatV) # also looks good

## And here is the watershed that does not work with terra::vect:

# convert the 6th watershed to a sf polygon object:

v <- st_cast(st_boundary(l_DA_6[[6]]$geometry), "POLYGON")

# plot it: 

plot(v) # looks good

# now convert the sf object to SpatVector:

SpatV<-vect(v)

# and plot the SpatVector:

plot(SpatV) # just showing a rectangle which is wrong

The only thing I can tell that is different between the two sf objects is that the second one is much larger. I was having this issue with the CRAN version of terra and the developer version terra 1.6.17. Any help is much appreciated.

Solution

All sf objects seem to be invalid, once passed through st_make_valid() and after unifying classes, conversion to SpatVector looks fine.

Made an edit as resulting SpatVectors ended up being lines, not polygons (1st rev). Which led to some additional manipulations and a few small polygons were dropped, so you should probably validate if this is indeed an acceptable result.

Part of the underlying issue here is probably related to unprojected CRS and the use of s2 (note the sf_use_s2() is TRUE when loading sf). So perhaps consider changing to projected CRS during some earlier step.

library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(terra)
#> terra 1.7.39

# check validity
lapply(l_DA_6, st_is_valid, reason = TRUE) |> str()
#> List of 6
#>  $ BLACK CREEK AT CHURCHVILLE NY    : chr "Loop 0: Edge 0 is degenerate (duplicate vertex)"
#>  $ OATKA CREEK AT GARBUTT NY        : chr "Loop 0: Edge 78 is degenerate (duplicate vertex)"
#>  $ ALLEN CREEK NEAR ROCHESTER NY    : chr "Loop 0: Edge 11 is degenerate (duplicate vertex)"
#>  $ NORTHRUP CREEK AT NORTH GREECE NY: chr "Loop 0: Edge 59 is degenerate (duplicate vertex)"
#>  $ EAST BROOK EAST OF WALTON NY     : chr "Loop 0: Edge 17 is degenerate (duplicate vertex)"
#>  $ BLACK RIVER AT WATERTOWN NY      : chr "Loop 1: Edge 7 is degenerate (duplicate vertex)"


# make valid, check resulting classes
l2 <- lapply(l_DA_6, st_make_valid)
lapply(l2, st_geometry) |> 
  lapply(class) |>
  str()
#> List of 6
#>  $ BLACK CREEK AT CHURCHVILLE NY    : chr [1:2] "sfc_MULTIPOLYGON" "sfc"
#>  $ OATKA CREEK AT GARBUTT NY        : chr [1:2] "sfc_MULTIPOLYGON" "sfc"
#>  $ ALLEN CREEK NEAR ROCHESTER NY    : chr [1:2] "sfc_MULTIPOLYGON" "sfc"
#>  $ NORTHRUP CREEK AT NORTH GREECE NY: chr [1:2] "sfc_POLYGON" "sfc"
#>  $ EAST BROOK EAST OF WALTON NY     : chr [1:2] "sfc_POLYGON" "sfc"
#>  $ BLACK RIVER AT WATERTOWN NY      : chr [1:2] "sfc_MULTIPOLYGON" "sfc"

# some shapes ended up as multipolygons, 
# cast all to polygons and drop small artefacts (100m^2 area threshold)
l3 <- lapply(l2, st_geometry) |>
  lapply(st_cast, "POLYGON") |> 
  lapply(\(g) g[units::drop_units(st_area(g)) > 100])
lapply(l3, class) |> str()
#> List of 6
#>  $ BLACK CREEK AT CHURCHVILLE NY    : chr [1:2] "sfc_POLYGON" "sfc"
#>  $ OATKA CREEK AT GARBUTT NY        : chr [1:2] "sfc_POLYGON" "sfc"
#>  $ ALLEN CREEK NEAR ROCHESTER NY    : chr [1:2] "sfc_POLYGON" "sfc"
#>  $ NORTHRUP CREEK AT NORTH GREECE NY: chr [1:2] "sfc_POLYGON" "sfc"
#>  $ EAST BROOK EAST OF WALTON NY     : chr [1:2] "sfc_POLYGON" "sfc"
#>  $ BLACK RIVER AT WATERTOWN NY      : chr [1:2] "sfc_POLYGON" "sfc"

# convert to SpatVector, check validity
l_vect <- lapply(l3, vect)
lapply(l_vect, is.valid) |> str()
#> List of 6
#>  $ BLACK CREEK AT CHURCHVILLE NY    : logi TRUE
#>  $ OATKA CREEK AT GARBUTT NY        : logi TRUE
#>  $ ALLEN CREEK NEAR ROCHESTER NY    : logi FALSE
#>  $ NORTHRUP CREEK AT NORTH GREECE NY: logi TRUE
#>  $ EAST BROOK EAST OF WALTON NY     : logi TRUE
#>  $ BLACK RIVER AT WATERTOWN NY      : logi TRUE

# lets deal with that as well 
l_vect <- lapply(l_vect, makeValid) 

# check resulting geomtypes
lapply(l_vect, geomtype) |> str()
#> List of 6
#>  $ BLACK CREEK AT CHURCHVILLE NY    : chr "polygons"
#>  $ OATKA CREEK AT GARBUTT NY        : chr "polygons"
#>  $ ALLEN CREEK NEAR ROCHESTER NY    : chr "polygons"
#>  $ NORTHRUP CREEK AT NORTH GREECE NY: chr "polygons"
#>  $ EAST BROOK EAST OF WALTON NY     : chr "polygons"
#>  $ BLACK RIVER AT WATERTOWN NY      : chr "polygons"

# plot
par(mfrow=c(2,3))
lapply(names(l_vect), \(n) plot(l_vect[[n]], main = n))

^{Created on 2023-07-18 with reprex v2.0.2}