# Module `nemo.polygons`

Support code for the 43 polygons of the AEMO study.

Expand source code
``````# Copyright (C) 2014, 2015, 2016 The University of New South Wales
#
# This file is free software; you can redistribute it and/or modify it
# the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.

"""Support code for the 43 polygons of the AEMO study."""

import json
from math import atan2, sin, cos, radians, sqrt
import numpy as np

from dijkstra import dijkstra
from nemo import regions

# The fraction of a region's load in each polygon.
regions.nsw.polygons = {21: 0, 22: 0, 23: 0, 24: .05, 28: 0, 29: 0,
30: 0, 31: .8, 33: 0, 34: 0, 35: .05, 36: .1}
regions.qld.polygons = {1: .04, 2: 0, 3: 0, 4: .11, 5: 0, 6: 0,
7: .27, 8: 0, 9: 0, 10: 0, 11: .02, 14: 0, 15: 0,
16: .14, 17: .42}
regions.sa.polygons = {12: 0, 13: 0, 18: 0, 19: 0, 20: 0, 25: 0,
26: 0, 27: .1, 32: .9}
regions.snowy.polygons = {}
regions.tas.polygons = {40: .2, 41: .2, 42: 0, 43: .6}
regions.vic.polygons = {37: .2, 38: .1, 39: .7}

# Ensure all weights sum to one.
for r in regions.All:
if r.polygons:
assert round(sum(r.polygons.values())) == 1

# Useful for testing
wildcard = 31

