英文:
terra::vect does not work for larger polygons
问题
我正在尝试使用R包streamstats
划定的分水岭边界,将其转换为SpatVector
对象,使用terra::vect
。我将一个包含6个我正在处理的sf
对象分水岭子集的Rdata文件上传到Google Drive,以便制作一个可重现的示例,下面是显示列表中的最后一个分水岭(也比其他5个大得多)在vect
中不起作用的代码:
library(sf)
library(terra) # terra 1.6.17
## 导入数据:
load('C:/PhD/Research/Code/Manuscript/top6.Rdata')
## 这是一个与terra::vect一起使用的分水岭的示例:
# 将第5个分水岭转换为sf多边形对象:
v <- st_cast(st_boundary(l_DA_6[[5]]$geometry), "POLYGON")
# 绘制它:
plot(v) # 看起来不错
# 现在将sf对象转换为SpatVector:
SpatV<-vect(v)
# 并绘制SpatVector:
plot(SpatV) # 也看起来不错
## 这是在terra::vect中不起作用的分水岭:
# 将第6个分水岭转换为sf多边形对象:
v <- st_cast(st_boundary(l_DA_6[[6]]$geometry), "POLYGON")
# 绘制它:
plot(v) # 看起来不错
# 现在将sf对象转换为SpatVector:
SpatV<-vect(v)
# 并绘制SpatVector:
plot(SpatV) # 只显示一个矩形,这是错误的
我唯一能看出两个sf
对象之间不同的是第二个对象要大得多。我在CRAN版本的terra
和开发者版本的terra 1.6.17
上都遇到了这个问题。任何帮助都将不胜感激。
英文:
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.
答案1
得分: 3
以下是翻译好的部分:
"All sf
objects seem to be invalid, once passed through st_make_valid()
and after unifying classes, conversion to SpatVector
looks fine."
"一旦通过st_make_valid()
传递,然后统一类别,所有的sf
对象似乎都是无效的,转换为SpatVector
后看起来很好。"
"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."
"由于生成的SpatVectors最终变成了线而不是多边形,我进行了编辑(第1版)。这导致了一些额外的操作,丢弃了一些小的多边形,因此您可能应该验证这是否确实是可接受的结果。"
"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."
"这里潜在问题的一部分可能与未投影的CRS以及在加载sf
时使用s2有关(注意在加载sf
时sf_use_s2()
为TRUE)。因此,也许可以考虑在较早的步骤中更改为投影的CRS。"
(以下为代码部分,不提供翻译。)
英文:
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))
<!-- -->
<sup>Created on 2023-07-18 with reprex v2.0.2</sup>
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论