1

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

6
  • 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.
    – Vince
    Jul 22, 2022 at 13:06
  • 2
    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).
    – Spacedman
    Jul 22, 2022 at 13:54
  • 1
    If you are using sf classes you should use st_voronoi.
    – Spacedman
    Jul 22, 2022 at 13:54
  • 1
    @JeffreyEvans I've got my eye in for typos at the moment, so you mean sf::st_distance
    – Spacedman
    Jul 22, 2022 at 13:57

1 Answer 1

4

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.

3
  • how can I change the color, because two stations are in the green area?
    – user204173
    Jul 24, 2022 at 9:45
  • for a better distinguish
    – user204173
    Jul 24, 2022 at 9:45
  • 1
    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.
    – Spacedman
    Jul 24, 2022 at 10:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.