# Finding the side of a square that a point intersects with R (sp)

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.

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
• 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
• what if it is a corner point? – Edzer Pebesma Feb 23 '17 at 21:33
• 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

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

## 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