[R-sig-Geo] Problems graphing shapefiles with detailed boundaries

[Roger Bivand]

On Thu, 19 Oct 2006 matt.pettis at thomson.com wrote:

> 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