R – How to Calculate Distance from Centroid to Border of Spatial Polygon in R

Question

I have 3409 Electoral Divisions from the Irish census for the Republic of Ireland. I want to calculate 400 metre distances from the centroid of every electoral division to the spatial polygon boundary in the direction of north, south, east and west. I've figured out how to calculate 400 metres in each direction with the output as long/lat coordinates, but not how to iteratively do this until the boundary is reached.

Here's an example of the code calculating 400 metres from the centroid going north:

library(geosphere) 
pn <- cbind(DF$long, DF$lat) 
bn <- 360 
dn <- 400 
an <- 6378137 
fn <- 1/298.257223563 
DF$North <- destPoint(pn, bn, dn, an, fn)

How can I do this using R to code for iteratively calculating long/lat coordinates of 400 metres from the centroid until the boundary of each 3409 electoral divisions is reached? Please note the columns that I have in the spatial polygon dataframe is centroid_x, centroid_y, long_centroid, lat_centroid, shape_area, shape_length.

Any tips or pointers for code would be most appreciative.

Answer 1

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.