 
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.Bivand@nhh.no writes:

On Sat,  13 Jun 2009, Barry Rowlingson wrote:

> On Sat, Jun 13, 2009 at 5:31  PM, <BRWIN338@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]]

