1

enter image description hereI am working on an application where a user can select a land plot (sf_POLYGON) by clicking on it. Once selected (highlighted plot no. 6 in the figure), I then need to calculate some "offsets" and plot the resulting polygon. These offsets are dependent on the adjacent polygons. For example, if the selected plot has a road adjacent to it, then the offset is say 6m inside the plot. The side adjacent to a plot will have an offset of 4m. The same is displayed in the inscribed polygon. The blue lines are for sides touching the road and the orange lines are for the sides touching other plots.

The way I was thinking to solve this was :

  1. Find the vertices of the selected polygon.
  2. Draw lines between the vertices.
  3. Use st_intersect or st_touches to find out the polygons that the line touches.

I tried to convert the polygon into LINESTRING using st_cast(test_poly, 'LINESTRING') to figure out if I can get the lines automatically and then see which polygons they intersect/touch. However st_coordinates on the polygon cast to LINESTRING reveals coordinates which are not all vertices. So my question is, how do I extract the vertices of the selected polygon? Or if there is a better approach to solve this problem?

The second problem is of creating the offset polygon. Not sure how to do that, but i am guessing buffer is something that would help?

library(sf)
library(leaflet)
library(dplyr)[![enter image description here][1]][1]
test_data <- readRDS('test_data.Rds')
leaflet() %>% 
  addTiles() %>% 
  addPolygons(data = test_data, fillOpacity = 0.8, stroke = TRUE,
              weight = 1,
              color = 'black',
              label = test_data$Plot_Numbe,
              layerId = test_data$OBJECTID,
              fillColor = test_data$color)

test_poly <- test_data %>% filter(OBJECTID == 1977)
test_line <- st_cast(test_poly, 'LINESTRING')
st_coordinates(test_line)

Can't seem to attach the data. Here's a link to it.

6
  • 1
    The points are all vertices, they are vertices where the angle is very close to (or equal to) 180 degrees. You can remove these with careful use of st_simplify, taking care not to remove "vertices" that you want to keep.
    – Spacedman
    Jun 22, 2020 at 11:50
  • @Spacedman thanks for this. Trying to figure out a good value for dTolerance for different polygons. I suppose I don't need to cast the polygon to a linestring, right?
    – Dhiraj
    Jun 22, 2020 at 12:44
  • No, st_coordinates will pull out polygon coordinates just fine. Check it works with MULTIPOLYGON objects in case you have any of those...
    – Spacedman
    Jun 22, 2020 at 13:40
  • Thanks, checking it now. Do you think the approach that I am hoping to implement is ok, or is there a better way of doing it?
    – Dhiraj
    Jun 22, 2020 at 13:44
  • You've asked two questions and we like one question per post here. I think it might be acceptible to ask something like how to create a variable width internal buffer in R as a new post.
    – Spacedman
    Jun 22, 2020 at 15:12

1 Answer 1

1

So I finally figured it out, all thanks to @Spacedman for showing me the way to find out my 'vertices' or perhaps corners would be a better terminology. Given the data I have, I wrote this function to identify not only the sides which border a road, but also identify corners, given that the road is one single polygon with no breaks. Struggled a lot to figure out how st_relate, st_intersects work and ended up learning a lot about DE-9IM. This might help someone facing similar issues with land plots/parcels:

find_neighbors <- function(id) {
  test_poly <- final %>% filter(OBJECTID == id)
  test_line <- st_cast(test_poly, 'MULTILINESTRING')
  x <- st_touches(test_line, final)
  neighbors <- final[x[[1]], ] %>% filter(OBJECTID != id)
  #Find if neighbor contains a road
  road_nbr <- which(neighbors$Landuse == 'Road')
  #Reduce the vertices of test_poly
  test_poly <- st_simplify(test_poly, dTolerance = 0.00001)
  # Make linestrings from all sides of the test polygon
  poly_coords <- st_coordinates(test_poly)
  vertices <- st_multipoint(poly_coords) %>% st_zm() %>% 
    st_geometry() %>% st_cast('POINT')
  # Number of total sides to be created
  sides <- length(vertices) - 1
  # Build linestrings
  linestrings <- lapply(X = 1:sides, FUN = function(x) {
    pair <- st_combine(c(vertices[x], vertices[x + 1]))
    line <- st_cast(pair, "LINESTRING") %>% st_set_crs(st_crs(neighbors))
    return(line)
  })
  multilinestring <- st_multilinestring(do.call("rbind", linestrings))
  #Find sides bordering roads if roads are in the neighbors list
  if(length(road_nbr) > 0) {
    intersects_road <- logical()
    for(i in 1:length(road_nbr)) {
      for(j in 1:length(linestrings)) {
        buffered_poly <- st_buffer(multilinestring[[j]], dist = 0.00001,
                                   singleSide = TRUE,
                                   endCapStyle = 'FLAT') 
        buffer_again <- st_buffer(buffered_poly, dist = -0.000001,
                                  singleSide = TRUE,
                                  endCapStyle = 'FLAT')
        z <- st_intersects(buffer_again, neighbors[road_nbr[1],])
        intersects_road <- c(intersects_road, !is_empty(z[[1]]))
      }
    }
    is_corner <- with(rle(intersects_road), sum(lengths[values] >= 2))
    return(list(intersects_road, is_corner))
  }
}

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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