Re: [R-sig-Geo] thinning a SpatialPolygonsDataFrame
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) 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 density4; or proc gmap map=mymap(where=(density4)); ... 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 -- *** David PINAUD Ingénieur de Recherche Analyses spatiales Centre d'Etudes Biologiques de Chizé - CNRS UPR1934 79360 Villiers-en-Bois, France poste 485 Tel: +33 (0)5.49.09.35.58 Fax: +33 (0)5.49.09.65.26 http://www.cebc.cnrs.fr/ *** __ Information from ESET Mail Security, version of virus signature database 4621 (20091119) __ The message was checked by ESET Mail Security. http://www.eset.com attachment: pinaud.vcf___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] thinning a SpatialPolygonsDataFrame
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 - 50
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 density4;
or
proc gmap map=mymap(where=(density4));
...
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
Re: [R-sig-Geo] thinning a SpatialPolygonsDataFrame
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 - 50
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
[R-sig-Geo] thinning a SpatialPolygonsDataFrame
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 Streethttp://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
Re: [R-sig-Geo] thinning a SpatialPolygonsDataFrame
Hi Michael, maybe you should try the function dp() in the package shapefiles that is an implementation of the Douglas-Peucker polyLine simplification algorithm. 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 -- *** David PINAUD Ingénieur de Recherche Analyses spatiales Centre d'Etudes Biologiques de Chizé - CNRS UPR1934 79360 Villiers-en-Bois, France poste 485 Tel: +33 (0)5.49.09.35.58 Fax: +33 (0)5.49.09.65.26 http://www.cebc.cnrs.fr/ *** __ Information from ESET Mail Security, version of virus signature database 4618 (20091118) __ The message was checked by ESET Mail Security. http://www.eset.com attachment: pinaud.vcf___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
