How to combine sf_intersects and sf_touches within a function

Viewed 42

Hi All,

I'm trying to develop a function which maps point coordinates in the UK to the UK local authority areas in which they reside in.

I've developed a function using the st_intersects function of the sf package and it works to identify cases where the centroid intersects with the shapefile of a given area.

However, I also want to identify cases where the centroid actually touches the LA area boundary line.

I'm aware that where LA areas share a boundary this would map the centroid to two LAs but that's ok for my purposes.

I tried to adapt my function to include st_touches however it produces an empty table.

Would be very grateful if someone could take a look and suggest where I might be going wrong.

Reprex below:

The LA area shapefile is publicly available here (the shapefile from the downloads section):

https://geoportal.statistics.gov.uk/datasets/97a614bdcc6043bd9a3cbfbba8a1f302_0/explore?location=55.230176%2C-3.316939%2C6.92

The set of long/lat coordinates to map:


data = structure(list(longitude = c(-2.87481491108804, -2.87127368026435, 
-2.96479539064797, -2.92511948140495, -3.01836788419, -2.98402990798746, 
-2.9077264274875, -2.89347586558879, -2.97850780890271, -3.02924916756829, 
-3.05411812798672, -2.93694678850526, -2.88306663657761, -2.9845903484732, 
-3.05984423623939, -2.90909675996297, -2.92181750227365, -3.02464633153093, 
-3.02123789197677, -3.09890063032602, -2.96844903903745, -2.86785459942595, 
-2.9321154435321, -3.10888230637966, -2.97222428832271, -2.89647987410202, 
-3.03027882189321, -2.97866205698605, -2.99429382433499, -3.05246646242959, 
-2.89616224682703, -2.98563264041128, -2.95302953391751, -2.91151041579422, 
-2.97953110815078, -3.10796669571098, -2.94570764658262, -2.97301831593517, 
-2.89323494964768, -2.97063026715387, -2.92198857953645, -3.03675565980253, 
-2.83891728241267, -2.91861697258726, -3.06892905484705, -2.99398588944906, 
-2.93551831989608, -2.82716052265689, -3.01724727594163, -2.92682901324458, 
-2.91803014985187, -2.99284666045525, -2.97786793943296, -2.9230445053825, 
-2.95481455034888, -2.93676036122089, -2.96585002429569, -2.86998042564442, 
-2.93970334000959, -2.98376964052552, -2.87448694998545, -2.95369805041689, 
-2.877566853416, -2.99224636258483, -3.09026039652281, -2.95817525261748, 
-2.95017575249327, -2.94359797634677, -2.96710539009414, -3.07792734952214, 
-3.17993626061202, -2.97741036119629, -2.95222841522805, -2.95965218909241, 
-2.9459062849359, -2.86310994503288, -3.03906859820527, -3.02177819639329, 
-2.97166265836844, -3.08145848394861, -3.03808462354296, -2.96263651468074, 
-2.95851898803173, -3.02221082538059, -2.89362064147485, -2.8993882625834, 
-2.92770961205626, -2.96369118745873, -2.89571920584541, -3.09314291426038, 
-3.0919937117765, -3.03611523902419, -2.93117415659881, -2.97460696743252, 
-2.85036432827623, -2.94987226425819, -2.95390181653767, -2.91019481240696, 
-2.99966589456369, -2.93169619777761, -3.08411731249804, -3.1706199707398, 
-2.91896769518433, -2.97090655977713, -3.11395127242533, -3.04781510226649, 
-3.18378788354865, -2.90075209953503, -2.93948388571496, -2.98732033573349, 
-2.96611890948035, -2.87803392597212, -3.02717324311628, -2.90642650501013, 
-3.0181126772633, -2.98716889544258, -2.95609545040814, -2.89457583588837, 
-2.9183348068514, -2.88807533717785, -3.02947911076491, -2.98582946684398, 
-2.90956913387858, -2.92612896591959, -2.8618485532865, -2.93181934586936, 
-2.99086525128266, -2.88799987536917, -3.18188458893135, -2.92530589829918, 
-2.84517946221428, -3.02066308251878, -2.943131290944, -2.9316549262618, 
-2.89416059133544, -3.03093156077094, -3.03954865354919, -3.02892322932484, 
-3.15157580297919, -2.97019350603118, -2.95926518021845, -3.12995261693527, 
-3.0178786559463, -2.96998299026556, -2.92983501387733, -2.91526010044323, 
-2.98893750456704, -2.93303968342649, -2.97826195863277, -3.10231870679243, 
-2.86918162475871, -2.94860479427078, -2.89593919501248, -3.04120959095772, 
-3.08974718532665, -3.05135478822517, -2.98129768711627, -3.03665798566578, 
-2.89429626072004, -2.92008555617904, -2.87339557195832, -2.89374834459854, 
-2.97879863154011, -2.95902409390154, -2.8950914531408, -2.96902081781808, 
-2.95976145407803, -2.92255923707144, -3.05150511302466, -3.05903671872711, 
-2.95262009289705, -3.00863929575895, -3.0547892173777, -3.03276242335601, 
-2.84477534615073, -3.11863250997746, -3.03915412086859, -2.93123288923061, 
-3.11376765617702, -3.15438194594307, -3.02283334806073, -2.95561312619315, 
-3.02508648595712, -3.15483442747855, -2.96632534070746, -2.98706717919529, 
-2.86455911629215, -2.97919924773756, -3.14477357952329, -2.91786092204208, 
-3.048045603786, -2.92396930828526, -2.87962098053312, -2.93684151727313, 
-3.11262802664107, -2.90201460498108, -2.94757682014421, -2.97896151794228, 
-2.98981135039931, -2.95833097373606), latitude = c(53.4379570197908, 
53.4201221364593, 53.3942986189115, 53.3842844842061, 53.3617791885161, 
53.4151651024398, 53.4508620255715, 53.3939173214826, 53.4064101823484, 
53.3718255079932, 53.4378662280985, 53.4375126675335, 53.4206467370139, 
53.4029535199428, 53.4259950343029, 53.4140285695349, 53.3888028546485, 
53.3879553701277, 53.383130581584, 53.3980816902538, 53.4094946698524, 
53.4280841646817, 53.4251728245991, 53.4065738999879, 53.4338598436095, 
53.4259310929223, 53.3889496585032, 53.3412480513966, 53.4088327587453, 
53.4282695483082, 53.3714698780691, 53.3992414758791, 53.4283426004897, 
53.3694659558339, 53.4034617561608, 53.3852333876622, 53.4307199421458, 
53.4145901432611, 53.3443452440954, 53.3896653425554, 53.408541145569, 
53.4215685101068, 53.3869120986626, 53.3861128798233, 53.4347861533497, 
53.4001430176595, 53.3902962809302, 53.3444171680078, 53.3626786684076, 
53.3737795084725, 53.4356820437825, 53.4061211555453, 53.4062807459519, 
53.4440809390267, 53.4300901595633, 53.4454602665698, 53.4136326339338, 
53.3744591704793, 53.4266774115989, 53.4103300396462, 53.4245895519719, 
53.4252630671781, 53.4012893604736, 53.4258121557628, 53.3283517163451, 
53.4296137958026, 53.3894613832297, 53.3911944946134, 53.4311598021066, 
53.4019404403309, 53.3660259655574, 53.4066620280597, 53.4052301130645, 
53.4610898898285, 53.3834008004254, 53.396693596265, 53.3732415469181, 
53.3830899845843, 53.4145203148397, 53.384011307051, 53.385933567503, 
53.4418746881753, 53.3893676837379, 53.3754186100903, 53.3477042237358, 
53.4268260826889, 53.4529736611373, 53.4009301384742, 53.4255412928395, 
53.3799946476186, 53.3832592444011, 53.4280280487181, 53.4371892773805, 
53.3873993179818, 53.387265880647, 53.4311902629893, 53.4000835375925, 
53.416299149301, 53.3472942931141, 53.4335267442322, 53.3780004976281, 
53.3901906552793, 53.3856427565351, 53.4064274548061, 53.3915868687672, 
53.3988915899388, 53.38243801477, 53.3582036174684, 53.4310567146736, 
53.4071737491274, 53.4613164150686, 53.4509761373123, 53.3767513171219, 
53.4386111541956, 53.3813234955064, 53.4115730142442, 53.4690375638476, 
53.3925877169023, 53.4149064630925, 53.4378681311582, 53.363283937878, 
53.3439845549594, 53.3919170930087, 53.4003921305836, 53.3729720867359, 
53.3971837533247, 53.4043577837836, 53.4356933993701, 53.3940528400041, 
53.4497743100287, 53.3896667782373, 53.4087450666446, 53.4421561637244, 
53.4419585659703, 53.3476642223286, 53.3773841273259, 53.4318096029357, 
53.359549289032, 53.3954996343164, 53.4603395616878, 53.3790873492044, 
53.3995806276861, 53.3923100412117, 53.4610154343059, 53.4344581137143, 
53.4181570758271, 53.4185223637193, 53.3972820875749, 53.402475229261, 
53.4081268716824, 53.4237868466879, 53.3959009588377, 53.4045147184146, 
53.4117047911735, 53.3747671044637, 53.4265465099582, 53.3173604267677, 
53.3853886811948, 53.4071671211117, 53.4458146150033, 53.4509562137809, 
53.3545258760886, 53.4081608055445, 53.4639623704057, 53.4066875108909, 
53.4326544056812, 53.3958027859543, 53.378011861023, 53.4143471677853, 
53.3916913284024, 53.3831046999105, 53.3526121691896, 53.424242035722, 
53.4170165452222, 53.3426757282892, 53.3946264472929, 53.4299953680367, 
53.3730800431989, 53.3787971976028, 53.3974951074273, 53.3894123258524, 
53.4079890811827, 53.3898479198151, 53.4014279834277, 53.4445415555092, 
53.4092342886507, 53.3858335487274, 53.4380366162007, 53.4083295907901, 
53.4145595413981, 53.4247867968646, 53.4332901497374, 53.4017057070297, 
53.4238324512685, 53.3875542058095, 53.3639110429771, 53.4288174617707, 
53.3844646095652, 53.4258501205813, 53.4004076184837)), row.names = c(NA, 
-200L), class = c("tbl_df", "tbl", "data.frame"))