# Vertices of the closed polygons (nb. must be closed)
_polygons = {
1: ((144.602, -13.838), (145.602, -13.838), (148.447, -18.282), (146.646, -19.041), (145.107, -18.792), (143.459, -18.480),
(144.602, -13.838)),
2: ((140.949, -18.099), (143.459, -18.480), (143.701, -20.715), (140.949, -20.612), (140.949, -18.099)),
3: ((143.459, -18.480), (145.107, -18.792), (146.646, -20.797), (143.701, -20.715), (143.459, -18.480)),
4: ((148.447, -18.282), (150.754, -21.545), (147.810, -22.513), (146.646, -20.797), (145.107, -18.792), (146.646, -19.041),
(148.447, -18.282)),
5: ((140.949, -20.612), (143.701, -20.715), (143.987, -23.665), (140.949, -23.665), (140.949, -20.612)),
6: ((143.701, -20.715), (146.646, -20.797), (147.810, -22.513), (148.887, -23.665), (146.382, -23.665), (143.987, -23.665),
(143.701, -20.715)),
7: ((150.754, -21.545), (152.974, -24.137), (149.639, -24.817), (148.887, -23.665), (147.810, -22.513),
(150.754, -21.545)),
8: ((140.949, -23.665), (143.987, -23.665), (144.141, -25.958), (140.999, -25.996), (140.949, -23.665)),
9: ((143.987, -23.665), (146.382, -23.665), (147.458, -25.958), (144.141, -25.958), (143.987, -23.665)),
10: ((146.382, -23.665), (148.887, -23.665), (149.639, -24.817), (150.529, -25.958), (147.458, -25.958),
(146.382, -23.665)),
11: ((152.974, -24.137), (154.457, -25.958), (150.529, -25.958), (149.639, -24.817), (152.974, -24.137)),
12: ((135.183, -25.999), (137.933, -25.997), (137.933, -29.075), (135.187, -29.075), (135.183, -25.999)),
13: ((137.933, -25.997), (140.999, -25.996), (140.999, -28.999), (137.933, -29.075), (137.933, -25.997)),
14: ((140.999, -25.996), (144.141, -25.958), (144.232, -28.999), (141.001, -28.999), (140.999, -25.996)),
15: ((144.141, -25.958), (147.458, -25.958), (147.550, -28.999), (144.232, -28.999), (144.141, -25.958)),
16: ((147.458, -25.958), (150.529, -25.958), (150.688, -28.999), (147.550, -28.999), (147.458, -25.958)),
17: ((154.457, -25.958), (154.852, -29.075), (150.688, -28.999), (150.529, -25.958), (154.457, -25.958)),
18: ((131.199, -29.075), (135.187, -29.075), (135.187, -31.325), (131.199, -31.325), (131.199, -29.075)),
19: ((135.187, -29.075), (137.933, -29.075), (137.933, -31.325), (135.187, -31.325), (135.187, -29.075)),
20: ((137.933, -29.075), (140.999, -28.999), (141.001, -30.999), (141.001, -31.354), (137.933, -31.325),
(137.933, -29.075)),
21: ((141.001, -28.999), (144.232, -28.999), (144.141, -31.109), (141.001, -30.998), (141.001, -28.999)),
22: ((144.232, -28.999), (147.550, -28.999), (146.843, -31.840), (144.141, -31.766), (144.141, -31.109),
(144.232, -28.999)),
23: ((147.550, -28.999), (150.688, -28.999), (149.359, -31.878), (146.843, -31.840), (147.550, -28.999)),
24: ((154.852, -29.075), (153.798, -32.639), (149.359, -31.878), (150.688, -28.999), (154.852, -29.075)),
25: ((131.199, -31.325), (135.187, -31.325), (133.638, -34.053), (131.199, -32.658), (131.199, -31.325)),
26: ((135.187, -31.325), (137.933, -31.325), (137.900, -33.852), (137.856, -36.985), (135.700, -36), (133.638, -34.053),
(135.187, -31.325)),
27: ((137.933, -31.326), (141.001, -31.354), (141.002, -33.311), (141.003, -33.982), (140.964, -33.981),
(137.900, -33.852), (137.933, -31.326)),
28: ((141.001, -30.998), (144.141, -31.109), (144.141, -31.766), (143.954, -33.303), (141.002, -33.311),
(141.001, -31.354), (141.001, -30.998)),
29: ((144.141, -31.766), (146.843, -31.840), (146.250, -34.053), (143.943, -34.016), (143.954, -33.303),
(144.141, -31.766)),
30: ((146.843, -31.840), (149.359, -31.878), (148.315, -34.107), (146.250, -34.053), (146.843, -31.840)),
31: ((153.798, -32.639), (152.260, -34.724), (148.315, -34.107), (149.359, -31.878), (153.798, -32.639)),
32: ((137.900, -33.852), (140.964, -33.981), (140.964, -33.990), (140.963, -33.990), (140.962, -34.110),
(140.966, -35.237), (140.966, -35.389), (140.964, -35.749), (140.974, -37.359), (140.971, -37.791),
(140.966, -38.056), (140.966, -38.568), (137.856, -36.985), (137.900, -33.852)),
33: ((141.002, -33.311), (143.954, -33.303), (143.943, -34.016), (143.811, -35.299), (140.966, -35.237),
(140.962, -34.108),
(140.963, -33.990), (140.964, -33.990), (140.964, -33.981), (141.003, -33.982), (141.002, -33.311)),
34: ((143.943, -34.016), (146.250, -34.053), (145.698, -35.989), (143.811, -35.299), (143.943, -34.016)),
35: ((146.250, -34.053), (148.315, -34.107), (147.118, -36.510), (145.698, -35.989), (146.250, -34.053)),
36: ((152.260, -34.724), (150.590, -37.831), (147.118, -36.510), (148.315, -34.107), (152.260, -34.724)),
37: ((140.966, -35.237), (143.811, -35.299), (143.483, -39.249), (140.968, -38.568), (140.966, -38.056),
(140.971, -37.791), (140.975, -37.359), (140.963, -35.742), (140.966, -35.389), (140.966, -35.237)),
38: ((146.426, -40.581), (145.698, -35.989), (147.118, -36.510), (150.590, -37.831), (146.426, -40.581)),
39: ((145.698, -35.989), (146.426, -40.581), (143.481, -39.249), (143.811, -35.299), (145.698, -35.989)),
40: ((143.483, -39.249), (146.426, -40.581), (146.426, -42.033), (144.097, -42.033), (143.483, -39.249)),
41: ((146.426, -40.581), (149.052, -38.849), (149.052, -42.033), (146.426, -42.033), (146.426, -40.581)),
42: ((144.097, -42.033), (146.426, -42.033), (146.426, -44), (144.097, -43.747), (144.097, -42.033)),
43: ((146.426, -42.033), (149.052, -42.033), (149.052, -43.747), (146.426, -44), (146.426, -42.033)),
}

