# Practical 9 Code Solution

``````"""
Understanding GIS: Practical 9
@author jonnyhuck

Resolve the North-South divide using a Genetic Algorithm

References:
https://projectionwizard.org/#

New Topics:
Lambda Functions
"""

from sys import exit
from shapely.ops import split
from fiona.errors import DriverError
from matplotlib.patches import Patch
from numpy.random import randint, uniform
from matplotlib_scalebar.scalebar import ScaleBar
from shapely.geometry import LineString, MultiPolygon
from matplotlib.pyplot import subplots, savefig, Line2D

def group_polygons(polys, individual):
"""
Convert list of polygons into two multipolygons, grouped above and below the cutline.
"""
# loop through all polygons
top = []
bottom = []
for poly in polys.geoms:

# if the max y value, is it above the cutline, then it is part of the top 'half'
if poly.bounds[3] > individual['y']:
top.append(poly)

# otherwise, it is part of the bottom 'half'
else:
bottom.append(poly)

# return list of two multipolygons
return [MultiPolygon(top), MultiPolygon(bottom)]

def get_fitness(population):
"""
Calculate the fitness of a position on the y-axis by splitting the polygon at
the desired position on the y-aixs and comparing the area of the two resulting
polygons.
"""
# calculate fitness of each individual
for individual in population:

# split the polygon into two along the cutline
polys = group_polygons(split(gb, LineString([(gb.bounds[0], individual['y']), (gb.bounds[2], individual['y'])])), individual)

# overwrite polys with a sorted version of itself
polys = sorted(polys, key=lambda poly: poly.area)

# update fitness as ratio of smallest to largest
individual['fitness'] = polys[0].area / polys[1].area

# return the resulting array
return population

def crossover(parents, offspring_size):
"""
* Single point crossover function for decimal numbers
"""
# loop enough times to make the offspring size
offspring = []
for i in range(offspring_size):

# get binary representations of parents' y values
parent_1 = list(bin(parents[i % len(parents)]['y'])[2:])
parent_2 = list(bin(parents[(i+1) % len(parents)]['y'])[2:])

# swap some random chromasomes (bits) in the binary strings
for r in randint(len(parent_1)-1, size=len(parent_1) // 2):
parent_1[r] = parent_2[r]

# convert back to number and store in a dictionary
offspring.append({'y': int("".join(parent_1), 2), 'fitness': None})

# return the next generation
return offspring

def mutation(population, mutation_probability, max_mutation):
"""
* Mutate a value by +/- max_mutation
"""
# mutation changes a single gene in each offspring randomly.
for i in range(len(population)):

# does this child want to mutate?
if (uniform() < mutation_probability):

# apply the random value as a mutation to the child
population[i]['y'] += int(uniform(-max_mutation, max_mutation))

# return the resulting offspring
return population

# define CRS
uk_ea = '+proj=tcea +lon_0=-3 +datum=WGS84 +units=m +no_defs'

# Load the shapefile and project
try:

# if the file does not exist, warn and exit
except DriverError:
print("Warning, invalid filepath. Exiting.")
exit()

# extract the geometry of the UK
uk = world[(world.ISO_A3 == 'GBR')]['geometry'].iloc[0]

# extract the largest polygon (mainland Great Britain)
gb = sorted(uk.geoms, key=lambda poly: poly.area, reverse=True)[0]

# settings
pop_size                = 50       # population size
num_parents_mating      = 10       # mating pool size (how many of the pop get to breed)
threshold               = 0.999    # the desired precision of the result
mutation_probability    = 0.1      # probability of a child mutating
max_mutation            = 100000.0 # 10km max mutation

# create the initial population (array of dictionaries) then calculate the fitness for each individual
population = [{'y': int(y), 'fitness': None} for y in uniform(low=gb.bounds[1], high=gb.bounds[3], size=pop_size)]
population = get_fitness(population)

# initialise loop varibles
generation = 0
best_fit = 0
previous_best_fit = None

# loop until we either find a solution to within the threshold, or the solutions stop improving
while best_fit < threshold and best_fit != previous_best_fit:

# select the best parents in the population for mating
parents = sorted(population, key=lambda individual: individual['fitness'], reverse=True)[:num_parents_mating]

# get the next generation, mutate and update fitness values
population = get_fitness(mutation(crossover(parents, pop_size), mutation_probability, max_mutation))

# get the current best individual
best_match = sorted(population, key=lambda individual: individual['fitness'], reverse=True)[0]
previous_best_fit = best_fit
best_fit = best_match['fitness']

# increment generation counter and report current fitness (1 = perfect)
generation += 1
print(f"\tgeneration {generation}: {best_fit}")

# report the best match and output
print(f"Best solution : {best_match['y']} (fitness: {best_match['fitness']} generations: {generation})")

# construct the final answer as a linestring
cutline = LineString([(gb.bounds[0], best_match['y']), (gb.bounds[2], best_match['y'])])

# split GB to get the two polygons
polys = group_polygons(split(gb, cutline), best_match)

# setup figure for output
fig, my_ax = subplots(1, 1, figsize=(16, 10))
my_ax.axis('off')
my_ax.set_title("The North-South Divide")

GeoSeries(polys[0], crs=uk_ea).plot(
ax = my_ax,
color = '#ccebc5',
edgecolor = '#4daf4a',
linewidth = 0.3
)
GeoSeries(polys[1], crs=uk_ea).plot(
ax = my_ax,
color = '#b3cde3',
edgecolor = '#377eb8',
linewidth = 0.3
)
GeoSeries(cutline, crs=uk_ea).plot(
ax = my_ax,
color = '#e41a1c',
linewidth = 1
)

# manually draw a legend
my_ax.legend([
Patch(facecolor='#ccebc5', edgecolor='#4daf4a', label='The North'),
Patch(facecolor='#b3cde3', edgecolor='#377eb8', label='The South'),
Line2D([0], [0], color='#e41a1c', lw=1)],
['The North', 'The South', 'Cutline'], loc='lower right')