1

我有一组保存为 SpatVector 的多边形，我正在尝试使用 terra 包计算多边形的面积。问题是区域的数量少于多边形。这种行为是不可预测的，并且随机发生在远程机器上！

#A reproducible example
 
library(sf)
library(sfheaders)
library(terra)

polygons_df<-data.frame(geom=c(1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,5,5,
                                    5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,8,8,8,8,
                                    9,9,9,9,9,10,10,10,10,10,10,10,10,10,10,10,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,
                                    11,11,11,11,11,11,11,11,11,11,11,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,13,13,
                                    13,13,13,14,14,14,14,14,14,14,15,15,15,15,15,15,15,15,15,15,15,15,15),
                                    part=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
                                                    1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
                                                    1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
                                                    1,1,1),
                                    x=c(22.2,22.3,22.3,21.9,21.9,22.0,22.0,22.1,22.1,22.2,22.2,23.6,23.9,23.9,23.7,23.7,23.6,23.6,25.3,25.4,25.4,25.3,25.3,25.5,25.6,
                                            25.6,25.5,25.5,25.4,25.5,25.5,25.4,25.4,20.1,20.2,20.2,20.1,20.1,27.4,27.6,27.6,27.5,27.5,27.4,27.4,27.8,28.2,28.2,27.8,27.8,
                                            27.7,27.7,27.9,27.9,27.8,27.8,30.7,31.0,31.0,31.3,31.3,31.4,31.4,31.1,31.1,30.9,30.9,30.8,30.8,30.7,30.7,30.6,30.6,30.7,30.7,
                                            20.6,20.9,20.9,20.6,20.6,25.4,25.5,25.5,25.4,25.4,25.7,25.9,25.9,26.1,26.1,26.0,26.0,25.9,25.9,25.7,25.7,29.9,30.1,30.1,30.3,
                                            30.3,30.6,30.6,31.2,31.2,31.4,31.4,31.5,31.5,31.8,31.8,31.7,31.7,31.5,31.5,31.0,31.0,30.9,30.9,30.8,30.8,30.7,30.7,30.1,30.1,
                                            29.7,29.7,29.8,29.8,29.9,29.9,25.8,25.9,25.9,26.2,26.2,26.4,26.4,26.3,26.3,26.2,26.2,26.0,26.0,25.9,25.9,25.8,25.8,25.7,25.7,
                                            25.5,25.5,25.6,25.6,25.7,25.7,25.8,25.8,28.5,28.7,28.7,28.5,28.5,25.6,25.9,25.9,25.7,25.7,25.6,25.6,25.5,26.0,26.0,25.8,25.8,
                                            25.6,25.6,25.5,25.5,25.4,25.4,25.5,25.5),
                                    y=c(34.0,34.0,33.6,33.6,33.7,33.7,33.8,33.8,33.9,33.9,34.0,35.0,35.0,34.8,34.8,34.9,34.9,35.0,35.5,35.5,35.4,35.4,35.5,35.5,35.5,
                                            35.4,35.4,35.5,35.4,35.4,35.3,35.3,35.4,36.8,36.8,36.6,36.6,36.8,37.0,37.0,36.9,36.9,36.8,36.8,37.0,37.2,37.2,36.9,36.9,36.8,
                                            36.8,37.0,37.0,37.1,37.1,37.2,37.7,37.7,37.6,37.6,37.4,37.4,37.3,37.3,37.2,37.2,37.3,37.3,37.4,37.4,37.5,37.5,37.6,37.6,37.7,
                                            37.6,37.6,37.5,37.5,37.6,37.6,37.6,37.5,37.5,37.6,38.1,38.1,38.0,38.0,37.9,37.9,37.8,37.8,37.9,37.9,38.1,38.8,38.8,38.7,38.7,
                                            38.6,38.6,38.5,38.5,38.4,38.4,38.5,38.5,38.4,38.4,38.3,38.3,38.2,38.2,38.1,38.1,38.0,38.0,38.1,38.1,38.2,38.2,38.3,38.3,38.4,
                                            38.4,38.7,38.7,38.6,38.6,38.8,39.1,39.1,39.0,39.0,38.9,38.9,38.8,38.8,38.7,38.7,38.8,38.8,38.5,38.5,38.3,38.3,38.2,38.2,38.3,
                                            38.3,38.5,38.5,38.7,38.7,39.0,39.0,39.1,38.4,38.4,38.3,38.3,38.4,39.4,39.4,39.2,39.2,39.3,39.3,39.4,40.0,40.0,39.5,39.5,39.6,
                                            39.6,39.7,39.7,39.8,39.8,39.9,39.9,40.0))

polygons_df$hole<-rep(0,nrow(polygons_df))

polygons_cc <- sfheaders::sf_polygon(
    obj = polygons_df
    , x = "x"
    , y = "y"
    , polygon_id = "geom"
)
sf::st_crs(polygons_cc) <- "+proj=longlat +datum=WGS84 +no_defs"

polygons_cc<-vect(polygons_cc)

plot(polygons_cc)

t<- 5005000

#prepare results data frame 
result_df<-data.frame(timestep=rep(t,length(terra::values(polygons_cc))),
                                            polygon_name=paste0(sprintf("%.0f",t),"_",
                                                                                    as.vector(unlist(terra::values(polygons_cc)))),
                                            geom=unname(terra::values(polygons_cc)))

