Find the perimeter of a subset of a near-regular grid of points

Question

Let us consider a set of near-regular grids in 2-D. These grids are adjacent (neighbouring grids have one or more same vertices) to the neighbouring grids. Here are the sample of 10 grids with the coordinates of the vertices (longitude,latitude) are as follows

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

The scatter-plot of the above sample set of points (vertices):

The grids are constructed as in the following picture.

My question is how to get the perimeter (red contour passing through all boundary points)??

Note that: The red encircled points (1,3,7,9,10,13,17,18,19,16,15,12,11,6,5,2) in figure 1 are the boundary points.

Note: It is observed the sides of the grids are less than 6000 metres and length of diagonals of all grids are more than 6000 metres.

I am using distHaversine from the geosphere package function in R to find the distance between two points.

Answer 1

In outline: all pairs of points closer than 6000m form a graph in the form of grid squares. Construct that graph, and then find all subgraphs isomorphic to a square (a loop of size 4). The external edges will appear less often than internal edges since they are only part of one square (internal edges are shared by multiple squares). Hence find all the internal edges and drop them, then traverse the resulting graph which should be a simple loop.

Code:

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)
}

Usage:

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

Results:

Warnings:

This is untested on grid fragments with holes or where grid cells are only touching at a single corner points. Maybe that can't happen. Also I'm not sure what the efficiency of the algorithms are, but it seems pretty quick.