[R-sig-Geo] thinning a SpatialPolygonsDataFrame

[Roger Bivand]

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
>>>>> 
>>>> 
>>>> 
>>> 
>> 
>> 
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no