modifying shpfiles in maptools
On Thu, 24 May 2007, Markus Loecher wrote:
Hello,
I am trying to manipulate a shpfile (actually subsetting it) but am
not getting anywhere. I would like to plot the census blocks of San
Francisco county (ID = 075) only but the shpfile I have is for all of CA.
So I am executing the following lines of code, where shpfile.CA is
the big one for CA:
shpfile.CA <- read.shape("bg06_d00.shp", dbf.data=TRUE, verbose=TRUE);
#keep only those entries for SF county:
shpfile.CA$att.data[,"COUNTY"] <- as.character(shpfile.CA$att.data[,"COUNTY"]);
i.keep <- shpfile.CA$att.data[,"COUNTY"] == "075";
shpfile$att.data <- shpfile.CA$att.data[i.keep,];
shpfile$Shapes <- shpfile.CA$Shapes[i.keep];
To me this seems like a correct subsetting of a list, but when I try
to plot the new, smaller shpfile, I get the error message:
Error in 1:attr(theMap$Shapes, "nshps") : NA/NaN argument
What foolish mistake am I committing ?
Also, it seems that plot.Map is deprecated. I wished I had an example
of how to use plot.Spatial instead, with a shpfile...
Yes, the Map object structure is too close to the input shapefile to be
easy to manipulate. Try rather:
shpfile.CA <- readShapePoly("bg06_d00.shp")
shpfile.CA$COUNTY <- as.character(shpfile.CA$COUNTY)
shpfile.SF <- shpfile.CA[shpfile.CA$COUNTY == "075",]
should work if shpfile.CA$COUNTY is correct - I'd do a
table(shpfile.CA$COUNTY) first to check that the values came out right.
The examples are in maptools (and rgdal) for reading and writing, and in
sp for manipulating and plotting - for easy plots see spplot().
Roger
Thanks! markus
_______________________________________________ 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