R中多边形的长度

Question

我想计算每个多边形的长度。 - 在每个多边形周围我创建了点（st_sample）， - 从点的组合我创建了所有可能的折线， - 对于多边形内的折线，我计算了长度， - 最长的折线是我的结果（多边形的最大长度）。

我编写的代码让我得到了结果，但它真的很慢。你有一些改进我的代码的解决方案吗？我知道有两个循环我不能指望速度上的奇迹，但我不知道如何以另一种方式获得结果。

如果没有别的，我至少有一些替代解决方案可以在没有循环的情况下一步一步从一个多边形的点组合创建所有折线？ :)

谢谢

library(sf)
library(data.table)

poly=st_read(system.file("shape/nc.shp", package="sf"))
poly=poly[1:10,]
poly=st_cast(poly,"POLYGON")
poly$max_length=0

##Combination of 10 points, withot repetiton
aa=CJ(1:10,1:10)
aa=aa[!duplicated(t(apply(aa[,.(V1, V2)], 1, sort))),][V1!=V2]

##for each polygon create sample of coordinates along line, from them I create polyline and calculated length for linestring which are inside polygon

for (ii in 1:nrow(poly)){
  ncl=st_cast(poly[ii,],"LINESTRING")
  ##sample of point along line
  ncp=st_cast(st_sample(ncl,10, type="regular", exact=T),"POINT")
  
  ##create empty sf 
  aaa=st_sf(st_sfc())
  st_crs(aaa)="NAD27"
  
  ##for each combination of points create linestring and calculate length only for polylines which are inside polygon
  for (i in 1:nrow(aa)){
 aaa=rbind(aaa,st_sf(geometry=st_cast(st_union(ncp[t(aa[i])]),"LINESTRING")))
  }
poly$max_length[ii]=as.numeric(max(st_length(aaa[unlist(st_contains(poly[ii,],aaa)),])))
}

第二次尝试在 data.table 中运行函数。少一个循环，但问题可能是第二个循环。

poly=st_read(system.file("shape/nc.shp", package="sf"))
poly=poly[1:10,]
poly=st_cast(poly,"POLYGON")
poly$max_length=0

##Combination of 10 points, withot repetiton
aa=CJ(1:10,1:10)
aa=aa[!duplicated(t(apply(aa[,.(V1, V2)], 1, sort))),][V1!=V2]

overFun <- function(x){
  
  ncl=st_cast(x[,geometry],"LINESTRING")
  ##sample of point along line
  ncp=st_cast(st_sample(ncl,40, type="regular", exact=T),"POINT")
  
  ##create empty sf 
  aaa=st_sf(st_sfc())
  st_crs(aaa)="NAD27"
  
  ##for each combination pof points create linestring and calculate length
  for (i in 1:nrow(aa)){
    aaa=rbind(aaa,st_sf(geometry=st_cast(st_union(ncp[t(aa[i])]),"LINESTRING")))
  }
  
  as.numeric(max(st_length(aaa[unlist(st_contains(x[,geometry],aaa)),])))}  

setDT(poly)

##run function inside data.table
poly[,max_length:=overFun(poly), by=seq(nrow(poly))]

编辑：我为我的问题找到了一些足够快的解决方案来满足我的需求。 在 data.table 中使用并行库，其功能也适用于 data.table。仍然存在问题，为什么使用函数 st_contains 排除了某些折线（见上图）。也许精度有问题？

library(sf)
library(data.table)

poly=st_read(system.file("shape/nc.shp", package="sf"))
poly=st_cast(poly,"POLYGON")
setDT(poly)

##Combination of 10 points, withot repetiton
aa=CJ(1:10,1:10)
aa=aa[!duplicated(t(apply(aa[,.(V1, V2)], 1, sort))),][V1!=V2]

