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