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
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo