Problems graphing shapefiles with detailed boundaries

Hi,

I'm loading in a shapefile with all of the voting precincts in
Minnesota.  Each precinct belongs to a Senate district, which has a
variable named "SEN" in the .dbf  ("Precinct" is also in the .dbf).

I am including a script that should allow anyone to replicate my
problems, which I will lay out here (on Windows XP, R 2.3.1):

1.  I limit myself to the first 2 SEN districts.  If I subset to each
district and plot them with spplot, they individually get drawn and
colored just fine.  When I try to draw them together, I get a
checkerboard effect in the colorization (same with just drawing the
internal boundaries).
This is a simple mistake in subsetting, you have to specify the chosen 
values in an appropriate way:

x12 <- vtd[vtd[["SEN"]] == 1 | vtd[["SEN"]] == 2, "SEN"]
spplot(x12)

2.  When I try to remove internal boundaries from the SEN districts, it
works fine on SEN == 1, but I get some phantom lines in SEN == 2.
I'm not sure whether the above clears up the problems, I'd also use 
vtd$SEN rather than vtd[["SEN"]], less error prone unless you need to 
insert the fiel/column name as a string.

I can also see the internal debris in vtd$SEN == 2, those are slivers 
from the original shapefile:

x2 <- vtd[vtd[["SEN"]]  == 2,]
x <- unionSpatialPolygons(x2, x2[["SEN"]])
plot(x)
summary(sapply(slot(slot(x, "polygons")[[1]], "Polygons"), function(x) 
  slot(x, "area")))

is:

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
7.485e-05 2.075e-02 9.230e-02 2.593e+08 3.406e-01 1.996e+10 

so pruning the output object to leave only the large polygon is an option.

Hope this helps,

Roger

I wager these two problems might be related, but I am not sure.  Any
ideas?

thanks,
Matt

########################################################################
####
#----------------------------------------------------------------------
# Config
#----------------------------------------------------------------------
# When I need to turn off proxy
# Sys.putenv("no_proxy"=True)
library(sp)
library(maptools)
library(spgpc)
# needed for unionSpatialPolygons?  Looks like not...
#library(gpclib)

#----------------------------------------------------------------------
# Load Shapefile
#----------------------------------------------------------------------
# This shapefile is the shapefile of all of the voting precincts in
Minnesota, 2006
# All Precincts are contained within a SEN group
# shapefile source:
ftp://ftp.commissions.leg.state.mn.us/pub/gis/shape/vtd2006.zip
# unzipped to pdir\Data\vtd2006
pdir <- "C:\\Documents and Settings\\U0008998\\My
Documents\\Projects\\GISVote";
vtd <- readShapePoly(paste(pdir,"\\Data\\vtd2006\\vtd2006.shp",sep=""))

#----------------------------------------------------------------------
# Plot the whole map of Minnesota
#----------------------------------------------------------------------
plot(vtd)

#----------------------------------------------------------------------
# Shading regions
#----------------------------------------------------------------------
# Show and color SEN == 1, the Northwest corner of Minnesota
x <- vtd[vtd[["SEN"]] == c(1), "SEN"]
spplot(x,col.regions=2)

# Show and color SEN == 2, east and south of SEN == 1
x <- vtd[vtd[["SEN"]] == c(2), "SEN"]
spplot(x,col.regions=3)

# Plot them both together
x <- vtd[vtd[["SEN"]] == c(1,2), "SEN"]
spplot(x,col.regions=2:3)

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# Note that the whole regions are not colored
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#----------------------------------------------------------------------
# Plotting SEN districts together
#----------------------------------------------------------------------
# Subset to SEN == 1 and plot
x <- vtd[vtd[["SEN"]]  == 1,names(vtd)]
plot(x)

# Subset to SEN == 2 and plot
x <- vtd[vtd[["SEN"]]  == 2,names(vtd)]
plot(x)

# Plot SEN == 1 and 2 together
x <- vtd[vtd[["SEN"]]  == 1:2,names(vtd)]
plot(x)
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# Note that the whole regions are not plotted
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#----------------------------------------------------------------------
# Removing internal Precinct boundaries from SEN regions
#----------------------------------------------------------------------
# Subset to SEN == 1 and plot
x <- vtd[vtd[["SEN"]]  == 1,names(vtd)]
x <- unionSpatialPolygons(x, x[["SEN"]])
plot(x)

# Subset to SEN == 2 and plot
x <- vtd[vtd[["SEN"]]  == 2,names(vtd)]
x <- unionSpatialPolygons(x, x[["SEN"]])
plot(x)
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# Note that some internal boundaries are not removed from SEN == 2
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

# Subset to SEN == 1 and 2 and plot
x <- vtd[vtd[["SEN"]]  == 1:2,names(vtd)]
x <- unionSpatialPolygons(x, x[["SEN"]])
plot(x)
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# Note that most internal boundaries are not removed
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

########################################################################
####

	[[alternative HTML version deleted]]

_______________________________________________
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

Problems graphing shapefiles with detailed boundaries

Thread (5 messages)