0

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.

6
  • One way to do it might be to generate a "template" set of points that are known to be big enough to span any electoral division, and then shift that template to the centroid of each division and do a point-in-polygon test. If I understand correctly the output for each division is four lines of points going in the four cardinal directions, spaced out by 400m, from the centroid?
    – Spacedman
    Aug 5, 2021 at 13:09
  • Do you want these points to stop when the boundary is first reached, so that if you have a really irregular concave-shaped wiggly boundary you might get only a point 400m N of the centroid but not 800m (even though 800m is inside the division) because the boundary cuts across? To do that you could use the "template" approach above but create line segments and test for intersection.
    – Spacedman
    Aug 5, 2021 at 13:12
  • Ideally I would like every 400m points to the boundary so going west from the centroid to the border could be 1km therefore 2 points would be generated, while for the electoral division going north could be 500m so 1 point would be generated. Would the template work for this? In my head I have a solution but I can't appear to find a method to do it. My idea is to select long/lat point on the boundary by direction (n,e,w,s) of an electoral division. If I had those 4 points on the boundary I could calculate the distance from the centroid and have that distance as the constraint.
    – cmdublin
    Aug 6, 2021 at 14:44
  • Do you need the precision of spherical geometry? Or could you work with a projected coordinate system? I don't think the difference will be much for a 400m length of a typical area about 1/3400 the size of Ireland.
    – Spacedman
    Aug 7, 2021 at 10:27
  • I've got 3409 boundaries from census.cso.ie/censusasp/saps/boundaries/ED%20Disclaimer1.htm and some areas have centroids that are outside the boundary. How do you want to deal with those?
    – Spacedman
    Aug 7, 2021 at 10:57

1 Answer 1

0

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)

enter image description here

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)

enter image description here

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

enter image description here

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.

6
  • This is exactly it. Thank you so much. I'll apply to my dataset and let you know how I get on. I used the centroids that were provided by the OSI for each electoral division. I have found Ireland tricky for spatial analysis, particularly working with electoral divisions (see Navan as an example). The most northern point is perfectly fine as these points will be used to calculate coverage within the radius.
    – cmdublin
    Aug 7, 2021 at 17:22
  • Also in areas that are less than 400m, that's ok as the centroid covers the area in terms of radius.
    – cmdublin
    Aug 7, 2021 at 17:24
  • I've been stuck today taking the list further. It worked perfectly on my data but I'm trying to unnest the list so I can extract all of the geometry points and convert to lon/lat coordinates. I need to create a separate dataframe with these coordinates so I can apply to mclp. Can you share the best package to unnest the list? I can manually save each item minus the areas that were too small for coordinates which I've done on a few electoral divisions.
    – cmdublin
    Aug 8, 2021 at 15:26
  • do.call(rbind,list_of_things) will make a data frame from a list just as if you did rbind(d1,d2,d3, etc....). Not sure if the empty ones will make this fail though, and you might want to add a column to each data frame to show which polygon it came from if you are going to drop some from the unnested data frame.
    – Spacedman
    Aug 8, 2021 at 16:14
  • Thanks for that.Hadn't come across do.call before, great function.I deleted the empty electoral divisions from the list which returned 3292. I applied do.call using tset<-do.call(rbind, list(pts)) & it returned a matrix. I've unlisted it by df<- data.frame(matrix(unlist(tset), nrow=137662, byrow=TRUE, ncol = 2),stringsAsFactors=FALSE) returns geometry coordinates. Separately, I've tested code on an isolated feature to transform the coordinates to long/lat as it won't work on the entire pts. test <- st_transform(pts[[1]], CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")) I'm nearly close..
    – cmdublin
    Aug 8, 2021 at 20:12

Your Answer

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

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