numpolygons = len(_polygons)

# Table mapping polygon number to region.
_region_table = [None] * (numpolygons + 1)
for rgn in [regions.nsw, regions.qld, regions.sa, regions.tas, regions.vic]:
for poly in rgn.polygons:
_region_table[poly] = rgn

def region(polygon):
"""Return the region a polygon resides in.

>>> region(1)
QLD1
>>> region(40)
TAS1
"""
return _region_table[polygon]

wind_limit = [None, 80.3, 0, 36.9, 6.5, 15.6, 1.5, 6.9, 2.6, 0, 4.1,
1.5, 2.1, 0.9, 30.3, 0, 0, 40.5, 0.2, 0, 49.1, 2.3, 0,
1.7, 116.3, 3.3, 71.9, 128.3, 11.7, 0.5, 0.6, 52.5,
20.0, 0, 0, 0.9, 101.0, 9.15, 10.2, 15.6, 11.4, 14.1,
0.5, 29.1]

pv_limit = [None, 133, 1072, 217, 266, 1343, 1424, 287, 1020, 657,
175, 47, 488, 749, 1338, 1497, 1093, 243, 558, 647, 639,
921, 1310, 1182, 125, 81, 493, 689, 937, 736, 522, 31,
527, 535, 618, 339, 26, 670, 78, 347, 13, 21, 0.21, 5]

cst_limit = [None, 102, 822, 166, 204, 1030, 1092, 220, 782, 504, 134,
36, 374, 574, 1026, 1148, 838, 186, 428, 496, 490, 706,
1004, 906, 96, 62, 378, 528, 718, 564, 400, 24, 404, 410,
474, 260, 20, 514, 60, 266, 10, 16, 0.16, 4]

def dumps():
"""Dump the polygon data in GeoJSON format.

>>> assert len(dumps()) > 0
"""
polys = []
for i in range(1, numpolygons):
geometry = {'type': 'Polygon', 'coordinates': [_polygons[i]]}
properties = {'Polygon #': i, 'CST limit (GW)': cst_limit[i], 'PV limit (GW)': pv_limit[i],
'Wind limit (GW)': wind_limit[i]}
feature = {'type': 'Feature', 'geometry': geometry, 'properties': properties}
polys.append(feature)
top = {'type': 'FeatureCollection', 'features': polys}
return json.dumps(top)

def _centroid(vertices):
"""Find the centroid of a polygon."""

# Ensure the polygon is closed
assert vertices[0] == vertices[-1]
thesum = 0
vsum = (0, 0)
for i in range(len(vertices) - 1):
v1 = vertices[i]
v2 = vertices[i + 1]
cross = v1[0] * v2[1] - v1[1] * v2[0]
thesum += cross
vsum = (((v1[0] + v2[0]) * cross) + vsum[0], ((v1[1] + v2[1]) * cross) + vsum[1])
z = 1. / (3. * thesum)
return (vsum[0] * z, vsum[1] * z)

def path(poly1, poly2):
"""
Return a path from polygon 1 to polygon 2.

>>> path(1, 30)
[(1, 4), (4, 10), (10, 16), (16, 23), (23, 30)]
>>> path(23, 43)
[(23, 30), (30, 35), (35, 38), (38, 41), (41, 43)]
"""
return connections[(poly1, poly2)]

def subset(p, polysuperset):
"""
Are all polygons in path P present in superset?

>>> subset([(1,2), (2,3)], [1,2,3])
True
>>> subset([(1,4), (4,3)], [1,2,3])
False
"""
# Flatten the list of pairs into one long list.
polylist = [i for sub in p for i in sub]
# Now for a simple set operation.
return set(polylist) <= set(polysuperset)

def direct_p(poly1, poly2):
"""
Return True if polygon 1 and polygon 2 are directly connected.

>>> direct_p(1, 2)
True
>>> direct_p(1, 40)
False
"""
return len(path(poly1, poly2)) <= 1

