Attaching data to SpatialPolygonsDataFrame

Question

As this is a more GIS related question (and maybe solvable with ogr/gdal functions I'm not aware of), I try it here first, although working with R.

What I want to achieve: I have some data, based on that data I want to create Voronoi/Thiessen polygons, but within specific geographical margins (so basically just an intersection issue). This works and is not a problem.

However, when intersecting the SpatialPolygon and change it to a SpatialPolygonsDataFrame, I have only the geometry, but no data attached to it. When trying to attach the data, I get an error:

Object length mismatch:foo has 107 Polygons objects, but data has 92 rows

which is strange, as I used this data to create the polygons. The only thing I can come up with is that due to the intersection some polygons at the margins have been splitted in multiple parts (yet just a guess).

So here is some sample code to reproduce the problem (you can download the shp here). Most of the code is working, the issues arise only in the last comment section, so you may only pay attention to the last couple of lines:

# get data and projection
library(rgdal)
X <- c(579443.1,507821.9,527521.1,580656.9,462417.3,598714.2,528403.5,590470.2,621899.1,552101.1,579211.2,524459.0,622132.6,568279.1,551119.6,549101.4,626778.5,570531.6,498509.8,510324.2,501678.7,497111.0,566879.1,628653.7,558457.7,504368.8,482406.3,493979.3,591583.4,580834.5,579033.3,482417.6,594981.1,640817.4,583138.0,552565.4,480051.3,538693.4,595036.1,601637.4,587406.8,513506.1,511508.2,528562.7,508318.7,505132.4,560761.3,533665.2,569308.3,524761.9,538409.6,573477.6,618846.6,608033.6,601203.8,552204.2,515622.8,582403.8,493001.6,481086.7,563902.5,487437.6,626197.1,525136.5,510164.3,551818.6,486836.1,507096.5,556712.9,639675.3,553121.2,551446.2,568682.8,576450.7,490634.5,545741.7,503393.7,591073.8,561552.8,545749.3,524766.0,566276.8,578132.3,535480.5,610635.8,513964.5,520131.4,534375.5,568933.5,509446.1,475644.0,516545.3)
Y <- c(5324384,5413164,5302565,5268202,5320375,5308655,5399552,5279586,5400961,5315579,5347046,5314407,5338724,5314637,5364789,5405192,5387508,5279402,5401240,5370434,5315436,5380213,5303072,5356320,5320272,5426778,5309010,5285152,5256444,5300410,5380396,5292246,5271598,5326458,5388330,5310118,5382563,5296763,5351272,5397560,5334196,5300293,5301863,5372164,5394113,5300187,5379644,5352554,5337499,5340189,5396554,5298737,5381064,5319755,5372751,5290355,5293599,5286506,5337546,5330379,5338980,5304273,5334241,5381226,5356570,5333053,5402336,5374597,5411038,5348579,5294778,5374055,5325943,5270214,5316524,5335903,5324051,5359391,5327128,5380892,5320120,5273887,5356543,5308660,5361716,5400485,5405815,5417928,5264594,5330510,5324040,5383757)
Z <- c(339,590,924,681,495,233,568,419,214,1034,303,468,152,488,196,464,176,509,498,787,398,680,478,151,429,571,416,1492,736,811,236,495,685,145,237,524,919,751,232,183,318,1430,807,666,575,611,333,495,276,284,466,480,192,198,280,619,1076,361,486,401,244,402,149,554,904,298,539,586,451,145,601,217,855,870,634,285,338,170,583,429,469,1100,177,480,183,526,498,460,1743,263,343,595)
data <- data.frame(X,Y,Z)

data.xy <- data[c("X","Y")]
coordinates(data) <- data.xy
proj4string(data)=CRS("+init=epsg:32633")
borderorig <- readOGR('border.shp', 'border')
border <- spTransform(borderorig, CRS("+init=epsg:32633"))

# create sp voronoi/thiessen diagram 
voronoipolygons <- function(x, poly) {
  require(deldir)
  if (.hasSlot(x, 'coords')) {
    crds <- x@coords  
  } else crds <- x
  bb = bbox(poly)
  rw = as.numeric(t(bbox(poly)))
  z <- deldir(crds[,1], crds[,2],rw=rw)
  w <- tile.list(z)
  polys <- vector(mode='list', length=length(w))
  require(sp)
  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)

  SpatialPolygonsDataFrame(
    SP, data.frame(x=crds[,1], y=crds[,2], 
                   row.names=sapply(slot(SP, 'polygons'), 
                                    function(x) slot(x, 'ID'))))  
}

v <- voronoipolygons(data, border)
proj4string(v) <- proj4string(border)
plot(v)

# intersect and export data
###### PROBLEM STARTS HERE #######
library(rgeos)
thiessen <- gIntersection(border, v, byid=TRUE) # may take a few seconds
thiessenspdf <- as(thiessen, "SpatialPolygonsDataFrame") # change SpatialPolygon to SpatialPolygonDataFrame
notworking <- SpatialPolygonsDataFrame(thiessen, data@data) # not working attempt to change to spdf and attach data to polygon feature

plot(thiessenspdf)
points(data, pch=20)

