2

I'm trying to use the method proposed by the @ialm in the Question "How to speed up the plotting of polygons in R?"

There, he loads the maps using getData from the raster package, I tested and it works fine.

I downloaded files from the Brazilian Institute of Geography and Statistics - IBGE and one of the ShapeFiles could be this one, then I use readOGR from rgdal.

I'm facing an error on method 2

  # Get the main polygons, will determine by area.
getSmallPolys <- function(poly, minarea=0.01) {
  # Get the areas
  areas <- lapply(poly@polygons, 
                  function(x) sapply(x@Polygons, function(y) y@area))

  # Quick summary of the areas
  print(quantile(unlist(areas)))

  # Which are the big polygons?
  bigpolys <- lapply(areas, function(x) which(x > minarea))
  length(unlist(bigpolys))

  # Get only the big polygons
  for(i in 1:length(bigpolys)){
    if(length(bigpolys[[i]]) >= 1 && bigpolys[[i]][4] >= 1){
      poly@polygons[[i]]@Polygons <- poly@polygons[[i]]@Polygons[bigpolys[[i]]]
      poly@polygons[[i]]@plotOrder <- 1:length(poly@polygons[[i]]@Polygons)
    }
  }
  return(poly)
}

Particularly here if(length(bigpolys[[i]]) >= 1 && bigpolys[[i]][4] >= 1)

Error in if (length(bigpolys[[i]]) >= 1 && bigpolys[[i]][4] >= 1) { : 
  missing value where TRUE/FALSE needed

As for me bigpolys[[i]][4] is always NA, so I started to suspect the structure of files are not the same.

So far, as starter in this area, I have tried a combination of query methods ( slots, head, str, attributes, class, unclass ) trying to uncover the differences that would justify the error in my case.

> class(mex)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
> class(mapa)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

> slotNames(mapa)
[1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"
> slotNames(mex)
[1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"

I need help in 2 directions, what are the best methods to compare 2 shapefiles like in the example above and detect their structural differences, and if possible an appropriate solution for clipping the small polygons.

user1873410
  • 121
  • 1
  • 1
    I would use rgeos::gArea(mex, byid = TRUE) to get a vectorized area that you can use as an index against the spdf rows. I.e. mex[rgeos::gArea(mex, byid = TRUE) < 1, ]. If you need to decompose those objects first use ?disaggregate. – mdsumner Apr 24 '16 at 16:02
  • Any chance of a full answer @mdsumner? That was a really useful tip – Phil May 09 '16 at 12:59
  • Remind me again on Wednesday, or maybe the week after . . . more seriously we should put up content like this on a github project, use the issues and wiki tools to capture real content – mdsumner May 09 '16 at 14:42

0 Answers0