Loading the shapefile

library(dplyr)
library(data.table)
library(sf)
library(purrr)

shapefile <- sf::read_sf("downloads/LAD_MAY_2021_UK_BFC.shp")

The function with st_intersects only (works):

getBuildings <- function(data, shapefile) {
  
  shapefile <- sf::st_transform(shapefile, crs = 4326)
  
  pnts_sf <- sf::st_as_sf(data, coords = c('longitude', 'latitude'), crs = sf::st_crs(shapefile))
  
  locations_by_area<- pnts_sf %>%
    dplyr::mutate(
      intersection = as.integer(sf::st_intersects(geometry, shapefile$geometry))
      , LAD21NM = dplyr::if_else(is.na(intersection), '', shapefile$LAD21NM[intersection])
    ) %>%
    data.table::setDT() %>%
    dplyr::mutate(latitude = unlist(purrr::map(geometry,2)),
                  longitude = unlist(purrr::map(geometry,1)))
  
  locations_by_area
  
  
}

results = getBuildings(data = data, shapefile = shapefile)

The function with st_touches (doesn't provide any results):

getBuildings <- function(data, shapefile) {
  
  shapefile <- sf::st_transform(shapefile, crs = 4326)
  
  pnts_sf <- sf::st_as_sf(data, coords = c('longitude', 'latitude'), crs = sf::st_crs(shapefile))
  
  locations_by_area<- pnts_sf %>%
    dplyr::mutate(
      intersection = as.integer(sf::st_intersects(geometry, shapefile$geometry))
      touch = as.integer(sf::st_touches(geometry, shapefile$geometry))
      , LAD21NM = dplyr::if_else(is.na(intersection), '', shapefile$LAD21NM[intersection])
      , LAD21NM2 = dplyr::if_else(is.na(touch), '', shapefile$LAD21NM[touch])
    ) %>%
    data.table::setDT() %>%
    dplyr::mutate(latitude = unlist(purrr::map(geometry,2)),
                  longitude = unlist(purrr::map(geometry,1)))
  
  locations_by_area
  
  
}

results = getBuildings(data = data, shapefile = shapefile)


The end result I'm looking for is a table of the original coordinates with another column indicating which LA area they reside in and another column indicating which LA the centroid touches the border of (I'm guessing there would be two rows for some centroids where they touch a border shared by two LA areas).

Any help would be greatly appreciated, please and thank you.

0 Answers
Related