[R-sig-Geo] thinning a SpatialPolygonsDataFrame

[Michael Friendly]

Thanks very much for doing the heavy lifting here, Roger.

For use in the Guerry package, I've turned this into a function (rather 
than saving other versions of gfrance).

If you include something like this in sp, I'll withdraw it from Guerry.

thinnedSpatialPoly <- function(SP, tolerance, minarea) {
  if (!require(shapefiles)) stop("shapefiles package is required")
  stopifnot(inherits(SP, "SpatialPolygons"))

    # TODO: determine and set defaults for tolerance, minarea 
    # TODO: suppress warnings: "In Polygon(crds_s) : Non-finite label 
point detected and replaced"
  pls <- slot(SP, "polygons")
  pls_dp <- vector(mode="list", length=length(pls))
  for (i in 1:length(pls)) {
    Pls <- slot(pls[[i]], "Polygons")
    Pls_dp <- vector(mode="list", length=length(Pls))
    for (j in 1:length(Pls)) {
      crds <- slot(Pls[[j]], "coords")
      crds_s <- dp(list(x=crds[,1], y=crds[,2]), tolerance=tolerance)
      crds_s <- do.call("cbind", crds_s)
      if(!identical(crds_s[1,], crds_s[nrow(crds_s),]))
        crds_s <- rbind(crds_s, crds_s[1,])
      Pls_dp[[j]] <- Polygon(crds_s)
    }
    Keep <- logical(length(Pls_dp))
    for (j in 1:length(Pls_dp)) {
      Keep[j] <- TRUE
      if (slot(Pls_dp[[j]], "area") < minarea) Keep[j] <- FALSE
    }
    Pls_dp <- Pls_dp[Keep]
    pls_dp[[i]] <- Polygons(Pls_dp, ID=slot(pls[[i]], "ID"))
  }
    SP_dp <- SpatialPolygons(pls_dp, proj4string=slot(SP, "proj4string"))
  if(inherits(SP, "SpatialPolygonsDataFrame")) {
    data <- slot(SP, "data")
    SP_dp <- SpatialPolygonsDataFrame(SP_dp, data=data)
  }
  SP_dp
}

best,
-Michael


