# Combining postal codes into sales territories by value using R and ggplot2?

I have the shapefile for Canadian census subdivisions. For each subdivision I calculated some value (someValue). I want to combine census subdivisions into territories in such a way that territories have similar values (total sum of someValue for subdivisions within a territory) and that these territories have no gaps/breaks in them.

What I tried was:

• use distance matrix
• create list of candidate census subdivision based on proximity to some seed value
• iterate until the target value for territory is reached

I don't believe that this an optimal solution

Below is a sub-optimal solution based on simple x-y coordinates. I've since updated the code to check whether or not each potential candidate is adjacent to an already existing territory, but it did not improve the solution. Two criteria that I have is terr.size: I want the size of an average territory to be small, and terr.sddev that measures dispersion of territories. I create many runs of the code below and pick the iteration with the smallest terr.size and terr.sddev

``````library(ggplot2)

lon<-rep(1:10,each=10)
lat<-rep(1:10,10)
someValue<-rnorm(100, mean = 20, sd = 5)
dataset<-data.frame(lon,lat,someValue)
dataset<-cbind(CCSD=row.names(dataset),dataset)

ref.dataset<-dataset

dmatrix<-as.matrix(dist(dataset[,2:3]))
rownames(dmatrix)<-dataset\$CCSD

#total control units (CCSD)
tCCSD=1000
#assume 100 territories
nter=100
#likely average number of CCSD by territory plus 20 percent
aCCSD=1.2*(tCCSD/nter)
#value of of a territory
vter=200
#range of 10%: sufficient to consider vter is close to 200
rter=.1

#terr list
terr.out<-list()
#start at random value, i.e. 14
firstCCSD<-1

for (i in 1:nter)

{
#list of candidates, adjust the lenght SO NOT OUT OF RANGE
terr.range<-min(nrow(dataset),ceiling(aCCSD)) # to account that there might not be enough CCSDs left at the end
#print(terr.range)
terr.index<-firstCCSD
#print(terr.index)
candidates<-as.numeric(names(sort(dmatrix[,firstCCSD])))[1:terr.range]
#print(candidates)
#print(length(candidates))
terr.value<-subset(dataset,CCSD==candidates,4)
#terr.value<-dataset[candidates,3]
#print(terr.value)

for (j in 1:terr.range){
if(terr.value<=vter*(1-rter)){
#print(terr.value)
terr.value<-terr.value+subset(dataset,CCSD==candidates[j+1],4)
#print(terr.value)
terr.index<-c(terr.index,candidates[j+1])
#print(terr.index)
nextCCSD<-candidates[j+2]
#print(nextCCSD)
j=j+1
}
}

#create indexies and add to list
terr.value #total value of the territory
terr.index #list of CCSDs in territory
terr.out[[i]]<-list(terr.index,terr.value,length(terr.index))

#exclude used CCSD for the next step
firstCCSD<-nextCCSD
#print(firstCCSD)
dmatrix<-dmatrix[!dataset\$CCSD %in% terr.index,]
#print(dmatrix)
dataset<-dataset[!dataset\$CCSD %in% terr.index,]
#print(dataset)
i=i+1
}

result<-data.frame(CCSD=1,territory=1)
result<-result[-1,]

for (index in 1:length(terr.out)){
CCSD<-unlist(sapply(terr.out,function(x) x)[index])
territory<-rep(index,length(CCSD))
result<-rbind(result,cbind(CCSD,territory))
index=index+1
}

ref.dataset1<-merge(ref.dataset,result)

#SUMMARY STATS
#distro of territories
aggregate(ref.dataset1\$someValue,list(ref.dataset1\$territory),FUN=sum)
#stdev to be used as selection criteria
terr.sddev<-
sd(aggregate(ref.dataset1\$someValue,list(ref.dataset1\$territory),FUN=sum)\$x)
#range for width and length of territories
terr.dims<-rowMeans(sapply(unique(ref.dataset1\$territory),function(x) {

lat_temp<-max(ref.dataset1[ref.dataset1\$territory==x,]\$lat)-min(ref.dataset1[ref.dataset1\$territory==x,]\$lat)
lon_temp<-max(ref.dataset1[ref.dataset1\$territory==x,]\$lon)-min(ref.dataset1[ref.dataset1\$territory==x,]\$lon)
result<-c(lat_temp, lon_temp)
#result<-lat_temp
}
)
)

#print(terr.sddev)
terr.size<-sqrt(terr.dims^2+terr.dims^2)
#print(terr.size)

#REPLACE THE NAME
ref_temp<-data.frame(
ref.dataset1,
TerrBalance=rep(terr.sddev, nrow(ref.dataset1)),
TerrSize=rep(terr.size, nrow(ref.dataset1))
)
#table(result\$CCSD)

qplot(lon, lat, data=ref.dataset1,  color=as.factor(territory),size=5)
``````
• Define "optimal". Fastest? One that creates most compact territories (define "compact")? Do you have code for a sub-optimal solution that you can share? – Spacedman Jun 13 '17 at 15:54
• @Spacedman sub-optimal code and criteria for success are added to the original post. My hope that there is a "k-means"-type solution instead of iterative approach that I'm attempting – Konstantin Mingoulin Jun 13 '17 at 17:57