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