How to calculate near neighbours in sf object with polygon geometry in R?

I have polygon shape data for German postcodes. For those postcode polygon I like to calculate various nearest neighbour measures. I have seen that procedures working with the sp package (using coordinates(), like knearneigh(coordinates(GER), k = 4)). I opt for sf spatial objects in R and are confused on how to implement neighbours here. Thank you

library(sf)
library(dplyr)
library(leaflet)

URL <- "https://downloads.suche-postleitzahl.org/v2/public/plz-5stellig.shp.zip"

# use GDAL virtual file systems to load zipped shapefile from remote url
GER_postcode <- paste0("/vsizip//vsicurl/", URL) %>%  read_sf()

# country outline from giscoR
GER_outline <- giscoR::gisco_get_countries(country = "DE")

# subsample 
GER_postcode_subsample <- GER_postcode %>% filter(substr(plz, 1, 1) %in% c(0, 1, 7))

# k nearest neighbours for sf dataframe

Solution

I found the answer in the spdep package which contains the speaking function poly2nb(). Don't know why I havn't found this earlier.

library(spdep)
queens <- poly2nb(GER_postcode_subsample, 
                     queen = TRUE, # a single shared boundary point meets the contiguity condition
                     snap = 1) # we consider points in 1m distance as 'touching'
summary(queens)

Neighbour list object:
Number of regions: 2624 
Number of nonzero links: 13698 
Percentage nonzero weights: 0.1989434 
Average number of links: 5.220274 
9 regions with no links:
275 284 554 616 922 947 1889 2328 2329
Link number distribution:

  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14 
  9  28 117 307 513 574 468 310 168  78  28  11   2  10   1 
28 least connected regions:
112 297 298 474 529 809 843 852 896 917 921 946 951 1027 1050 1147 1524 1687 1884 2068 2271 2291 2314 2327 2343 2367 2368 2509 with 1 link
1 most connected region:
1349 with 14 links