Roger Bivand wrote:
> On Thu, 19 Nov 2009, Pinaud David wrote:
>
>> Hi Michael,
>> Roger is fully right, this function does not preserve the topology, 
>> so be aware of that, some problems can occur...
>> If you want to use shapefiles::dp() just for raw plotting and visual 
>> simplification, you can try (on the fly):
>>
>> library(rgdal)
>> library(shapefiles)
>>
>> fr <- readOGR("polygonFRA_WGS84.shp", "polygonFRA_WGS84") # a 
>> shapefile of France with complex topology (holes, islands...) in 
>> WGS84 coordinates
>> pp <- slot(fr, "polygons") # take the polygons
>> cf <- coordinates(slot(pp[[39]], "Polygons")[[1]])  # extract the 
>> coordinates of the main polygon ("continental France")
>> pf <- list(x=cf[,1], y=cf[,2]) # list of coordinates, as dp() needs a 
>> list and not a matrix or dataframe...
>> cf1 <- dp(pf, 0.1) # simplification, with a bandwith of 0.1 decimal 
>> degree
>> plot(fr)  # to see the result
>> points(cf1, col="red", t="l")
>
> This seems to work OK given that slivers are not important:
> library(sp)
> library(Guerry)
> # installed from R-Forge
> data(gfrance)
> object.size(gfrance)
> pls <- slot(gfrance, "polygons")
> pls_dp <- vector(mode="list", length=length(pls))
> require(shapefiles)
> tol <- 2500
> minarea <- 500000
> for (i in 1:length(pls)) {
>   Pls <- slot(pls[[i]], "Polygons")
>   Pls_dp <- vector(mode="list", length=length(Pls))
>   for (j in 1:length(Pls)) {
>     crds <- slot(Pls[[j]], "coords")
>     crds_s <- dp(list(x=crds[,1], y=crds[,2]), tolerance=tol)
>     crds_s <- do.call("cbind", crds_s)
>     if(!identical(crds_s[1,], crds_s[nrow(crds_s),]))
>       crds_s <- rbind(crds_s, crds_s[1,])
>     Pls_dp[[j]] <- Polygon(crds_s)
>   }
>   Keep <- logical(length(Pls_dp))
>   for (j in 1:length(Pls_dp)) {
>     Keep[j] <- TRUE
>     if (slot(Pls_dp[[j]], "area") < minarea) Keep[j] <- FALSE
>   }
>   Pls_dp <- Pls_dp[Keep]
>   pls_dp[[i]] <- Polygons(Pls_dp, ID=slot(pls[[i]], "ID"))
> }
> gfrance_dp <- SpatialPolygonsDataFrame(SpatialPolygons(pls_dp), 
> data=slot(gfrance, "data"))
> object.size(gfrance_dp)
>
> If this was a function, the input would be an object inheriting from 
> SpatialPolygons, the tolerance, and a minimum polygon area. The 
> removal of small Polygon objects needs protecting from removing all 
> belonging to a given Polygons object. The CRS also needs copying 
> across. The objects are rebuilt to correct areas, centroids, plot 
> orders, and the bounding box, all of which may change. For figures, 
> the companion thread on R-help is relevant, PDF output for choropleth 
> maps is often a good deal larger in file size than that of the 
> equivalent PNG device.
>
> I'll ask the shapefiles authors whether I can copy dp() to maptools 
> and include suitable methods - this will help in benchmarking the 
> future GEOS version.
>
> Hope this helps,
>
> Roger
>
>>
>> HTH
>> David
>>
>> Michael Friendly a écrit :
>>> I think, for my application, I'd be happy with the D-P polyline 
>>> simplification algorithm, because that is what is used
>>> in SAS (worked well), and I don't think there are unusual topologies 
>>> in my map of France that would be so severely
>>> out of whack as to lead to *significant* visual artifacts.  In fact, 
>>> you might well expect some artifacts from any visual
>>> thinning, but it's a matter of the tradeoff in the way the thinned 
>>> map is used in a visualization.  Mark Monmonier's
>>> US Visibility Map might be an extremely thinned, but highly useful 
>>> example.
>>>
>>> For R spatial analysis, I think this is worth pursuing and 
>>> integrating into sp methods.  In SAS, proc greduce works simply
>>> by adding another variable, density, to the (x,y) coordinates of the 
>>> spatial polygons, density %in% 1:5, where density==5
>>> is the full map.  It is then a simple matter to subset the polygon 
>>> outlines by saying
>>>
>>> data smallmap;
>>>    set mymap;
>>>    where density<4;
>>>
>>> or
>>>
>>> proc gmap map=mymap(where=(density<4));
>>> ...
>>>
>>> Meanwhile, I can't see easily how I could use shapefiles::dp() to 
>>> thin my Guerry::gfrance maps, because the documentation is,
>>> shall we say, somewhat thin.
>>> -Michael
>>>
>>>
>>>
>>> Roger Bivand wrote:
>>>> On Wed, 18 Nov 2009, Pinaud David wrote:
>>>>
>>>>> Hi Michael,
>>>>> maybe you should try the function dp() in the package shapefiles 
>>>>> that is an implementation of the Douglas-Peucker polyLine 
>>>>> simplification algorithm.
>>>>
>>>> Note that its help page does warn that it is not 
>>>> topology-preserving, that is that lines are generalised, but that 
>>>> coincident lines (boundaries of neighbouring polygons) may be 
>>>> generalised differently. GEOS offers a topology-preserving line 
>>>> generalisation facility, which ought to take longer but do better 
>>>> than dp(), because it will not lead to visual artefacts 
>>>> (overlapping polygons, interpolygon slivers, etc.).
>>>>
>>>> Roger
>>>>
>>>>> HTH
>>>>> David
>>>>>
>>>>> Michael Friendly a écrit :
>>>>>> The Guerry package contains two maps of france (gfrance, 
>>>>>> gfrance85) which are quite detailed and large in size (~900K).
>>>>>> In writing a vignette for the package, there are quite a few 
>>>>>> figures that use the map multiple times in a layout, and
>>>>>> consequently result in huge file sizes for the .PDF files 
>>>>>> created.  For these purposes, the map need not be nearly
>>>>>> so detailed.
>>>>>>
>>>>>> I'm wondering if there is a facility to "thin" the map by drawing 
>>>>>> it at a lower density of lines in the polygon regions.
>>>>>> When I was working with SAS, there was a GREDUCE procedure that 
>>>>>> did this nicely.
>>>>>>
>>>>>> thanks,
>>>>>> -Michael
>>>>>>
>>>>>
>>>>>
>>>>
>>>
>>>
>>
>>
>


-- 
Michael Friendly     Email: friendly AT yorku DOT ca 
Professor, Psychology Dept.
York University      Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Street    http://www.math.yorku.ca/SCS/friendly.html
Toronto, ONT  M3J 1P3 CANADA