# 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()