R中多边形的长度

问题描述

我想计算每个多边形的长度。 -围绕我创建的每个多边形点（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",][V1!=V2]

overFun <- function(x){
  
  ncl=st_cast(x[,geometry],40,"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,"LInesTRING")))
  }
  
  as.numeric(max(st_length(aaa[unlist(st_contains(x[,])))}  

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,][V1!=V2]

overFun <- function(x){
  
ncl=st_cast(poly[1,"LInesTRING")
##sample of point along line
ncp=st_cast(st_sample(ncl,"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,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)

解决方法

如果您在多边形的周围，请考虑以下代码：

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()移除它们。

data.table data.table r r sf