发生错误的地方：


> # calculate area of polygons / objects 
> result_df$area<-terra::project(polygons_cc,"+proj=laea +lon_0=45 +lat_0=20 +datum=WGS84 +units=m +no_defs")%>%
+   expanse(unit="km")

> Error in `$<-.data.frame`(`*tmp*`, "area", value = c(1027.28688050304,  : 
  replacement has 15 rows, data has 18

> #another error from different shapefile: 
> Error in `$<-.data.frame`(`*tmp*`, "area", value = c(123.089596546252,  : 
  replacement has 78 rows, data has 81

以下是机器信息：

> sessionInfo()
R version 4.1.1 (2021-08-10)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server 7.9 (Maipo)

Matrix products: default
BLAS/LAPACK: /home/user/.conda/envs/snowflakes/lib/libopenblasp-r0.3.18.so

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=en_GB.UTF-8    
 [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=en_GB.UTF-8   
 [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] data.table_1.14.2     reticulate_1.22       Rfast_2.0.3          
 [4] RcppZiggurat_0.1.6    Rcpp_1.0.7            REdaS_0.9.3          
 [7] MASS_7.3-54           terra_1.4-22          dplyr_1.0.7          
[10] rgeos_0.5-8           ggplot2_3.3.5         stringr_1.4.0        
[13] sf_1.0-4              raster_3.5-2          sp_1.4-6             
[16] DescTools_0.99.44     spatstat_2.2-0        spatstat.linnet_2.3-0
[19] spatstat.core_2.3-2   rpart_4.1-15          nlme_3.1-153         
[22] spatstat.geom_2.3-0   spatstat.data_2.1-0   ncdf4.helpers_0.3-6  
[25] ncdf4_1.18           

loaded via a namespace (and not attached):
 [1] jsonlite_1.7.2        splines_4.1.1         assertthat_0.2.1     
 [4] expm_0.999-6          gld_2.6.3             lmom_2.8             
 [7] pillar_1.6.4          lattice_0.20-45       glue_1.5.0           
[10] polyclip_1.10-0       colorspace_2.0-2      Matrix_1.3-4         
[13] spatstat.sparse_2.0-0 pkgconfig_2.0.3       s2_1.0.7             
[16] purrr_0.3.4           mvtnorm_1.1-3         scales_1.1.1         
[19] tensor_1.5            rootSolve_1.8.2.3     spatstat.utils_2.2-0 
[22] tibble_3.1.6          proxy_0.4-26          mgcv_1.8-38          
[25] generics_0.1.1        ellipsis_0.3.2        withr_2.4.2          
[28] magrittr_2.0.1        crayon_1.4.2          deldir_1.0-6         
[31] fansi_0.4.2           class_7.3-19          tools_4.1.1          
[34] lifecycle_1.0.1       Exact_3.1             munsell_0.5.0        
[37] compiler_4.1.1        e1071_1.7-9           rlang_0.4.12         
[40] classInt_0.4-3        units_0.7-2           rstudioapi_0.13      
[43] goftest_1.2-3         boot_1.3-28           wk_0.5.0             
[46] gtable_0.3.0          codetools_0.2-18      abind_1.4-5          
[49] DBI_1.1.1             R6_2.5.1              utf8_1.2.2           
[52] KernSmooth_2.23-20    stringi_1.7.5         parallel_4.1.1       
[55] vctrs_0.3.8           png_0.1-7             tidyselect_1.1.1

我知道我可以sf单独在多边形上使用包或循环，但是，我有数百万个多边形，而且terra似乎只是在资源和计算时间方面有效。也许问题来自terra::values()？有没有办法expanse()使用多边形 ID 来检查问题出在哪里？

编辑：产生问题的shapefile之一是here

4

1 回答 1

0

很难说为什么会发生这种情况，因为我没有得到那个错误，所以它可能与使用不同版本的包有关（对于 sf，可能是 GEOS 库）。无论哪种方式，我认为如果你这样做会很好

v <- vect(as.matrix(polygons_df), "polygons", crs="+proj=longlat")
p <- project(v, "+proj=laea +lon_0=45 +lat_0=20 +datum=WGS84 +units=m")
expanse(p, unit="km") |> length()

我相信当您第一次使用 . 创建 sf 对象时会出现错误sf_polygon。第三个几何有三个部分

polygons_df[polygons_df$geom==3, "part"] |> unique() |> length()
#[1] 3

但是（的 sf 版本）polygon_cc不是多面体

polygons_cc[3,]
#Simple feature collection with 1 feature and 1 field
#Geometry type: POLYGON
#Dimension:     XY
#Bounding box:  xmin: 25.3 ymin: 35.3 xmax: 25.6 ymax: 35.5
#CRS:           +proj=longlat +datum=WGS84 +no_defs
#  geom                       geometry
#3    3 POLYGON ((25.3 35.5, 25.4 3...

而且看起来不太好：

在此处输入图像描述

在当前工作流程中，您只需跳过该步骤，如上所示。

代码中另一个不必要的部分是将多边形投影到lcc（Lambert 保形圆锥），然后计算它们的面积。您可以并且通常应该使用经度/纬度数据进行计算。即使在这种情况下，差异非常小。实际上， 的默认行为terra是将数据转换回 longlat 以计算面积。

于 2022-01-18T16:34:34.350 回答