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
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@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.biv...@nhh.no
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
