# 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
# under the terms of the GNU General Public License as published by
# 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."""

from math import atan2, cos, radians, sin, sqrt

import numpy as np

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)

# Dictionary mapping polygon number to region.
_region_table = {}
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]

# These build limits come from the ROAM Consulting report on wind and
# solar modelling for the AEMO 100% Renewables project. See
# nemo.ozlabs.org for a link to this report.

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]

# Only these four polygons have been chosen for off-shore wind farm siting.
offshore_wind_limit = {31: 10, 36: 10, 38: 10, 40: 10}

def _centroid(vertices):
"""Find the centroid of a polygon."""
# pylint: disable=invalid-name
# 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 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
# pylint: disable=invalid-name
radius = 6371  # km
point1 = centroids[poly1]
point2 = centroids[poly2]
dlat = radians(point1[0] - point2[0])
dlon = radians(point1[1] - point2[1])
lat1 = radians(point1[0])
lat2 = radians(point2[0])
a = sin(dlat / 2) ** 2 + \
sin(dlon / 2) ** 2 * cos(lat1) * cos(lat2)
c = 2 * atan2(sqrt(a), sqrt(1 - a))
return int(radius * c)

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``````

## Functions

``` 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
# pylint: disable=invalid-name
radius = 6371  # km
point1 = centroids[poly1]
point2 = centroids[poly2]
dlat = radians(point1[0] - point2[0])
dlon = radians(point1[1] - point2[1])
lat1 = radians(point1[0])
lat2 = radians(point2[0])
a = sin(dlat / 2) ** 2 + \
sin(dlon / 2) ** 2 * cos(lat1) * cos(lat2)
c = 2 * atan2(sqrt(a), sqrt(1 - a))
return int(radius * c)``````
``` 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]``````