overFun <- function(x){
  
ncl=st_cast(poly[1,geometry],"LINESTRING")
##sample of point along line
ncp=st_cast(st_sample(ncl,10, type="regular", exact=T),"POINT")
  
df=data.table(ncp[aa[,V1]],ncp[aa[,V2]] )
  
df[,v3:=st_cast(st_union(st_as_sf(V1),st_as_sf(V2)),"LINESTRING"), by=seq(nrow(df))]
as.numeric(max(st_length(df[unlist(st_contains(poly[1,geometry], df$v3)),]$v3)))}

library(parallel)
cl <- makeCluster(detectCores() - 1)
clusterExport(cl, list("overFun","data.table","st_cast","CJ","poly","st_sample","st_sf","st_sfc","aa","st_length","st_union",
                       "st_as_sf","st_contains"))
system.time(poly[,c("max_length"):=.(clusterMap(cl, overFun, poly$geometry)),])
stopCluster(cl)

Answer 1

如果您追求多边形的周长，请考虑以下代码：

library(sf)
library(dplyr)

shape <- st_read(system.file("shape/nc.shp", package="sf")) # included with sf package

lengths <- shape %>% 
  mutate(circumference = st_length(.)) %>% 
  st_drop_geometry() %>% 
  select(NAME, circumference) 
   
head(lengths)   
         NAME circumference
1        Ashe  141665.4 [m]
2   Alleghany  119929.0 [m]
3       Surry  160497.7 [m]
4   Currituck  301515.3 [m]
5 Northampton  211953.8 [m]
6    Hertford  160892.0 [m]

如果您的内部有一些孔并且不希望它们包含在圆周中，请考虑通过nngeo::st_remove_holes() 将它们移除。

Answer 2

我也遇到过类似的问题，坦白说还没有找到现成的解决方案。

我将使用 sf 包中的同一个阿什县。

library(sf)
library(dplyr)

shape <- st_read(system.file("shape/nc.shp", package="sf")) %>% 
  dplyr::filter(CNTY_ID == 1825) %>% # Keep only one polygon
  st_transform(32617) # Reproject to WGS 84 / UTM zone 17N

解决方案 1

您可以仅使用dplyr、tidyr 和sf 将多边形转换为点并计算所有点之间的距离。从这个品种中，选择最大值。这将是您的示例图中的一条绿线。

library(tidyr)

shape %>% 
  st_cast("POINT") %>% # turn polygon into points
  distinct() %>% # remove duplicates
  st_distance() %>% # calculate distance matrix
  as.data.frame() %>% 
  gather(point_id, dist) %>% # convert to long format
  pull(dist) %>% # keep only distance column
  max()
#> Warning in st_cast.sf(., "POINT"): repeating attributes for all sub-geometries
#> for which they may not be constant
#> 45865.15 [m]

解决方案 2

您也可以使用Momocs 包。它是为 2D 形态测量分析而创建的。虽然在第一种情况下不必将我们的形状重新投影到 UTM（sf 可以处理地理坐标），但您的多边形应该在 Momocs 包的情况下进行投影。

library(Momocs)

shape %>% 
  st_cast("POINT") %>% # Polygon to points
  distinct() %>% # remove duplicates
  st_coordinates() %>% # get coordinates matrix
  coo_calliper() # calculate max length
#> Warning in st_cast.sf(., "POINT"): repeating attributes for all sub-geometries
#> for which they may not be constant
#> [1] 45865.15

评论

Momocs 包中还有其他几个功能。例如，您可以根据形状的惯性轴计算形状的长度，即与 x 轴对齐。 coo_length 将返回给您44432.02 [m]。

例如，可以将Momocs 包中的多个函数应用于坐标矩阵，如下所示：

point_matrix <- shape %>% 
  st_cast("POINT") %>% 
  distinct() %>% 
  st_coordinates() 
#> Warning in st_cast.sf(., "POINT"): repeating attributes for all sub-geometries
#> for which they may not be constant

funs <- list("length" = coo_length,
             "width" = coo_width,
             "elongation" = coo_elongation)

sapply(funs, function(fun, x) fun(x), x = point_matrix)
#>       length        width   elongation 
#> 4.443202e+04 3.921162e+04 1.174917e-01