[R-sig-Geo] Editing a shapefile

[Roger Bivand]

On Sat, 13 Jun 2009, Barry Rowlingson wrote:

> On Sat, Jun 13, 2009 at 5:31 PM, <BRWIN338 at 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
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.Bivand at nhh.no