0

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[1],4)
#terr.value<-dataset[candidates[1],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[1])[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[1]^2+terr.dims[2]^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
-1

We have an open source algorithm called Territorium available here that does this and it's also integrated into our open source desktop app here. Specifically the Territorium algorithm takes the number of territories as an input and optimises to (a) keep the sum of quantity within a min and max bounds for each territory and (b) minimise the travel cost between each element and its assigned territory centre (i.e. minimise territory size or dispersion). It's not guaranteed to keep 100% contiguous territories (yet) as it doesn't model adjacency / contiguity constraints directly however if the input problem isn't too tightly constrained on the min/max quantity bounds than objective (b) will usually ensure territories are contiguous (or very near to contiguous).

Our algorithm isn't guaranteed to be optimal either but a lot of work went into the search heuristics, so it will get reasonably good solutions. You could either study the code on github and maybe use some of the search heuristics as a basis for improving your own search heuristics, or just use it directly?

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.