I currently have two data sets concerning the Netherlands, one containing longitudes and latitudes of population cores (Hubs) along with other information, such as average house price, and the other containing the longitudes and latitudes of particular amenities (References).

I would like to create a Voronoi tessellation to map each area of the Netherlands to a particular population core. I then wish to create circles of varying radii for each amenity and create a dataframe of the areas that each Voronoi section is overlapped by each circle.

I have currently used a function from R: Voronoi tesselation from long/lat data? to generate a SpatialPolygonsDataFrame mapping the Voronoi tessellation. I then used Create a circle of defined radius around a point and then find the overlapping area on a shapefile to create a data.frame of polygons of radius 100m.

voronoipolygons = function(layer) {
    crds = layer@coords
    z = deldir(crds[,1], crds[,2])
    w = tile.list(z)
    polys = vector(mode='list', length=length(w))
    for (i in seq(along=polys)) {
        pcrds = cbind(w[[i]]$x, w[[i]]$y)
        pcrds = rbind(pcrds, pcrds[1,])
        polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
    SP = SpatialPolygons(polys)
    voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(dummy = 
    seq(length(SP)), row.names=sapply(slot(SP, 'polygons'), 
    function(x) slot(x, 'ID'))))

which is the function that performs the tessellation on the 'hubs' data. The remaining shape function is then created and compared:

Hubs = read.csv("Test_Hub_Data.csv", header=TRUE)
coordinates(Hubs) <- c("X.Coord", "Y.Coord")

#The function is applied to the Hubs data to create the dataframe
vp <- voronoipolygons(Hubs)

#The second data set is read in
References = read.csv("Test_Reference_Points.csv", header = TRUE)

#Coordinates are assigned to convert to a sf
coordinates(References) <- c("X.Coord", "Y.Coord")

#The X and Y names actually correspond to longitude and latitude of each point 
dat_hold <- st_as_sf(References, coords = c("X.Coord", "Y.Coord"), crs = 4326) %>% st_transform(3035)

#Circles are created at each of the points in References at 100m
dat_circles <- st_buffer(dat_hold, dist = 100)

#The Voronoi data is changed to a shape function
VoronoiMap <- st_as_sf(vp)

#The intersections of both shape functions are determined
dat_Intersec <- st_intersection(VoronoiMap, dat_circles) 

However, when running this function an error occurs stating that st_crs(x) == st_crs(y). I assume that this is a formatting problem with the inputs, but I am unsure how to fix it.

I've stated the entire problem, as I'm sure that there is an easier way of completing this using R functions, with the calculation of overlapping areas being a common method comparing a Voronoi Tessellation that I'm sure alternatives exist.

A sample test data set that I've been using is included here:

X Coord Y Coord Info 1  Info 2  Info 3
52.37639999 4.894589901 1   9   13
52.36989975 4.876679897 2   8   15
52.86299896 6.493730068 3   7   62
52.01011658 4.703330517 4   6   24
51.67409897 4.548630238 5   5   9
50.84980011 5.687580109 6   4   46
51.88669968 4.491449833 7   3   73
52.1048943  5.0528496   8   2   97
52.0746994  4.310130119 9   1   69

and the other data set corresponding to amenities:

X Coord Y Coord
51.88529968 4.493549824
52.3360136  4.8750552
52.37440109 4.895800114
51.92269897 4.47382021
51.9192276  4.481627464
51.43019867 5.451980114
52.0780896  5.1470653
51.90299988 4.458670139
52.04100037 4.318089962
50.98810196 5.78059721
  • Where in your code do you specify the coordinate system for the "Hubs"? – whuber Jan 29 '18 at 15:36
  • @whuber I tried to use the coordinates function just after initialising the sp library. If you're referring to the crs = 4326 %>% st_transform(3035) for the sf conversion, I tried using Vor_hold <- st_as_sf(vp, coords = coordinates(vp), crs = 4326) %>% st_transform(3035). Have I missed something obvious? – AlwaysInTheDark Jan 29 '18 at 15:58
  • I'm referring to the absence of any corresponding code to specify the coordinates of vp, unless somehow I'm overlooking it. That's what the error message is telling you, too: it says you're trying to compare two datasets with different coordinate systems. – whuber Jan 29 '18 at 16:02
  • What type of coordinates do you mean? I'm quite new to this, so when I run coordinates(vp) the coordinates of Hubs are returned in form c("X.Coord","Y.Coord"), which I thought meant that it was in a reasonable format. I thought the issue might be that I'm classifying the coordinates using the crs = 4326 and 3035 systems, but I don't know how to fix that. Was this the system you were referring to? – AlwaysInTheDark Jan 29 '18 at 16:10
  • The coordinates are the numbers. They are meaningless for most geographical analysis until you explain (to the software) what they actually refer to: lat-lon, meters, feet? What coordinate system? The software is expecting you to provide consistent descriptions; otherwise it figures you might be making a mistake and it complains. The solution is to tell it explicitly what coordinate system is associated with each dataset. – whuber Jan 29 '18 at 19:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.