# Finding which region is nearest to which stations in R

I have a point dataset of 16 places in Switzerland that looks like these and is called stations

``````stations <- read_csv("station.csv)
view(stations)
``````
Cityid x y
1 611285.549 267664.109
2 600454.843 199648.553
3 754069.768 225837.065
4 716492.326 112276.605
5 782974.270 186304.740
6 500013.259 117819.944
7 554010.979 217120.979
8 538122.750 152357.507
9 704699.521 113740.950
10 717157.313 95810.122
11 666156.604 211382.117
12 734356.986 277041.847
13 594091.196 119770.519
14 683125.786 247005.951
15 81150.847 720009.498
16 690036.302 283496.466

I have also a shapefile municipality in Switzerland and their geometry

``````ms_shp <- st_read("CH.shp")
``````

Now I would like to create Thiessen polygons to find out which regions belongs to which station according to the shortest distance and to have a list which stations are the nearest to which stations.

I used that code but got an error

``````library('dismo')
points <- matrix(c(stations\$x,stations\$y), ncol=2)
voronois <- voronoi(points)
spplot(voronoi)
``````

The CRS of the datapoints and the shapefile is equal.

The error message is

``````argument not used
``````

Is there another method, that I can find which regions are nearest to which stations?

I thought Thiessen polygon is a good method, but maybe another method in R is more useful

• Your goal and Thiessen are incompatible. The best you'll get is polygons overlapping polygons. Calculating the distance of each point to each polygon and choosing the closest is the solution. Doing that in R is left to the arcane-tolerant. Jul 22, 2022 at 13:06
• You could just use `sf::st_disiance` Jul 22, 2022 at 13:07
• You've done `voronois <- voronoi(points)` and `spplot(voronoi)` not `spplot(voronois)` (with an extra `s`). Jul 22, 2022 at 13:54
• If you are using `sf` classes you should use `st_voronoi`. Jul 22, 2022 at 13:54
• @JeffreyEvans I've got my eye in for typos at the moment, so you mean `sf::st_distance` Jul 22, 2022 at 13:57

Make `stations` into an sf points data frame with same crs as polygons:

``````stations = st_as_sf(stations, coords=c("x","y"), crs=st_crs(ms_shp))
``````

create new column in the polygons of the index of the nearest station feature:

``````ms_shp\$nearest_station = st_nearest_feature(ms_shp, stations)
``````

plot:

``````plot(ms_shp\$geom, col=ms_shp\$nearest_station, border="white")
• Use a different palette (see `help(palette)`) or use `tmap` for maps. You really need a palette with 16 colours that are distinguishable by human eye, see the `pals` package for lots of possibles. Ask a new question if you get stuck. Jul 24, 2022 at 10:29