R simplify shapefile

Question

I have a shapefile with amazonic big rivers. The shapefile alone has 37.9 MB, together with attribute table it goes up to 42.1 MB. I'm generating PNG images of all brazilian Amazon, of 1260x940 pixels each, and all these data in the shapefile only slows down the drawing of each map, so I wanna simplify it.

gSimplify function, in rgeos package, seem to only simplify each polygon, not to get rid of the smaller ones. I tried it wih tolerance of 0.1 and 1000, and always I get length(shp@polygons) the same value: 27633. And the final plot takes almost the same time to draw. I need a function to which I tell that the final raster will be 1260x940 pixels, so it could remove every unnecessary point. Is there a function to do that?

Thanks in advance.

Answer 1

Pretty comprehensive solution here: http://www.r-bloggers.com/simplifying-polygon-shapefiles-in-r/

In summary, you need to get areas of your polygons:

area <- lapply(rivers@polygons, function(x) sapply(x@Polygons, function(y) y@area))

where rivers is you shapefile object in R.

Then you figure out the large polygons and keep them:

    sizeth <- 0.001 #size threshold of polygons to be deleted
    mainPolys <- lapply(area, function(x) which(x > sizeth))

    rivers@data <- rivers@data[-c(1:2),] 
    rivers@polygons <- rivers@polygons[-c(1:2)] 
    rivers@plotOrder <- 1:length(rivers@polygons)
    mainPolys <- mainPolys[-c(1:2)]

    for(i in 1:length(mainPolys)){   if(length(mainPolys[[i]]) >= 1 &&
    mainPolys[[i]][1] >= 1){
         rivers@polygons[[i]]@Polygons <- rivers@polygons[[i]]@Polygons[mainPolys[[i]]]
         rivers@polygons[[i]]@plotOrder <- 1:length(rivers@polygons[[i]]@Polygons)   } }

This might not be good enough, and you might not want to delete any polygons, in which case the dp() function from the shapefiles package will do the trick:

res <- 0.01 #the argument passed to dp() which determines extent of simplification. Increase or decrease as required to simplify more/less    
for(i in 1:length(rivers@polygons)){
      for(j in 1:length(rivers@polygons[[i]]@Polygons)){
        temp <- as.data.frame(rivers@polygons[[i]]@Polygons[[j]]@coords)
        names(temp) <- c("x", "y")
        temp2 <- dp(temp, res)
        rivers@polygons[[i]]@Polygons[[j]]@coords <- as.matrix(cbind(temp2$x, temp2$y))
      }
    }