nearest polygon

# msoa: 7,264 Middle-Layer Super Output Census small areas in England and Wales (polygons)
# points: 1,531,286 active postcodes in England and Wales

library(data.table)
library(sf)

msoa <- RcensusUK::MSOA |> st_transform(27700)
points <- RpostcodesUK::postcodes[is_active == 1 & CTRY %in% c('ENG', 'WLS'), .(PCU, x_lon, y_lat)] gpoints <- st_as_sf(points, coords = c('x_lon', 'y_lat'), crs = 4326) |> st_transform(27700)
points[, MSOA := st_join(gpoints, msoa, join = st_within) |> subset(select = MSOA) |> st_drop_geometry()]
neighs <- spdep::poly2nb(msoa)
output <- data.table(PCU = character(0), nearest = character(0))
for(x in 1:nrow(msoa)){
    message('Processing ', x)
    xm <- msoa$MSOA[x]
    n_msoa <- msoa[neighs[[x]],]
    p_in_msoa <- points[MSOA == xm, .(PCU)]
    gp_in_msoa <- gpoints |> subset(PCU %in% p_in_msoa$PCU)
    lkps <- st_distance(gp_in_msoa, n_msoa) |> t() |> as.data.table()
    p_in_msoa[, nearest := rbindlist(
        lapply(
            names(lkps), 
            \(z) n_msoa[lkps[, .I[order(get(z))]][2],] |> st_drop_geometry()
        )
    )]
    output <- rbindlist(list(output, p_in_msoa))    
}

small example

library(sf)
library(data.table)

msoa <- RcensusUK::MSOA |> st_transform(27700)
points <- data.table(
            'x_lon' = c(-0.696258, -2.280421, -3.747959), 
            'y_lat' = c(52.750425, 52.100153, 50.771948)
)
gpoints <- st_as_sf(points, coords = c(1, 2), crs = 4326) |> 
                st_transform(27700)
points[, MSOA := st_join(gpoints, msoa) |> st_drop_geometry()]
lkps <- st_distance(gpoints, msoa) |> t() |> as.data.table()
points[, nearest := rbindlist(
                        lapply(
                            names(lkps), 
                            \(x) msoa[lkps[, .I[order(get(x))]][2],] |> st_drop_geometry()
                        )
)]
points

library(leaflet)
leaflet() |> 
    addTiles() |> 
    addMarkers(data = gpoints |> st_transform(4326)) |> 
    addPolygons(data = msoa |> subset(MSOA %in% points$MSOA) |> st_transform(4326)) |> 
    addPolygons(data = msoa |> subset(MSOA %in% points$nearest) |> st_transform(4326), color = 'red')