Le portail francophone de la géomatique

### loading data
bh <- readShapePoly(system.file("etc/shapes/bhicv.shp",
      package="spdep")[1])
### data standardized 
dpad <- data.frame(scale(bh@data[,5:8]))
 
### neighboorhod list
bh.nb <- poly2nb(bh)
 
### calculating costs
lcosts <- nbcosts(bh.nb, dpad)
 
### making listw
nb.w <- nb2listw(bh.nb, lcosts, style="B")
 
### find a minimum spanning tree
mst.bh <- mstree(nb.w,5)
 
### the mstree plot
par(mar=c(0,0,0,0))
plot(mst.bh, coordinates(bh), col=2,       
     cex.lab=.7, cex.circles=0.035, fg="blue")
plot(bh, border=gray(.5), add=TRUE)
 
### three groups with no restriction
res1 <- skater(mst.bh[,1:2], dpad, 2)
 
### groups size
table(res1$groups)
 
### the skater plot
par(mar=c(0,0,0,0))
plot(res1, coordinates(bh), cex.circles=0.035, cex.lab=.7)
 
### the skater plot, using other colors
plot(res1, coordinates(bh), cex.circles=0.035, cex.lab=.7,
     groups.colors=rainbow(length(res1$ed)))
 
### the Spatial Polygons plot
res1 <- skater(mst.bh[,1:2], dpad, 9)
 
plot(bh, col=(1:9)[res1$groups])
 
### EXPERT OPTIONS
 
### more one partition
res1b <- skater(res1, dpad, 1)
 
### length groups frequency
table(res1$groups)
table(res1b$groups)
 
### thee groups with minimum population 
res2 <- skater(mst.bh[,1:2], dpad, 2, 200000, bh@data$Pop)
 
### thee groups with minimun number of areas
res3 <- skater(mst.bh[,1:2], dpad, 2, 3, rep(1,nrow(bh@data)))
 
### thee groups with minimun and maximun number of areas
res4 <- skater(mst.bh[,1:2], dpad, 2, c(20,50), rep(1,nrow(bh@data)))
 
table(res2$groups)
table(res3$groups)
table(res4$groups)
 
### if I want to get groups with 20 to 40 elements
res5 <- skater(mst.bh[,1:2], dpad, 2,
   c(20,40), rep(1,nrow(bh@data))) ## DON'T MAKE DIVISIONS 
table(res5$groups)
 
### In this MST don't have groups with this restrictions
### In this case, first I do one division
### with the minimun criteria
res5a <- skater(mst.bh[,1:2], dpad, 1, 20, rep(1,nrow(bh@data))) 
table(res5a$groups)
 
### and do more one division with the full criteria
res5b <- skater(res5a, dpad, 1, c(20, 40), rep(1,nrow(bh@data)))
table(res5b$groups)
 
### and do more one division with the full criteria
res5c <- skater(res5b, dpad, 1, c(20, 40), rep(1,nrow(bh@data)))
table(res5c$groups)
 
### It don't have another divison with this criteria
res5d <- skater(res5c, dpad, 1, c(20, 40), rep(1,nrow(bh@data)))
table(res5d$groups)