R – How to Find Which Region is Nearest to Which Stations

proximityrvoronoi-thiessen

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

Best Answer

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")
plot(stations$geom, add=TRUE, pch=22, col="black", bg="red")

enter image description here

Note "nearest" here is using a distance measure of the nearest distance from point to any point in or on the polygon feature. Some parts of features may be nearer to other stations. If you want to divide the polygon features by this criterion then you do need to apply voronoi polygons and do an intersection.

Related Question