0

I have a polygon that I would like to split into two using lwgeom::st_split(). The polygon is a buffer that I created using st_buffer() from the sf package around a linestring, which is also the "blade" object with which I want to split the polygon.

The difficulty I am having is that the line does not fully pass through the polygon, which prevents successful splitting (see the reprex below). This follows from the buffer being made in the round end cap style, which is a choice that I am forced to make for my particular application.

Is there a way to successfully split a polygon with a linestring that does not fully pass through the polygon? Such as by first cutting off the round ends somehow? Or some other operation?

library(sf)
library(tmap)

## create the linestring 
line <- structure(list(structure(c(-792285.491251808,-792285.488817639,-792285.436522982,-792285.384427687,-792285.33267453,-792285.281405349,-792285.230760659,-792285.180879259,-792285.13189786,-792285.083950704,-792285.037169198,-792284.991681557,-790240.702683605,-790240.659579577,-750536.528415538,-743907.032905619,-743907.021055776,-705997.347638523,-704336.097999047,-704336.060652765,-664913.220854293,-659976.689106017,-659976.677471231,-640303.816873962,-640303.813512255,-620718.901662572,-616828.918835104,-610882.403771673,-610882.394685652,-595564.578577553,-578300.598930178,-576658.107865964,-576658.074277915,-537693.254220241,-534429.766727421,-534429.748408387,-518012.54009363,-485910.696784914,-485910.691773335,-452102.325903294,-448880.542010921,-448880.520954105,-442063.307672312,-442063.266728173,-442063.227500371,-442063.190094792,-442063.154612405,-442063.121148988,-442063.089794869,-442063.060634682,-442063.033747138,-442063.009204815,-442062.98707396,-442062.96741431,-442062.950278934,-442062.935714084,-442062.923759075,-442062.914446178,-442062.90780053,-442062.90384007,-442062.902575489,-442062.9040102,-442062.90814033,-442062.914954731,-442062.924435009,-442062.936555573,-442062.951283708,-442062.968579657,-442062.988396733,-442063.010681444,-442063.035373637,-442063.062406661,-442063.091707545,-442063.096064972,2918748.30210906, 2918748.30212301, 2918748.2996829, 2918748.29450937,2918748.28661659, 2918748.2760262, 2918748.26276722, 2918748.24687599,2918748.22839606, 2918748.20737808, 2918748.18387966, 2918748.15796519,2917511.54289283, 2917511.51528765, 2890626.4525922, 2886229.1088279,2886229.10084542, 2860296.14740504, 2859255.91763225, 2859255.89306825,2832044.60830091, 2828717.72270611, 2828717.71474522, 2815052.34456301,2815052.34221764, 2801328.70929097, 2798649.81543859, 2794555.39323204,2794555.38690149, 2783756.45528285, 2771592.78589961, 2770523.36937069,2770523.34651698, 2742836.82623387, 2740608.56736203, 2740608.55455154,2728852.61822133, 2705882.11145883, 2705882.10784934, 2681373.37633357,2679115.88414793, 2679115.86898237, 2674070.35896075, 2674070.32697836,2674070.29291262, 2674070.25685546, 2674070.21890423, 2674070.17916136,2674070.13773414, 2674070.09473438, 2674070.05027816, 2674070.00448547,2674069.95747993, 2674069.90938842, 2674069.86034074, 2674069.8104693,2674069.75990871, 2674069.70879544, 2674069.65726748, 2674069.6054639,2674069.55352455, 2674069.50158962, 2674069.44979929, 2674069.39829337,2674069.34721089, 2674069.29668973, 2674069.24686627, 2674069.19787498,2674069.14984812, 2674069.10291533, 2674069.05720328, 2674069.01283537,2674068.96993136, 2674068.96421305), dim = c(74L, 2L), class = c("XY","LINESTRING", "sfg"))), class = c("sfc_LINESTRING", "sfc"), precision = 0, bbox = structure(c(xmin = -792285.491251808,ymin = 2674068.96421305, xmax = -442062.902575489, ymax = 2918748.30212301), class = "bbox"), crs = structure(list(input = "+proj=robin +ellps=WGS84 +datum=WGS84 +units=m +no_defs", 
wkt = "PROJCRS[\"unknown\",\n    BASEGEOGCRS[\"unknown\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ID[\"EPSG\",6326]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8901]]],\n    CONVERSION[\"unknown\",\n        METHOD[\"Robinson\"],\n        PARAMETER[\"Longitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"False easting\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"), class = "crs"), n_empty = 0L)

## create the buffer
buff <- st_buffer(line, 80467.2, endCapStyle = "ROUND")

## the line does not fully pass through the buffer polygon
tm_shape(buff) +
  tm_fill() +
  tm_shape(line) +
  tm_lines()


## thus, splitting the polygon does not successfully create two separate polygons
split <- buff %>%
  lwgeom::st_split(line) %>%
  st_collection_extract("POLYGON") 

tm_shape(split) +
  tm_fill("coral") +
  tm_shape(line) +
  tm_lines()

Improve this question
4
  • How would you expect the spit to work? Given the round end caps there is no such thing as the closest point to the line's end... Flat caps would side step the issue, as you seem to be aware
    – Jindra Lacko
    Commented Oct 12, 2023 at 16:24
  • I wondered if there was a pre-split step to make the overall splitting process work. Flat end caps would avoid this issue, but with my specific data this runs into the problem of strange buffer shapes. This is discussed in my as-yet unresolved post [gis.stackexchange.com/questions/468208/….
    – jccborders
    Commented Oct 12, 2023 at 16:33
  • You seem to want to extend the ends in the general direction of the line feature. You could do this by taking the first and last point (as in an Answer below) but this would fail in many cases, for example if your line was a "C" shape. Are all your features similar to your example here? ie nearly collinear vertices with no bends in the feature?
    – Spacedman
    Commented Oct 14, 2023 at 6:54
  • The linestrings take a variety of shapes, including "C", in ways that would be problematic as you suggest. However, I came up with a way to sidestep the original issue while using round end caps on the buffers. Thanks for the help.
    – jccborders
    Commented Oct 23, 2023 at 14:50

1 Answer 1

Reset to default
0

Just extend the line out from both ends and split

# Extract the first and last points as numeric values
first_point <- as.numeric(line[[1]][1,])
last_point <- as.numeric(line[[1]][nrow(line[[1]]),])

# Compute the overall direction vector
direction <- last_point - first_point
direction <- direction / sqrt(sum(direction^2))

# Extend the points x/2 on each side
len <- as.numeric(st_length(line)[1])
new_first_point <- first_point - (len/2) * direction
new_last_point <- last_point + (len/2) * direction

# Create the new extended line
extended_line <- st_sfc(st_linestring(rbind(new_first_point, first_point, last_point, new_last_point)))
original_crs <- st_crs(line)
extended_line <- st_set_crs(extended_line, original_crs)

split <- buff %>%
  lwgeom::st_split(extended_line) %>%
  st_collection_extract("POLYGON") 

tm_shape(split[1]) +
  tm_fill("coral") +
tm_shape(split[2]) +
  tm_fill("cyan") +
tm_shape(extended_line) +
  tm_lines()

enter image description here

Improve this answer
1
  • This approach does not perfectly help in my application, because I need to use the linestrings as they originally are (without extension). However, I came up with a way to sidestep the original issue while using round end caps on the buffers. Thanks for the help.
    – jccborders
    Commented Oct 23, 2023 at 14:51

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.