# program to create Dorling circular cartograms
# written by Zachary Forest Johnson, January 2008
# http://indiemaps.com/blog
# zach.f.johnson@gmail.com
# based on Daniel Dorling's C code
# code is NOT very flexible (currently only works with the provided input file for British counties)
# but this should change very soon
# requires an input file in the following format
# 6713130 530646 180164 6 0 148618 29 34509 22 69378 43 83111 26 51429 11 17013
# Population X-centre Y-centre Number_of_neighbors [neighbour_number common_boundary_length]...(repeated)
# this file works with a file called 'county.in' which you hopefully have
# soon I will write a script to generate the necessary information directly from a shapefile
# stay tuned to indiemaps.com/blog for updates
import math
bodies = 64
friction = 0.25
ratio = 0.4
screen_scale = 0.001
tree = {}
def add_point(pointer, axis):
global end_pointer, tree
if tree[pointer]['id'] == 0:
tree[pointer]['id'] = body
tree[pointer]['left'] = 0
tree[pointer]['right'] = 0
tree[pointer]['xpos'] = x[body]
tree[pointer]['ypos'] = y[body]
else:
if axis == 1:
if x[body] >= tree[pointer]['xpos']:
if tree[pointer]['left'] == 0:
end_pointer += 1
tree[pointer]['left'] = end_pointer
add_point(tree[pointer]['left'], 3-axis)
else:
if tree[pointer]['right'] == 0:
end_pointer += 1
tree[pointer]['right'] = end_pointer
add_point(tree[pointer]['right'], 3-axis)
else:
if y[body] >= tree[pointer]['ypos']:
if tree[pointer]['left'] == 0:
end_pointer += 1
tree[pointer]['left'] = end_pointer
add_point(tree[pointer]['left'], 3-axis)
else:
if tree[pointer]['right'] == 0:
end_pointer += 1
tree[pointer]['right'] = end_pointer
add_point(tree[pointer]['right'], 3-axis)
def get_point(pointer, axis):
global number, mylist
if pointer > 0:
if tree[pointer]['id'] > 0:
if axis == 1:
if x[body]-distance < tree[pointer]['xpos']:
get_point(tree[pointer]['right'],3-axis)
if x[body]+distance >= tree[pointer]['xpos']:
get_point(tree[pointer]['left'],3-axis)
if axis == 2:
if y[body]-distance < tree[pointer]['ypos']:
get_point(tree[pointer]['right'],3-axis)
if y[body]+distance >= tree[pointer]['ypos']:
get_point(tree[pointer]['left'],3-axis)
if (x[body]-distance < tree[pointer]['xpos'] and
x[body]+distance >= tree[pointer]['xpos']):
if (y[body]-distance < tree[pointer]['ypos'] and
y[body]+distance >= tree[pointer]['ypos']):
mylist[number] = tree[pointer]['id']
number += 1
infile = open('county2.in','r')
outfile = open('county.out','w')
#set number of iterations to do etc
itters = 100
done = 0
global body, number, distance, x, y
#read in the data
t_dist = 0.0
t_radius = 0.0
people = {}
x = {}
y = {}
nbours = {}
border = {}
perimeter = {}
nbour = {}
aList = infile.readlines()
monkey = 0
for body in range(bodies):
line = aList[body].split(' ')
people[body] = int(line[0])
if people[body]==78153: print 'coot', body, people[body]
x[body] = int(line[1])
y[body] = int(line[2])
nbours[body] = int(line[3])
perimeter[body] = 0.0
border[body] = {}
nbour[body] = {}
for nb in range(nbours[body]):
nbour[body][nb] = int(line[4+(nb)*2])
boundary = int(line[5+(nb)*2])
border[body][nb] = boundary
if nbour[body][nb] != -999:
perimeter[body] += border[body][nb]
if nbour[body][nb] > 0:
if nbour[body][nb] < body:
monkey += 1
xd = (x[body] - x[nbour[body][nb]])
yd = (y[body] - y[nbour[body][nb]])
t_dist += math.sqrt(xd*xd+yd*yd)
t_radius += math.sqrt(people[body]/math.pi) + math.sqrt(people[nbour[body][nb]]/math.pi)
print monkey
print t_dist
print t_radius
scale = t_dist / t_radius
# print scale
widest = 0.0
radius = {}
xvector = {}
yvector = {}
mylist = {}
for body in range(bodies):
radius[body] = scale * math.sqrt(people[body]/math.pi)
if radius[body] > widest:
widest = radius[body]
xvector[body] = 0.0
yvector[body] = 0.0
print('Scaling by %lf widest is %lf\n' % (scale,widest))
#start the big loop creating the tree each iter
for itter in range(itters):
for body in range(1,bodies):
tree[body] = {}
tree[body]['id'] = 0
end_pointer = 1
for body in range(bodies):
add_point(1,1)
displacement = 0.0
#loop of independent body movements
for body in range(bodies):
#get of neighbours within into
number = 0
distance = widest + radius[body]
get_point(1,1)
#initialize a few vectors
xrepel = yrepel = 0.0
xattract = yattract = 0.0
closest = widest
#work out repelling force of overlapping neighbours
if number > 0:
for nb in range(number):
other = mylist[nb]
if other != body:
xd = x[other]-x[body]
yd = y[other]-y[body]
dist = math.sqrt(xd * xd + yd * yd)
if dist < closest:
closest = dist
overlap = radius[body] + radius[other] - dist
if overlap > 0.0:
if dist > 1.0:
xrepel = xrepel - overlap*(x[other]-x[body])/dist
yrepel = yrepel - overlap*(y[other]-y[body])/dist
#work out forces of attraction between neighbours
for nb in range(nbours[body]):
other = nbour[body][nb]
if other != 0:
xd = (x[body]-x[other])
yd = (y[body]-y[other])
dist = math.sqrt(xd * xd + yd * yd)
overlap = dist - radius[body] - radius[other]
if overlap > 0.0:
overlap = overlap * border[body][nb]/perimeter[body]
xattract = xattract + overlap*(x[other]-x[body])/dist
yattract = yattract + overlap*(y[other]-y[body])/dist
#now work out the combined effect of attraction and repulsion */
atrdst = math.sqrt(xattract * xattract + yattract * yattract)
repdst = math.sqrt(xrepel * xrepel + yrepel * yrepel)
if repdst > closest:
xrepel = closest * xrepel / (repdst + 1.0)
yrepel = closest * yrepel / (repdst + 1.0)
repdst = closest
if repdst > 0.0:
xtotal = (1.0-ratio) * xrepel + ratio*(repdst*xattract/(atrdst + 1.0))
ytotal = (1.0-ratio) * yrepel + ratio*(repdst*yattract/(atrdst + 1.0))
else:
if atrdst > closest:
xattract = closest * xattract/(atrdst+1.0)
yattract = closest * yattract/(atrdst+1.0)
xtotal = xattract
ytotal = yattract
xvector[body] = friction * (xvector[body] + xtotal)
yvector[body] = friction * (yvector[body] + ytotal)
displacement += math.sqrt(xtotal * xtotal + ytotal * ytotal)
#update the positions
for body in range(bodies):
x[body] += xvector[body] + 0.5
y[body] += yvector[body] + 0.5
done += 1
displacement = displacement / bodies
if done % 10 == 1:
print('displacement is %f after %i iterations\n'% (displacement, done))
#write out the new file
print('writing out the new file \n')
outfile.write('%i %i\n'% (itters, done))
outfile.write('%i %f\n'% (bodies,widest))
for body in range(bodies):
#print('%f %f %i %f %i %i %i %f '% (xvector[body], yvector[body], people[body],radius[body], x[body], y[body],nbours[body], perimeter[body]))
outfile.write('%f %f %i %f %i %i %i %f '% (xvector[body], yvector[body], people[body],radius[body], x[body], y[body],nbours[body], perimeter[body]))
for nb in range(nbours[body]):
outfile.write('%i %f '% (nbour[body][nb], border[body][nb]))
outfile.write('\n')
outfile.close()