#  Given .shp file (Arcview 8.x or earlier), this R code creates:

#    an areal adjacency matrix using maptools,  

#    an edge adjacency matrix (indicating which edges touch each other)

#    polygon or edge maps (by calling samplePlot.txt)

#  This sample code produces zip code level maps for northern Minnesota, 1999

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


#setwd("F:\\Boundary\\PlotCode\\GetEdgeNeig\\tutorial\\")

## output files will be saved under current directory 

## save "matchzip.exe" and "edgeneig.ext" under directory "edgefolder"  

library(maptools)
library(spdep)

map<-read.shape("shapefiles\\smdc.shp") 
polylist<- Map2poly(map)
plot.polylist(polylist,axes=F)


################################################
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%##
##%  Use maptools to get area adj. matrix    %##
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%##
################################################

nb.r<-poly2nb(polylist, queen=F)    #two areas are neighbors if share common edges with length >0
mat<-nb2mat(nb.r, style="B")        #mat is the 0/1 adjacency matrix
n.site<-dim(mat)[1]                 #n.site: number of areas
n.edge<-sum(mat)/2                  #n.edge: number of unique pairs (ij and ji only count once).

SEind1<-SEind2<-0
matmy<-mat
for(i in 1:(n.site-1)){
for(j in (i+1):n.site){
if (mat[i,j]>0) {SEind1<-c(SEind1,i)
        SEind2<-c(SEind2,j)
        matmy[i,j]<-matmy[j,i]<-length(SEind1)-1}
}
}
SEind1<-SEind1[-1]                  #edges are sorted by row of upper triangle of the Adj. matrix
SEind2<-SEind2[-1]                  #SEind1[k]=i and SEind2[k]=j => kth edge is edge ij

#dput(SEind1,"SEind1.txt")
#dput(SEind2,"SEind2.txt")
#dput(mat,"W.txt")

## creat adjacency information needed for WinBUGS
mkAdj<-function(W){
n<-nrow(W)
adj<-0
for(i in 1:n){
    for(j in 1:n){  if(W[i,j]==1){adj<-append(adj,j)} }
        }
adj<-adj[-1]
return(adj)
}
#dput(mkAdj(W),"Sadj.txt")
#dput(as.vector(rowSums(W)),"Snum.txt")


################################################
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%##
##%   Create adj. matrix for the edges       %##
##%   and prepare plotting code for edges    %##
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%##
################################################

## 1. prepare needed files and save them ##
##    under the directory where you have ##
##    matchzip.exe and edgeneig.ext      ##


    # Dump out the coordinates (by polygon) #
    ## Defult order for polygons is by first column of map$att.data ##
    for(i in 1:n.site){
    write(t(map$Shapes[[i]]$verts), paste("edgefolder\\",i,".txt",sep=""), ncolumns=2)
    }
    
    # Dump out SEind, the site-edge correspondence table #
    write(rbind(SEind1, SEind2), paste("edgefolder\\SEind.txt",sep=""), ncolumns=2 )


## 2.   Double click  "matchzip.exe". A dos window will pop up. Type in SEind.txt.##
##  This will produce many (n.edge) files at the current directory.           ##
##  Then use the following code to prepare the edge plotting code.            ##
##  "edgelines" should be dumped out and called in later when make edge plots.##

    edgelines<-vector(mode="list",length=n.edge)
    for(i in 1:n.edge){
    edgelines[[i]]<-read.table(paste("edgefolder\\output-",i,".txt",sep=""),header=F,na.strings="*")
    }

    dput(edgelines,"edgelines.txt")
    
    # Then use  "lines(as.matrix(edgelines[[i]]),col=y.shad[i],lwd=3)" to plot edges.
    # Two functions that can produce polygon maps and edge maps are given in "samplePlot.R".
    # See part 6 of this file for examples.


## 3. Double click "edgeneig.exe". A dos window will pop up. Type in SEind.txt ##
##  Two files will be produced (may take a while): ##
##  file 1 is the upper triangular of the W matrix for the edges;                        ##
##  file 2 is the number of 1s for each row of the upper triangular matrix.              ##


##  4. Produce the Wstar matrix which gives the neighborhood structure of the edges. ##

    tempn<-scan(paste("edgefolder\\file2.txt",sep=""))
    tempneig<-scan(paste("edgefolder\\file1.txt",sep=""))

    Wstar<-matrix(0,nrow=n.edge,ncol=n.edge)
    start<-1
    end<-0
    for(i in 1:n.edge){
    if (tempn[i]>0) { 
    start<-end+1
    end<-start+tempn[i]-1
    neig<-tempneig[start:end]
    l<-tempn[i]
    for (k in 1:l) {Wstar[i,neig[k]]<-Wstar[neig[k],i]<-1}
    }
    } 

    #dput(Wstar,"Wstar.txt")

## 5. Remove files not needed any more. ##

   for(i in 1:n.edge){
   unlink(paste("edgefolder\\",i,".txt",sep=""))
   unlink(paste("edgefolder\\output-",i,".txt",sep=""))
    }
    
   unlink(paste("edgefolder\\file1.txt",sep=""))
   unlink(paste("edgefolder\\file2.txt",sep=""))
   unlink(paste("edgefolder\\SEind.txt",sep=""))


## 6. Examples for plotting polygons and edges. ##

	source("samplePlot.txt")
	x<-map$att.data$POP1999 
	PolyPlot(map, "Example map: Population 1999", x , 4, bk="c", cuts=c(min(x),100, 1000, 3500, max(x)))

	delta.calculate<-function(x,ind1, ind2,n){
	delta<-rep(0,n)
	for (k in 1:n){
	i<-ind1[k]
	j<-ind2[k]
	delta[k]<-abs(x[i]-x[j])}
	return(delta)
	}

	dx<-delta.calculate(x, SEind1, SEind2, n.edge)
	fig<-expression(paste(Delta[x[ij]]," = |",x[i]," - ", x[j], "|",sep=""))
	EdgePlot(map, edgelines,fig,dx, 4, bk="c", quantile(dx, prob=c(0,0.5,0.7,0.9,1)))