writeOGR(thiessenspdf, "B:/testR", "thiessenspdf", driver="ESRI Shapefile", overwrite_layer=TRUE)

So what I have problems with is the notworking <- SpatialPolygonsDataFrame(thiessen, data@data) part (which is kinda self explenatory).

However, I imported the shp into QGIS and used the Points in Polygon function with my data points which is exactly what I want to achieve. Although there seem to be a mismatch between number of polygons and number of points (which I can't see at all on the map), it worked perfectly and I had my point data added as an attribute to each polygon - and this is how it should look like:

I need to do this in R, however, as I need to automatise this for a huge amount of different data points.

Answer 1

Your example is not reproducible (it depends on 'border.shp').

To make a voronoi diagram:

 library(dismo)
 v <- voronoi(xy) 

To intersect while keeping attributes

 library(raster)
 x <- intersect(v, border)

After downloading the shp file:

library(rgdal)
library(dismo)

X <- c(579443.1,507821.9,527521.1,580656.9,462417.3,598714.2,528403.5,590470.2,621899.1,552101.1,579211.2,524459.0,622132.6,568279.1,551119.6,549101.4,626778.5,570531.6,498509.8,510324.2,501678.7,497111.0,566879.1,628653.7,558457.7,504368.8,482406.3,493979.3,591583.4,580834.5,579033.3,482417.6,594981.1,640817.4,583138.0,552565.4,480051.3,538693.4,595036.1,601637.4,587406.8,513506.1,511508.2,528562.7,508318.7,505132.4,560761.3,533665.2,569308.3,524761.9,538409.6,573477.6,618846.6,608033.6,601203.8,552204.2,515622.8,582403.8,493001.6,481086.7,563902.5,487437.6,626197.1,525136.5,510164.3,551818.6,486836.1,507096.5,556712.9,639675.3,553121.2,551446.2,568682.8,576450.7,490634.5,545741.7,503393.7,591073.8,561552.8,545749.3,524766.0,566276.8,578132.3,535480.5,610635.8,513964.5,520131.4,534375.5,568933.5,509446.1,475644.0,516545.3)
Y <- c(5324384,5413164,5302565,5268202,5320375,5308655,5399552,5279586,5400961,5315579,5347046,5314407,5338724,5314637,5364789,5405192,5387508,5279402,5401240,5370434,5315436,5380213,5303072,5356320,5320272,5426778,5309010,5285152,5256444,5300410,5380396,5292246,5271598,5326458,5388330,5310118,5382563,5296763,5351272,5397560,5334196,5300293,5301863,5372164,5394113,5300187,5379644,5352554,5337499,5340189,5396554,5298737,5381064,5319755,5372751,5290355,5293599,5286506,5337546,5330379,5338980,5304273,5334241,5381226,5356570,5333053,5402336,5374597,5411038,5348579,5294778,5374055,5325943,5270214,5316524,5335903,5324051,5359391,5327128,5380892,5320120,5273887,5356543,5308660,5361716,5400485,5405815,5417928,5264594,5330510,5324040,5383757)
Z <- c(339,590,924,681,495,233,568,419,214,1034,303,468,152,488,196,464,176,509,498,787,398,680,478,151,429,571,416,1492,736,811,236,495,685,145,237,524,919,751,232,183,318,1430,807,666,575,611,333,495,276,284,466,480,192,198,280,619,1076,361,486,401,244,402,149,554,904,298,539,586,451,145,601,217,855,870,634,285,338,170,583,429,469,1100,177,480,183,526,498,460,1743,263,343,595)
points <- data.frame(X,Y,Z)

coordinates(points) <- ~ X + Y
proj4string(points) <- CRS("+init=epsg:32633")
borderorig <- shapefile('border.shp')
border <- spTransform(borderorig, CRS("+init=epsg:32633"))

v <- voronoi(points)

thiessen <- intersect(border, v)

, byid=TRUE) # may take a few seconds
thiessenspdf <- as(thiessen, "SpatialPolygonsDataFrame") # change SpatialPolygon to SpatialPolygonDataFrame
notworking <- SpatialPolygonsDataFrame(thiessen, data@data) # not working attempt to change to spdf and attach data to polygon feature

plot(thiessen)
points(points, pch=20)

which works, but some polygons have disappeared. Because of:

> thiessen <- intersect(border, v)
## Warning messages:
##1: In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
##  Ring Self-intersection at or near point 492420.43743714999 5371917.22911005
##2: In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
##  Ring Self-intersection at or near point 626630.43743726006 5315517.2291098097
##3: In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
##  Ring Self-intersection at or near point 598680.43743717996 5400017.2291099904
##4: In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
##  Ring Self-intersection at or near point 543250.43743715005 5416107.22911005

They are easy to find by doing

plot(border)

followed by sequences of (for a particular point)

points(492420.43743714999, 5371917.22911005, col='red')
zoom(border)

So you may need want to clean up your polygons first

Answer 2

There are places where your regions go in and out of the voronoi polygon edges, creating more features than you started with. EG in the the top:

The intersection will create some small features where the border line repeatedly crosses the straight Voronoi line.

Quite how you want to handle this depends on your application.