def dist(poly1, poly2):
"""Return the distance between two polygon centroids.

>>> dist(1,1)
0
>>> dist(1,43)
2910
>>> dist(1,43) == distances[1,43]
True
"""
# Code adapted from Chris Veness
point1 = centroids[poly1]
point2 = centroids[poly2]
a = sin(dlat / 2) ** 2 + \
sin(dlon / 2) ** 2 * cos(lat1) * cos(lat2)
c = 2 * atan2(sqrt(a), sqrt(1 - a))

def pathlen(p):
"""Return the total length of a path P.

>>> pathlen([])
0
>>> x = dist(1,4)
>>> assert pathlen([(1, 4)]) == x
>>> y = dist(4,7)
>>> assert pathlen([(1, 4), (4, 7)]) == x + y
"""
total = 0
for poly1, poly2 in p:
total += dist(poly1, poly2)

centroids = {}
for _i, _vertices in _polygons.items():
_lon, _lat = _centroid(_vertices)
centroids[_i] = (_lat, _lon)

# A proposed transmission network.

net = {1: {2: dist(1, 2), 3: dist(1, 3), 4: dist(1, 4)},
2: {1: dist(2, 1), 3: dist(2, 3), 5: dist(2, 5)},
3: {1: dist(3, 1), 2: dist(3, 2), 4: dist(3, 4), 6: dist(3, 6)},
4: {1: dist(4, 1), 3: dist(4, 3), 6: dist(4, 6), 7: dist(4, 7), 10: dist(4, 10)},
5: {2: dist(5, 2), 6: dist(5, 6), 8: dist(5, 8)},
6: {3: dist(6, 3), 4: dist(6, 4), 5: dist(6, 5), 9: dist(6, 9)},
7: {4: dist(7, 4), 10: dist(7, 10), 11: dist(7, 11), 16: dist(7, 16)},
8: {5: dist(8, 5), 9: dist(8, 9), 14: dist(8, 14)},
9: {6: dist(9, 6), 8: dist(9, 8), 10: dist(9, 10), 15: dist(9, 15)},
10: {4: dist(10, 4), 7: dist(10, 7), 9: dist(10, 9), 16: dist(10, 16)},
11: {7: dist(7, 11), 16: dist(16, 11), 17: dist(11, 17)},
12: {13: dist(12, 13), 19: dist(12, 19)},
13: {12: dist(13, 12), 14: dist(13, 14), 20: dist(13, 20)},
14: {8: dist(14, 8), 13: dist(14, 13), 15: dist(14, 15), 21: dist(14, 21)},
15: {9: dist(15, 9), 14: dist(15, 14), 16: dist(15, 16), 22: dist(15, 22)},
16: {7: dist(16, 7), 10: dist(16, 10), 11: dist(16, 11), 15: dist(16, 15), 17: dist(16, 17), 23: dist(16, 23),
24: dist(16, 24)},
17: {11: dist(17, 11), 16: dist(17, 16), 24: dist(17, 24)},
18: {19: dist(18, 19), 25: dist(18, 25)},
19: {12: dist(19, 12), 18: dist(19, 18), 20: dist(19, 20), 26: dist(19, 26)},
20: {13: dist(20, 13), 19: dist(20, 19), 21: dist(20, 21), 27: dist(20, 27)},
21: {14: dist(21, 14), 20: dist(21, 20), 22: dist(21, 22), 29: dist(21, 29)},
22: {15: dist(22, 15), 21: dist(22, 21), 23: dist(22, 23), 29: dist(22, 29)},
23: {16: dist(23, 16), 22: dist(23, 22), 24: dist(23, 24), 30: dist(23, 30)},
24: {16: dist(24, 16), 17: dist(24, 17), 23: dist(24, 23), 31: dist(23, 31)},
25: {18: dist(25, 18), 26: dist(25, 26)},
26: {19: dist(26, 19), 25: dist(26, 25), 27: dist(26, 27)},
27: {20: dist(27, 20), 26: dist(27, 26), 28: dist(27, 28), 32: dist(27, 32), 33: dist(27, 33)},
28: {21: dist(28, 21), 27: dist(28, 27), 29: dist(28, 29), 33: dist(28, 33)},
29: {22: dist(29, 22), 28: dist(29, 28), 30: dist(29, 30), 34: dist(29, 34)},
30: {23: dist(30, 23), 29: dist(30, 29), 31: dist(30, 31), 35: dist(30, 35)},
31: {24: dist(31, 24), 30: dist(31, 30), 36: dist(31, 36)},
32: {27: dist(32, 27), 33: dist(32, 33), 37: dist(32, 37)},
33: {27: dist(33, 27), 28: dist(33, 28), 32: dist(33, 32), 34: dist(33, 34), 37: dist(33, 37), 39: dist(33, 39)},
34: {29: dist(34, 29), 33: dist(34, 33), 35: dist(34, 35), 39: dist(34, 39)},
35: {30: dist(35, 30), 34: dist(35, 34), 36: dist(35, 36), 38: dist(35, 38), 39: dist(35, 39)},
36: {31: dist(36, 31), 35: dist(36, 35), 38: dist(36, 38)},
37: {32: dist(37, 32), 33: dist(37, 33), 39: dist(37, 39)},
38: {35: dist(38, 35), 36: dist(38, 36), 39: dist(38, 39), 41: dist(38, 41)},
39: {33: dist(39, 33), 34: dist(39, 34), 35: dist(39, 35), 37: dist(39, 37), 38: dist(39, 38), 40: dist(39, 40)},
40: {39: dist(40, 39), 41: dist(40, 41), 42: dist(40, 42)},
41: {38: dist(41, 38), 40: dist(41, 40), 43: dist(41, 43)},
42: {40: dist(42, 40), 43: dist(42, 43)},
43: {41: dist(43, 41), 42: dist(43, 42)}}

