Here's what I have. First a function to generate the points going out in the cardinal directions:
cardinal <- function(n, sep){
### generate n points separated by sep in the NSEW directions from (0,0)
delta = seq(sep, length.out=n, by=sep)
N = cbind(0, delta)
S = cbind(0, -delta)
E = cbind(delta, 0)
W = cbind(-delta, 0)
rbind(N,S,E,W)
}
Then this function which shifts a set of points to be centred at a given point. We'll use this to move a set of points generated by cardinal
to each centroid:
addoff <- function(pt, offs){
### add the offsets in `offs` to the point in `pt`:
offs[,1] = offs[,1]+pt[,1]
offs[,2] = offs[,2]+pt[,2]
offs
}
Then the main function that takes an `sf` spatial data frame and a pattern of cardinal points:
nsew <- function(geoms, cps, messages=FALSE){
allpts = lapply(1:nrow(geoms), function(i){
if(messages)message("region ",i)
geom = st_geometry(geoms)[i]
pts = addoff(st_coordinates(st_centroid(geom)), cps)
pts = st_as_sf(data.frame(pts), coords=1:2, crs=st_crs(geoms))
inside = which(lengths(st_intersects(pts, geom))==1)
if(length(inside)==0){
pts = numeric(0)
}else{
pts = pts[inside,]
}
pts
})
allpts
}
To use, you first have to create a set of points in the cardinal directions to span the largest of your polygons. If its too small you'll miss points, but it will be slower if you have too many. This function should compute the number of steps of size w
needed to span the biggest bounding box in the set, which should be sufficient in the worst case scenario that one of your regions is a rectangle:
sizeneeded <- function(regions, w){
boxes = data.frame(t(sapply(st_geometry(regions), st_bbox)))
widths = boxes$xmax - boxes$xmin
heights = boxes$ymax - boxes$ymin
maxbox = max(c(widths, heights))
n = round(1 + (maxbox/w))
n
}
Okay that's the setup. Now...
# read the electoral divisions
ed = st_read("./electoral_divisions.shp")
# how many points max do we need?
ns = sizeneeded(ed,400)
print(ns) # its 61
# create the points
cps = cardinal(ns,400)
# loop over all EDs and find the cardinal points in `cps` that
# centred on the centroid fall in the ED. This may take a few minutes:
pts = nsew(ed, cps)
That is a list of the same length as the number of features in ed
, so to plot an ED and then show the points:
> plot(ed$geometry[543])
> plot(pts[[543]], add=TRUE)

Note that some EDs have no points in them. Region 80 is very small:
> plot(ed$geometry[80])
> axis(1)
> axis(2)
> pts[[80]]
numeric(0)

Note how with this method points don't stop when the first NSEW direction meets the border (region 1209):

The northern points skip one point where the edge jinks in, and the points continue. Not sure if this is what you want or not. If not, then some rewriting is needed.