Roger Thank you for the code. I will be experimenting with it this week. Joe Barry I have downloaded and experimented with Quantum GIS. I am impressed but have not quite been able to delete the desired shapes yet. I have been able to edit and delete entire polygons but am struggling with deleting just the holes. I tried the delete vertex function but have not yet been able to get the desired results. I plan to continue learning more about Quantum GIS. Thank you for your prompt reply to my posting. Joe In a message dated 6/14/2009 10:51:50 A.M. Mountain Daylight Time, roger.biv...@nhh.no writes:
On Sat, 13 Jun 2009, Barry Rowlingson wrote: > On Sat, Jun 13, 2009 at 5:31 PM, <brwin...@aol.com> wrote: >> Good Morning: >> I would like to edit the county warning area boundaries shapefile >> "w_15jl09.shp" downloaded from: >> _http://www.weather.gov/geodata/catalog/wsom/html/cwa.htm_ >> (http://www.weather.gov/geodata/catalog/wsom/html/cwa.htm) >> >> Several of the shapefile's polygons contain extremely small "holes" that I >> would like to exclude from my plots. For example, I would like to remove >> slot 2 from the polygon for area 'RIW' (see below) and save the results in >> a new shapefile. > > You could either spend ages fiddling about manipulating sp spatial > data objects (which are lists of objects which are lists of vectors of > lists of coordinates, or similar), or you could download and install a > desktop GIS package and do it with a few clicks. Well, doing things programmatically rather than by editing vector data by hand does have its attractions, as in this case with hundreds of small holes in a global US data set. The recipe isn't visually attractive, but just uses [ls]apply() on lists, and slot() to get at the slots. I've chosen an area cutoff as 1.0-e07 here, which isn't very satisfactory, as these are geographical coordinates, so the area metric varies with latitude. In addition, I haven't checked back to make sure that there were no islands in the lakes that we're draining. Deleting all polygons less than this area will possibly remove "strategic assets" in the Pacific and elsewhere, but for display may be justified - do you need all those islands? Another consideration is that the rings defined not to be holes might be holes, or artefacts of some kind - you'd be pushed to land a seagull on an island of area 6.530e-12? # read in data and establish baseline: library(rgdal) w_15jl09 <- readOGR(".", "w_15jl09") summary(w_15jl09) polypick <- which(w_15jl09$WFO == "RIW") holes <- sapply(slot(w_15jl09, "polygons"), function(i) sapply(slot(i, "Polygons"), function(j) slot(j, "hole"))) areas <- sapply(slot(w_15jl09, "polygons"), function(i) sapply(slot(i, "Polygons"), function(j) slot(j, "area"))) by(unlist(areas), unlist(holes), summary) table(unlist(holes)) areas[[polypick]] holes[[polypick]] # the actual work begins here # get the top-level list of Polygons objects pls <- slot(w_15jl09, "polygons") # run lapply along it, checking the "Polygons" slot of each for # candidates for deletion; if any, remove: pls1 <- lapply(pls, function(i) { ipls <- slot(i, "Polygons") holes <- sapply(ipls, slot, "hole") areas <- sapply(ipls, slot, "area") ipls1 <- !(holes & (areas < 1.0e-07)) if(all(ipls1)) res <- i else res <- Polygons(ipls[ipls1], ID=slot(i, "ID")) res } ) # reassemble oSP <- SpatialPolygons(pls1, proj4string=CRS(proj4string(w_15jl09))) out <- SpatialPolygonsDataFrame(oSP, data=as(w_15jl09, "data.frame")) # export out with writeOGR, end of work # check where we landed - does it look OK? oholes <- sapply(slot(out, "polygons"), function(i) sapply(slot(i, "Polygons"), function(j) slot(j, "hole"))) oareas <- sapply(slot(out, "polygons"), function(i) sapply(slot(i, "Polygons"), function(j) slot(j, "area"))) by(unlist(oareas), unlist(oholes), summary) table(unlist(oholes)) oareas[[polypick]] oholes[[polypick]] - for RIW certainly, but the flattened representation used for Polygon objects within Polygons objects doesn't let us drop lakes together with their possibly contained islands (but not, I think, would Simple Features). If need be, the Polygon objects for dropping could be checked to see whether other Polygon objects' label points (centroids) fall within them, not bullet-proof, but could be done. I've tried to anticipate line breaks in the code example, it runs with this formatting for me on the specified file, but mailers do odd things, so line breaks may need adjusting. [ls]apply() and slot() only look forbidding until they save you hours, then they become really good friends! Roger > > With Quantum GIS for example, you can load a shapefile, enable > editing, click on points and delete them, even move them around > interactively, then save as a shapefile again. > > Quantum GIS is free, open-source, and crossplatform - www.qgis.org > for info. Other desktop GIS packages are available. > > Barry > **************An Excellent Credit Score is 750. See Yours in Just 2 Easy Steps! (http://pr.atwola.com/promoclk/100126575x1221322979x1201367215/aol?redir=http://www.freecreditreport.com/pm/default.aspx?sc=668072&hmpgID=62&bcd=Jun eExcfooterNO62) [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