distances = np.zeros((numpolygons + 1, numpolygons + 1))
# mark row 0 and column 0 as unused (there is no polygon #0)
distances[0] = np.nan
distances[::, 0] = np.nan
rows, cols = numpolygons + 1, numpolygons + 1
for p1 in range(1, rows):
for p2 in range(1, cols):
distances[p1, p2] = dist(p1, p2)

existing_net = np.zeros((numpolygons + 1, numpolygons + 1))
# mark row 0 and column 0 as unused (there is no polygon #0)
existing_net[0] = np.nan
existing_net[::, 0] = np.nan

for (p1, p2, limit) in \
[(7, 4, 1100), (7, 16, 1000), (7, 11, 1100), (11, 17, 1100),
(16, 7, 400), (17, 11, 200), (11, 7, 200), (16, 17, 3500),
(16, 24, 1200), (17, 24, 250), (24, 16, 500), (24, 17, 100),
(24, 31, 900), (31, 24, 1100), (31, 36, 2000), (36, 31, 2000),
(36, 35, 1500), (35, 36, 2800), (33, 34, 100), (34, 33, 100),
(32, 27, 100), (27, 32, 550), (32, 33, 150), (33, 32, 230),
(32, 37, 900), (37, 32, 900), (37, 39, 900), (39, 37, 900),
(33, 39, 500), (39, 33, 1300), (34, 39, 1000), (39, 34, 1300),
(39, 38, 2500), (38, 39, 6400), (38, 41, 450), (41, 38, 600)]:
assert p1 in list(net[p2].keys()), (p2, p1)
assert p2 in list(net[p1].keys()), (p1, p2)
existing_net[p1, p2] = limit

connections = {}
for dest in range(1, numpolygons + 1):
for src in range(1, numpolygons + 1):
shortest = list(dijkstra.shortestPath(net, src, dest))
pairs = []
for _i in range(len(shortest) - 1):
pairs.append((shortest[_i], shortest[_i + 1]))
connections[(src, dest)] = pairs``````

## Functions

``` def direct_p(poly1, poly2) ```

Return True if polygon 1 and polygon 2 are directly connected.

