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.

3
  • 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, 2016 at 16:02
  • Any chance of a full answer @mdsumner? That was a really useful tip
    – Phil
    May 9, 2016 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 9, 2016 at 14:42

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.