查找点的近似规则网格的子集的周长

Question

让我们考虑二维的一组近似规则的网格。这些网格与相邻网格相邻（相邻网格具有一个或多个相同的顶点）。这是10个网格的示例，其顶点坐标（经度，纬度）如下

A<-

        lon    lat
        [,1]     [,2]
  [1,] 85.30754 27.91250
  [2,] 85.32862 27.95735
  [3,] 85.34622 27.89880
  [4,] 85.36732 27.94364
  [5,] 85.34958 28.00202
  [6,] 85.38831 27.98830
  [7,] 85.38487 27.88508
  [8,] 85.40598 27.92991
  [9,] 85.42353 27.87134
 [10,] 85.44466 27.91616
 [11,] 85.42698 27.97456
 [12,] 85.46567 27.96081
 [13,] 85.48334 27.90239
 [14,] 85.50437 27.94703
 [15,] 85.48645 28.00502
 [16,] 85.52517 27.99123
 [17,] 85.52198 27.88862
 [18,] 85.54302 27.93325
 [19,] 85.56384 27.97745

以上样本点集（顶点）的散点图：

网格的结构如下图所示。

我的问题是如何获得周长（红色轮廓线穿过所有边界点）？

注意：图1中的红色圆圈点（1、3、7、9、10、13、17、18、19、16、15、12、11、6、5、2）是边界点。

注意：观察到网格的侧面小于6000米，所有网格的对角线长度均大于6000米。

我正在使用R中distHaversine的geospherepackage函数来查找两点之间的距离。

Answer 1

概述：所有距离小于 6000m 的点对以网格方块的形式形成一个图形。构造该图，然后找到与正方形同构的所有子图（大小为 4 的循环）。外部边缘比内部边缘出现的频率要低，因为它们只是一个方块的一部分（内部边缘由多个方块共享）。因此，找到所有内部边并删除它们，然后遍历生成的图，这应该是一个简单的循环。

代码：

library(igraph); library(geosphere)

# main function
perimeterGrid <- function(pts, maxdist=6000, mindist=1){
    g = edgeP(makegrid(pts, maxdist=maxdist, mindist=mindist))
    loop = graph.dfs(minimum.spanning.tree(g),1)$order
    cbind(V(g)$x, V(g)$y)[loop,]
}

# haversine distance matrix
dmat <- function(pts){
    n=nrow(pts)
    do.call(rbind,lapply(1:n,function(i){distHaversine(pts[i,],pts)}))
}

# make the grid cells given a maxdist (and a mindist to stop self-self edges)    
makegrid <- function(pts, maxdist=6000,  mindist=1){
    dists = dmat(pts)
    g = graph.adjacency(dists<maxdist & dists>mindist,
        mode="undirected")
    V(g)$x=pts[,1]
    V(g)$y=pts[,2]
    g
}

# clever function that does the grid edge count etc
edgeP <- function(g){
    # find all the simple squares
    square=graph.ring(4)
    subs = graph.get.subisomorphisms.vf2(g,square)
    # expand all the edges
    subs = do.call(rbind, lapply(subs, function(s){
        rbind(s[1:2], s[2:3], s[3:4], s[c(4,1)])
    }))
    # make new graph of the edges of all the squares
    e = graph.edgelist(subs,directed=FALSE)
    # add the weight as the edge count
    E(e)$weight=count.multiple(e)

    # copy the coords from the source
    V(e)$x=V(g)$x
    V(e)$y=V(g)$y

    # remove multiple edges
    e=simplify(e)

    # internal edges now have weight 256.
    e = e - edges(which(E(e)$weight==256))
    # internal nodes how have degree 0
    e = e - vertices(degree(e)==0)
    return(e)
}

用法：

 plot(pts)
 polygon(perimeterGrid(pts),lwd=2)

结果：

警告：

这尚未在带有孔的网格片段或网格单元仅在单个角点处接触的网格片段上进行测试。也许这不可能发生。我也不确定算法的效率如何，但看起来相当快。