``````>>> direct_p(1, 2)
True
>>> direct_p(1, 40)
False
``````
Expand source code
``````def direct_p(poly1, poly2):
"""
Return True if polygon 1 and polygon 2 are directly connected.

>>> direct_p(1, 2)
True
>>> direct_p(1, 40)
False
"""
return len(path(poly1, poly2)) <= 1``````
``` def dist(poly1, poly2) ```

Return the distance between two polygon centroids.

``````>>> dist(1,1)
0
>>> dist(1,43)
2910
>>> dist(1,43) == distances[1,43]
True
``````
Expand source code
``````def dist(poly1, poly2):
"""Return the distance between two polygon centroids.

>>> dist(1,1)
0
>>> dist(1,43)
2910
>>> dist(1,43) == distances[1,43]
True
"""
# Code adapted from Chris Veness
point1 = centroids[poly1]
point2 = centroids[poly2]
a = sin(dlat / 2) ** 2 + \
sin(dlon / 2) ** 2 * cos(lat1) * cos(lat2)
c = 2 * atan2(sqrt(a), sqrt(1 - a))
``` def dumps() ```

Dump the polygon data in GeoJSON format.

``````>>> assert len(dumps()) > 0
``````
Expand source code
``````def dumps():
"""Dump the polygon data in GeoJSON format.

>>> assert len(dumps()) > 0
"""
polys = []
for i in range(1, numpolygons):
geometry = {'type': 'Polygon', 'coordinates': [_polygons[i]]}
properties = {'Polygon #': i, 'CST limit (GW)': cst_limit[i], 'PV limit (GW)': pv_limit[i],
'Wind limit (GW)': wind_limit[i]}
feature = {'type': 'Feature', 'geometry': geometry, 'properties': properties}
polys.append(feature)
top = {'type': 'FeatureCollection', 'features': polys}
return json.dumps(top)``````
``` def path(poly1, poly2) ```

Return a path from polygon 1 to polygon 2.

``````>>> path(1, 30)
[(1, 4), (4, 10), (10, 16), (16, 23), (23, 30)]
>>> path(23, 43)
[(23, 30), (30, 35), (35, 38), (38, 41), (41, 43)]
``````
Expand source code
``````def path(poly1, poly2):
"""
Return a path from polygon 1 to polygon 2.

>>> path(1, 30)
[(1, 4), (4, 10), (10, 16), (16, 23), (23, 30)]
>>> path(23, 43)
[(23, 30), (30, 35), (35, 38), (38, 41), (41, 43)]
"""
return connections[(poly1, poly2)]``````
``` def pathlen(p) ```

Return the total length of a path P.

``````>>> pathlen([])
0
>>> x = dist(1,4)
>>> assert pathlen([(1, 4)]) == x
>>> y = dist(4,7)
>>> assert pathlen([(1, 4), (4, 7)]) == x + y
``````
Expand source code
``````def pathlen(p):
"""Return the total length of a path P.

>>> pathlen([])
0
>>> x = dist(1,4)
>>> assert pathlen([(1, 4)]) == x
>>> y = dist(4,7)
>>> assert pathlen([(1, 4), (4, 7)]) == x + y
"""
total = 0
for poly1, poly2 in p:
total += dist(poly1, poly2)
``` def region(polygon) ```

Return the region a polygon resides in.

``````>>> region(1)
QLD1
>>> region(40)
TAS1
``````
Expand source code
``````def region(polygon):
"""Return the region a polygon resides in.

>>> region(1)
QLD1
>>> region(40)
TAS1
"""
return _region_table[polygon]``````
``` def subset(p, polysuperset) ```

Are all polygons in path P present in superset?

``````>>> subset([(1,2), (2,3)], [1,2,3])
True
>>> subset([(1,4), (4,3)], [1,2,3])
False
``````
Expand source code
``````def subset(p, polysuperset):
"""
Are all polygons in path P present in superset?

>>> subset([(1,2), (2,3)], [1,2,3])
True
>>> subset([(1,4), (4,3)], [1,2,3])
False
"""
# Flatten the list of pairs into one long list.
polylist = [i for sub in p for i in sub]
# Now for a simple set operation.
return set(polylist) <= set(polysuperset)``````