I have a kml file that describes four polygons(the kml file is called "core.kml", and I have put a copy at this link). I would like to choose random points within the area defined by the polygons. Ideally I would like to choose points evenly over the whole area, so choose more points from the larger polygons than the small ones, in proportion to their area. I have borrowed the code from answer 1 here, which works beautifully with the North Carolina polygon (at least when I update it by inserting nc <- st_geometry(nc)). But with my polygons, the st_sample step always gives this error:

st_as_s2(): dropping Z and/or M coordinate Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, :
Loop 0 is not valid: Edge 15 crosses edge 19

I have tried inserting sf_use_s2(FALSE), but I still get the same error. The plot command does work, and the polygons look as expected. My code:

core<-st_read("core.kml")
core<-st_geometry(core) 
sf_use_s2(FALSE) 
core.points<-st_sample(core,5) 
plot(core)

I have a kml file that describes four polygons(the kml file is called "core.kml", and I have put a copy at this link). I would like to choose random points within the area defined by the polygons. Ideally I would like to choose points evenly over the whole area, so choose more points from the larger polygons than the small ones, in proportion to their area. I have borrowed the code from answer 1 here, which works beautifully with the North Carolina polygon (at least when I update it by inserting nc <- st_geometry(nc)). But with my polygons, the st_sample step always gives this error:

st_as_s2(): dropping Z and/or M coordinate Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, :
Loop 0 is not valid: Edge 15 crosses edge 19

I have tried inserting sf_use_s2(FALSE), but I still get the same error. The plot command does work, and the polygons look as expected. My code:

core<-st_read("core.kml")
core<-st_geometry(core) 
sf_use_s2(FALSE) 
core.points<-st_sample(core,5) 
plot(core)

• r
• spatial

• 1 @Friede , it's st_make_valid() that generates a multipolygon when it tries to fix self-intersections. @chill, did you draw or generate those shapes yourself? If yes, you might have had some issues with shape editor or with correct point order - i.sstatic/08H7SuCY.png – margusl Commented Nov 19, 2024 at 12:10
• 1 Yes, that was meant like a warning. @chill, can you share where core.kml comes from? – Friede Commented Nov 19, 2024 at 12:22
• 1 @Friede I'm checking with my former student who created the files. I'm pretty sure she drew them in Google Earth Pro, saved them as kmz's, and then I converted to kml. But I shot her off an email and will update when I hear back. Neither of us came into this project with any experience with spatial analyses - very much learning as we go. Your solution was great! But yes, resolving what was wrong in the first place would be good, too. – chill Commented Nov 20, 2024 at 14:03
• 1 @Friede, yes, she made them in Google Earth Pro. I will zoom in on them and see if I can spot a problem at that level, by inspecting the edges. – chill Commented Nov 21, 2024 at 0:00
• 2 @margusl YES it turns out the original polygons had some folded, overlapping edges. Once I edited them then the simpler code from the answer I linked in my original question worked. Thank you both for guiding me to the place to look for the source of the errors I was getting. – chill Commented Nov 21, 2024 at 2:36

1 Answer 1

If it's about to ensure that there are points on each polygon, you might like to try something along the lines

 library(sf)
 sf_use_s2(FALSE) 
#> Spherical geometry (s2) switched off
 core = st_read("/Users/cara/Desktop/core.kml") |> st_geometry() |> st_make_valid() 

 core 
#> Geometry set for 4 features 
#> Geometry type: GEOMETRY
#> Dimension:     XY
#> Bounding box:  xmin: -79.06133 ymin: 33.69862 xmax: -78.90788 ymax: 33.8351
#> Geodetic CRS:  WGS 84
#> MULTIPOLYGON (((-79.014 33.80998, -79.00692 33....
#> POLYGON ((-79.06133 33.82931, -79.05514 33.8306...
#> POLYGON ((-78.99585 33.78962, -78.99496 33.7901...
#> POLYGON ((-78.91353 33.70372, -78.91197 33.7052...
 
 ppwa = \(sfc, pts) { 
   # points per weighted area 
   a = st_area(sfc)
   p = proportions(a)
   as.numeric(p * pts) 
 }
 
 # ensure sum(n) == pts, see below
 floor_preserve = \(x) {
   fx = floor(x)
   i = tail(order(x - fx), round(sum(x)) - sum(fx))
   fx[i] = fx[i] + 1
   fx
 }
 
 n = ppwa(core, 25) |> floor_preserve()
 n 
#> [1] 14  2  6  3
 sum(n)
#> [1] 25
 
 set.seed(1119)
 p = st_sample(core, n)
#> although coordinates are longitude/latitude, st_intersects assumes that they
#> are planar
#> ... gets repeated ...

 plot(core) 
 points(p, pch = 4, col = palette.colors()[7])
#> Warning in st_centroid.sfc(x, of_largest_polygon = of_largest_polygon):
#> st_centroid does not give correct centroids for longitude/latitude data

Notice that (1) we haven't corrected issues with core. Please have a look on the warnings. This might be a good discussion to start your research; islands are difficult; (2) the total number of points is not random, as we set it to 25. The number of points per polygon is therefore not random randon: it's 25 multiplied by the relative area of each polygon, i.e. weighted.