5

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

enter image description here

The grids are constructed as in the following picture.

enter image description here

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.

zx8754
  • 52,746
  • 12
  • 114
  • 209
Janak
  • 653
  • 7
  • 25

3 Answers3

3

One way of solving this problem is to compute the alpha hull around your points.

The alphahull package can do this. The package has excellent documentation, complete with animation, at http://yihui.name/en/2010/04/alphahull-an-r-package-for-alpha-convex-hull/

It is well worth taking a look at this documentation, in particular to understand the meaning of the alpha parameter.

First, replicate your data

A <- read.table(header=TRUE, text="
lon    lat
[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")

Now compute and plot the alpha hull. You have to provide the alpha value. I had to experiment to find a value small enough to capture the ouline.

library("alphahull")
hull <- with(A, ahull(lat, lon, alpha=0.033))
plot(hull)

enter image description here

Andrie
  • 176,377
  • 47
  • 447
  • 496
  • Thanks for your attention and efforts to answer my question. Since my data points are not exactly equi-spaced and shape of the point set are very irregular , I observed that "allphahull" frequently missed some boundary points. – Janak Feb 16 '15 at 10:50
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:

enter image description here

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.

Spacedman
  • 92,590
  • 12
  • 140
  • 224
  • You are right in pointing out that the grid fragments with hole or grid cells are only touching at a single corner point cannot occur in my data set due to some constraints in previous data processing stage. I have used your code for lots of samples. It is working fine. Grateful to you for your help. Thanks again. – Janak Feb 17 '15 at 05:43
  • I'm pretty sure there's an easier way. 1) drop all vertices which have 4 neighbours and the sum of the degree of the neighbours is >=12. 2) drop any edges where both vertices have degree 3. That replaces all the isomorphism stuff. 1) is easy, an easy method for 2) is eluding me... – Spacedman Feb 17 '15 at 17:03
  • I just gave a solution based on the Delaunay triangulation. Though I didn't include code. Seems a lot like what you are trying to suggest in your comment. – Nuclearman Feb 19 '15 at 02:39
  • Hi Spacedman. I am here again. Now I am trying to get the contour polygon and the hole inside the polygon is permitted. Could you please suggest the modification to the exiting code to accommodate the hole. – Janak May 24 '16 at 14:58
  • @Janak make a new question with this and try and give some sample data that illustrates the new problem. Also show how my code solution for this problem fails for that problem... – Spacedman May 24 '16 at 15:13
  • Sure. Here it is: http://stackoverflow.com/questions/37430604/extracting-the-outer-boundary-of-a-set-of-grid-points-with-some-missing-grids-h – Janak May 25 '16 at 07:38
1

Since the point set is a grid, or at least very close to a grid, seems like a good approach would be to make use of the Delaunay triangulation.

  1. Find the Delaunay triangulation.
  2. Find the two valid edge angles. If the grid was axis aligned, then these would be vertical and horizontal lines.
  3. Remove all edges that exceed error tolerances to be a valid edge on the grid based on the angle of the edge. An error tolerance of a few degrees is probably enough.
  4. Do a walk around remaining edges to find the perimeter. The left-most or right-most point is a good place to start the walk.

Performance should be O(N log N).

Nuclearman
  • 5,029
  • 1
  • 19
  • 35