[GIS] writeOGR with a spatialpolygon simplified by gSimplify

rrgdalshapefilesimplify

I'm using gSimplify (rgeos package) to simplify the geometries of a shapefile. The funcion works good, but now I can't write the output in a new shapefile.
I tried some ways:

writeOGR(simplyshape, file, driver="ESRI Shapefile", layer='test')

I get

obj must be a SpatialPointsDataFrame, SpatialLinesDataFrame or
SpatialPolygonsDataFrame

and with:

writePolyShape(simplyshape, file)

I get:

Error: is(x, "SpatialPolygonsDataFrame") is not TRUE

Best Answer

Coerce your object to the appropriate Spatial*DataFrame-class (Points/Lines/Polygons), e.g. for SpatialPolygons using as(x, "SpatialPolygonsDataFrame" ):

R> l <- readWKT("LINESTRING(0 7,1 6,2 1,3 4,4 1,5 7,6 6,7 4,8 6,9 4)")
R> x1 <- gSimplify(p, tol=10)
R> class(x1)
[1] "SpatialPolygons"
attr(,"package")
[1] "sp"
R> x2 <- as(x, "SpatialPolygonsDataFrame")
R> class(x2)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

Related Solutions

[GIS] Merging and dissolving shapefile in R

I think we should do this bit by bit. It's good to keep in mind that the simpler and more minimal you keep the example, the more likely you are to get responses. Here's my first shot. I don't see anything obviously wrong with it so far - it produces a plot. However, because you've provided only a subset of the data, the plot is not contiguous (see below). Why not try this with the full data set and see whether it works for you.

spanish counties

Note that the fortified data frame looks okay as well.

> str(foo)
'data.frame':   7876 obs. of  7 variables:
 $ long : num  168544 169128 169380 170315 170700 ...
 $ lat  : num  4814368 4814211 4813874 4814704 4816283 ...
 $ order: int  1 2 3 4 5 6 7 8 9 10 ...
 $ hole : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ piece: Factor w/ 49 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ group: Factor w/ 138 levels "2003.1","2003.2",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ id   : chr  "2003" "2003" "2003" "2003" ...

Here's the code.

library(rgeos)
library(ggplot2)
work.dir <- "your_dir_no_trailing_slash"
zones.sp <- readOGR(dsn = work.dir, layer = "SP0111")
zones.data <- read.csv("statsz.csv")
zones.spwd <- zones.sp
zones.spwd@data <- merge(zones.sp@data, zones.data, by = "PROVMUN", all = TRUE)
zones.IDs <- as.numeric(zones.spwd[["CODER"]])
# note that we only have a subset of the data so only a few zones...
which(!is.na(zones.IDs))
# this seems to work
regions.sp <- gUnaryUnion(zones.spwd, id = zones.IDs)
str(regions.sp)
plot(regions.sp)
foo <- fortify(regions.sp)

R – Dissolving/Unifying Ill-Behaved Polygons in R

Without your original data, I can't be sure this will work, but I thought it might help you out. I didn't bring it all the way there, this solution still likely needs some level of automation, but might give you a general way forward

First, I create some spatial polygons

polypoints1 <- matrix(c(1,2,2,1,1,2,2,1,1,2),ncol=2)
polypoints2 <- matrix(c(1,3,3,1,1,3,3,1,1,3),ncol=2)
polypoints3 <- matrix(c(1,2,2,1,1,2,2,1,1,2)+1.1,ncol=2)
polypoints4 <- matrix(c(1,2,2,1,1,2,2,1,1,2)+0.5,ncol=2)

p1 <- Polygon(polypoints1)
ps1 <- Polygons(list(p1),1)
sps1 <- SpatialPolygons(list(ps1))

p2 <- Polygon(polypoints2)
ps2 <- Polygons(list(p2),2)
sps2 <- SpatialPolygons(list(ps2))

p3 <- Polygon(polypoints3)
ps3 <- Polygons(list(p3),3)
sps3 <- SpatialPolygons(list(ps3))

p4 <- Polygon(polypoints4)
ps4 <- Polygons(list(p4),4)
sps4 <- SpatialPolygons(list(ps4))

I plotted them just to see

plot(sps2,col='green')
plot(sps1,add=T,col='blue')
plot(sps3,add=T,col='yellow')
plot(sps4,add=T,col='purple')

I merged them into an spdf

data=data.frame(c(x=rep(1,4)),row.names=c(1:4))
sps <- SpatialPolygons(list(ps1,ps2,ps3,ps4))
spdf <- SpatialPolygonsDataFrame(sps,data)

You can identify which polygon overlaps which like so:

gIntersects(spdf,spdf,byid =T)

From the above command you could create some kind of loops to do the overlapping combinations below (I'm just ignoring sps4 for brevity at this point)

poly2a <- gIntersection(spdf[2,],spdf[1,],drop_lower_td=T)
poly2a <- SpatialPolygonsDataFrame(poly2a,data.frame(c(x=1),row.names=c(1)))
plot(poly2a,add=T,col='red')

This time we need to change the ID since we're going to rbind these later

poly2b <- gIntersection(spdf[2,],spdf[3,],drop_lower_td=T)
poly2b <- spChFIDs(poly2b,"2")
poly2b <- SpatialPolygonsDataFrame(poly2b,data.frame(c(x=1),row.names=c(2)))
plot(poly2b,add=T,col='red')

Merge the overlapping polygons into another spdf

spdf_overlaps <- rbind(poly2a,poly2b)
poly2 <- unionSpatialPolygons(spdf_overlaps,rep(1,2))
plot(poly2,add=T,col='blue')

Now we have poly2 which is where we have 2 layers overlapping (except combinations with sps4) then to figure out 3 layers, we just have to check out where poly2 and spdf overlap (if you make a more automated version of this, you'll need to make sure that 'poly2' does not include sps4 as in this example)

gIntersects(poly2,spdf[4,],byid =T)

poly3 <- gIntersection(poly2,spdf[4,],drop_lower_td=T)
plot(poly3,add=T,col="red")

Check it out

gIsValid(poly2)
gIsValid(poly3)

Alternatively, you could always do a pseudo rasterization, much easier, but you loose some detail depending on your cell size:

First make the grid:

bb <- bbox(spdf)
cs <- c(0.1,0.1)  # cell size
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)


sp_grd <- SpatialGridDataFrame(grd,
                           data=data.frame(id=1:prod(cd)))

Then, make grid into a polygon which used for overlap

library(Grid2Polygons)
grid <- Grid2Polygons(sp_grd)
plot(grid)

Then count the number of polygons that overlap each grid cell

count <- apply(gContains(spdf,grid,byid=T),1,sum)

Finally, plot it!

plot(grid)
for(i in 1:length(grid)){
    plot(grid[i,],col=rev(heat.colors(3))[count[i]],add=T)
}

Best Answer

Related Solutions

[GIS] Merging and dissolving shapefile in R

R – Dissolving/Unifying Ill-Behaved Polygons in R

Related Question