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

Your Answer

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

Browse other questions tagged or ask your own question.