1

I have a point and a square polygon. Both are sp classes.

The point intersects one side of the polygon, and I need to know which side of the square it intersects (i.e. north, east, south, or west). The square almost lines up with the compass, so there are clear n, e, s, w edges, but it doesn't line up perfectly.

The diagram illustrates this; in this case the point intersects the northern edge.

Example of the problem - I need to know which edge the point intersects. In this case it intersects the northern edge

Any ideas how to work out which edge it intersects? On the surface it seems like a simple problem, but a solution has eluded me.

  • compare p.x with poly's bbox minx & maxx repeat for y – Ian Turton Feb 23 '17 at 11:46
  • 2
    Sample data would be nice, save us all having to set it up ourselves... – Spacedman Feb 23 '17 at 14:02
  • @iant That works only for perfectly isothetic squares, whereas here the OP has been careful to indicate that this is not going to be exactly the case. One should also be concerned that the point might not exactly intersect the square, making it possible for your proposed comparison to indicate the point lies on two sides. A good solution accommodate such possibilities. – whuber Feb 23 '17 at 17:17
  • 1
    what if it is a corner point? – Edzer Pebesma Feb 23 '17 at 21:33
  • 1
    does the square polygon consist of five points, and is it predictable where it starts and which direction it has? – Edzer Pebesma Feb 23 '17 at 21:35
3

Here's a solution - works on the assumption that your bounding boxes are square shaped in a projected CRS. Split into 2 chunks for clarity.

Data preparation

library(stplanr) # load sp + line_midpoint fun
d = SpatialPoints(coords = matrix(rnorm(100), ncol = 2)) # test data
b = bb2poly(bb = d) # create polygon of bb
p = raster::geom(b) # extract vertices of polygon
for(i in 1:4){ # the tricky bit - split into 4 lines
  if(i == 1)
    l = raster::spLines(rbind(p[i, c("x", "y")], p[i + 1, c("x", "y")])) else
      l = raster::bind(
        l,
        raster::spLines(rbind(p[i, c("x", "y")], p[i + 1, c("x", "y")]))
      )
}
l_points = line_midpoint(l) # midpoint - makes the nearest problem easier (only works for squares - use something else for rectangles)

Finding the edge + viz

plot(l)
points(l_points)
l_points$bearing = c("w", "n", "e", "s") # add names to edges
text(x = l_points@coords[,1] + 0.1, y = l_points@coords[,2] + 0.1, l_points$bearing)
p_to_detect = d[5,]
plot(p_to_detect, add = T)
(nearest_one = nabor::knn(data = coordinates(l_points), query = coordinates(p_to_detect), k = 1))
nearest_side = l_points$bearing[nearest_one$nn.idx]
text(x = p_to_detect@coords[,1] + 0.1, y = p_to_detect@coords[,2] + 0.1, nearest_side)

Result:

enter image description here

Discussion

I'd like to see this implemented in sf:

https://github.com/edzer/sfr

There are probably more efficient ways so more answers welcome (don't accept this as the right answer just yet!)

  • Thanks RobinLovelace, this seems to work well. My only query is whether it's safe to assume that the four sides of the polygon will be in the order west-north-east-south? (@EdzerPebesma queried this too). If I were making my own polygons this would not be a problem, but I can envisage circumstances when they will be read in from a shapefile with readOGR(). – Nick Malleson Feb 27 '17 at 11:28
  • Should be easy to deal with: just find the line whose midpoint is furthest north and then the one that's furthest east etc. – RobinLovelace Feb 28 '17 at 16:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.