class Towns:
    def get_config(self):
        config = configparser.RawConfigParser()
        if len(sys.argv) > 1:
            config_file = sys.argv[1]
        else:
            config_file = getModelFile('SIREN.ini')
        config.read(config_file)
        try:
            self.base_year = config.get('Base', 'year')
        except:
            self.base_year = '2012'
        parents = []
        try:
            parents = getParents(config.items('Parents'))
        except:
            pass
        try:
            self.bom_file = config.get('Files', 'bom')
            for key, value in parents:
                self.bom_file = self.bom_file.replace(key, value)
            self.bom_file = self.bom_file.replace('$USER$', getUser())
            self.bom_file = self.bom_file.replace('$YEAR$', self.base_year)
        except:
            self.bom_file = ''
        try:
            self.town_file = config.get('Files', 'towns')
            for key, value in parents:
                self.town_file = self.town_file.replace(key, value)
            self.town_file = self.town_file.replace('$USER$', getUser())
            self.town_file = self.town_file.replace('$YEAR$', self.base_year)
        except:
            self.town_file = ''

    def haversine(self, lat1, lon1, lat2, lon2):
        """
        Calculate the great circle distance between two points
        on the earth (specified in decimal degrees)
        """
     # convert decimal degrees to radians
        lon1, lat1, lon2, lat2 = list(map(radians, [lon1, lat1, lon2, lat2]))

     # haversine formula
        dlon = lon2 - lon1
        dlat = lat2 - lat1
        a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
        c = 2 * asin(sqrt(a))

     # 6367 km is the radius of the Earth
        km = 6367 * c
        return km

#   might extend capability to restrict towns to network boundary
    def town_in_boundary(self, x, y, poly):
        n = len(poly)
        inside = False
        p1x, p1y = poly[0]
        for i in range(n + 1):
            p2x, p2y = poly[i % n]
            if y > min(p1y, p2y):
                if y <= max(p1y, p2y):
                    if x <= max(p1x, p2x):
                        if p1y != p2y:
                            xints = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x
                        if p1x == p2x or x <= xints:
                            inside = not inside
            p1x, p1y = p2x, p2y
        return inside

    def __init__(self, ul_lat=None, ul_lon=None, lr_lat=None, lr_lon=None, remove_duplicates=True):
        """Initializes the data."""

        def read_town_csv(town_file):
            townfile = open(town_file)
            twns = csv.DictReader(townfile)
            if 'Town' in twns.fieldnames:
                name_field = 'Town'
            elif 'Name' in twns.fieldnames:
                name_field = 'Name'
            else:
                name_field = 'Site name'
            if remove_duplicates:
                for twn in twns:
                    for town in self.towns:
                        if twn[name_field] == town.name:
                            break
                        if len(twn[name_field]) < len(town.name) and town.name.find(' ') > 0:
                            if twn[name_field].title() == town.name[:len(twn[name_field])]:
                                break
                            if town.name[:6].title() == 'North ' and twn[name_field] == town.name[6:]:
                                break
                        if len(twn[name_field]) > len(town.name) and twn[name_field].find(' ') > 0:
                            if twn[name_field].title()[:len(town.name)] == town.name:
                                break
                    else:
                        if ul_lat is None or \
                          (float(twn['Latitude']) >= lr_lat and
                           float(twn['Latitude']) <= ul_lat and
                           float(twn['Longitude']) >= ul_lon and
                           float(twn['Longitude']) <= lr_lon):
                            if 'Lid' in twns.fieldnames and twn['Lid'] != '':
                                self.towns.append(Town(twn[name_field], twn['Latitude'], twn['Longitude'], twn['Lid'],
                                             twn['State'], twn['Country'], twn['Elev'], twn['Zone']))
                            else:
                                self.towns.append(Town(twn[name_field], twn['Latitude'], twn['Longitude']))
            else:
                for twn in twns:
                    if ul_lat is None or \
                      (float(twn['Latitude']) >= lr_lat and
                       float(twn['Latitude']) <= ul_lat and
                       float(twn['Longitude']) >= ul_lon and
                       float(twn['Longitude']) <= lr_lon):
                        if 'Lid' in twns.fieldnames and twn['Lid'] != '':
                            self.towns.append(Town(twn[name_field], twn['Latitude'], twn['Longitude'], twn['Lid'],
                                        twn['State'], twn['Country'], twn['Elev'], twn['Zone']))
                        else:
                            self.towns.append(Town(twn[name_field], twn['Latitude'], twn['Longitude']))
            townfile.close()

        self.get_config()
        self.towns = []
#   Process BOM stations first
        if os.path.exists(self.bom_file):
            if self.bom_file[-4:] == '.csv':
                read_town_csv(self.bom_file)
            else:
                c_closed = -1
                c_id = -1
                c_nme = -1
                c_state = -1
                c_lat = -1
                c_lon = -1
                c_elev = -1
                c_ctry = -1
                workbook = xlrd.open_workbook(self.bom_file)
                worksheet = workbook.sheet_by_index(0)
                num_rows = worksheet.nrows - 1
                num_cols = worksheet.ncols - 1
#     get column names
                curr_col = -1
                while curr_col < num_cols:
                    curr_col += 1
                    if worksheet.cell_value(0, curr_col) == 'Bureau of Meteorology Station Number' or \
                       worksheet.cell_value(0, curr_col) == 'Lid':
                        c_id = curr_col
                    elif worksheet.cell_value(0, curr_col) == 'Station Name' or \
                       worksheet.cell_value(0, curr_col) == 'Site Name' or \
                       worksheet.cell_value(0, curr_col) == 'Town':
                        c_nme = curr_col
                    elif worksheet.cell_value(0, curr_col) == 'Latitude to 4 decimal places, in decimal degrees' or \
                       worksheet.cell_value(0, curr_col) == 'Latitude':
                        c_lat = curr_col
                    elif worksheet.cell_value(0, curr_col) == 'Longitude to 4 decimal places, in decimal degrees' or \
                       worksheet.cell_value(0, curr_col) == 'Longitude':
                        c_lon = curr_col
                    elif worksheet.cell_value(0, curr_col) == 'State':
                        c_state = curr_col
                    elif worksheet.cell_value(0, curr_col) == 'Height of station above mean sea level in metres' or \
                       worksheet.cell_value(0, curr_col) == 'Elev':
                        c_elev = curr_col
                    elif worksheet.cell_value(0, curr_col) == 'Month/Year site closed. (MM/YYYY)':
                        c_closed = curr_col
                    elif worksheet.cell_value(0, curr_col) == 'Country':
                        c_ctry = curr_col
#     WMO (World Meteorological Organisation) Index Number
                curr_row = 0
                while curr_row < num_rows:
                    curr_row += 1
                    if c_closed >= 0:
                        if worksheet.cell_value(curr_row, c_nme) != '' \
                        and worksheet.cell_value(curr_row, c_closed).strip() == '':
                            if ul_lat is None or \
                              (worksheet.cell_value(curr_row, c_lat) >= lr_lat and
                               worksheet.cell_value(curr_row, c_lat) <= ul_lat and
                               worksheet.cell_value(curr_row, c_lon) >= ul_lon and
                               worksheet.cell_value(curr_row, c_lon) <= lr_lon):
                                self.towns.append(Town(worksheet.cell_value(curr_row, c_nme),
                                                       worksheet.cell_value(curr_row, c_lat),
                                                       worksheet.cell_value(curr_row, c_lon),
                                                       worksheet.cell_value(curr_row, c_id),
                                                       worksheet.cell_value(curr_row, c_state),
                                                       elev=worksheet.cell_value(curr_row, c_elev)))
                    else:
                        if worksheet.cell_value(curr_row, c_nme) != '':
                            if ul_lat is None or \
                              (worksheet.cell_value(curr_row, c_lat) >= lr_lat and
                               worksheet.cell_value(curr_row, c_lat) <= ul_lat and
                               worksheet.cell_value(curr_row, c_lon) >= ul_lon and
                               worksheet.cell_value(curr_row, c_lon) <= lr_lon):
                                self.towns.append(Town(worksheet.cell_value(curr_row, c_nme),
                                                       worksheet.cell_value(curr_row, c_lat),
                                                       worksheet.cell_value(curr_row, c_lon)))
                                if c_id >= 0:
                                    self.towns[-1].lid = worksheet.cell_value(curr_row, c_id)
                                if c_state >= 0:
                                    self.towns[-1].state = worksheet.cell_value(curr_row, c_state)
                                if c_elev >= 0:
                                    self.towns[-1].elev = worksheet.cell_value(curr_row, c_elev)
                                if c_ctry >= 0:
                                    self.towns[-1].country = worksheet.cell_value(curr_row, c_ctry)
# Process list of towns second
        if os.path.exists(self.town_file):
            read_town_csv(self.town_file)

    def SAM_Header(self, lat, lon):
#   SAM CSV (Solar)
#   Location ID,City,State,Country,Latitude,Longitude,Elevation,Time Zone,Source
        hdr = ''
        distance = 999999
        for town in self.towns:
            dist = self.haversine(lat, lon, town.lat, town.lon)
            if dist < distance:
                hdr = 'Location ID,City,State,Country,Latitude,Longitude,Time Zone,Elevation,Source\n' + \
                    '%s,%s,%s,%s,%s,%s,%s,%s,IWEC' % (town.lid, town.name, town.state, town.country,
                    lat, lon, town.zone, town.elev)
                distance = dist
        return hdr

    def SMW_Header(self, lat, lon, year=2014):
#   SAM SMW (Solar)
#   Location ID,Name,State,Time Zone,Latitude,Longitude,Elevation,Time Step,Start year,Start time
        hdr = ''
        distance = 999999
        for town in self.towns:
            dist = self.haversine(lat, lon, town.lat, town.lon)
            if dist < distance:
                hdr = '%s,"%s",%s,%s,%s,%s,%s,3600.0,%s,0:30:00' % (town.lid, town.name, town.state, town.zone,
                    lat, lon, town.elev, year)
                distance = dist
        return hdr

    def SRW_Header(self, lat, lon, year=2014):
#   SRW (Wind)
#   <location id>,<city>,<state>,<country>,<year>,<latitude>,<longitude>,<elevation>,<time step in hours>,<number of rows>
        hdr = ''
        distance = 999999
        for town in self.towns:
            dist = self.haversine(lat, lon, town.lon, town.lat)
            if dist < distance:
                hdr = '%s,%s,%s,%s,%s,%s,%s,%s,1,8760' % (town.lid, town.name, town.state, town.country,
                    year, round(lat, 4), round(lon, 4), town.elev)
                distance = dist
        return hdr

    def Nearest(self, lat, lon, distance=False):
        the_town = ''
        distnce = 999999
        for twn in self.towns:
            dist = self.haversine(lat, lon, twn.lat, twn.lon)
            if dist < distnce:
                the_town = twn.name
                distnce = dist
        for twn in self.towns:
            if twn.name == the_town:
                if distance:
                    return twn, distnce
                else:
                    return twn
        return

    def Stn_Location(self, lid):
        for town in self.towns:
            if town.lid == lid.lstrip('0'):
                return str(town.lat) + ' ' + str(town.lon)
        return ''

    def Get_Town(self, name):
        for town in self.towns:
            if town.name == name:
                return town
        return '' def getModelFile(*args):
    def set_models_locn(ini_file):
        app = QtWidgets.QApplication.instance()
        if app is None:
            app = QtWidgets.QApplication(sys.argv)
        if sys.platform == 'win32' or sys.platform == 'cygwin':
            fldr_div = '\\'
        else:
            fldr_div = '/'
        siren_dir = '.' + fldr_div
        newdir = QtWidgets.QFileDialog.getExistingDirectory(None,
                 'SIREN. Choose location for SIREN Model (Preferences) files',
                 '', QtWidgets.QFileDialog.ShowDirsOnly)
        if newdir != '':
            mydir = os.path.abspath(__file__)
            if sys.platform == 'win32' or sys.platform == 'cygwin':
                newdir = newdir.replace('/', fldr_div)
            mydir = mydir[:mydir.rfind(fldr_div)]
            that_len = len(commonprefix([newdir, mydir], chr=fldr_div))
            if that_len > 0:
                bits = newdir[that_len:].split(fldr_div)
                if that_len < len(mydir): # go up the tree
                    pfx = ('..' + fldr_div) * (len(bits) - 1)
                else:
                    pfx = ''
                siren_dir = pfx + newdir[that_len + 1:]
                if siren_dir == '':
                    siren_dir = '.'
                if siren_dir[-1] != fldr_div:
                    siren_dir += fldr_div
            else:
                siren_dir = newdir
            mf = open('siren_models_location.txt', 'w')
            mf.write(siren_dir)
            mf.close()
            updir = mydir[:mydir.rfind(fldr_div) + 1]
            if os.path.exists(updir + 'siren_sample') \
              and os.path.exists(mydir + fldr_div + 'SIREN.ini') \
              and not os.path.exists(siren_dir + 'SIREN.ini'): # has the sample
                msgbox = QtWidgets.QMessageBox()
                msgbox.setWindowTitle('SIREN - Copy File')
                msgbox.setWindowIcon(QtGui.QIcon('sen_icon32.ico'))
                msgbox.setText('SIREN sample_files folder found')
                msgbox.setInformativeText('Do you want to copy the sample SIREN Model (Y)?')
                msgbox.setDetailedText('It seems you have the siren sample files. ' + \
                    "If you reply 'Y'es the sample preferences file, SIREN.ini, will " + \
                    'be copied to '+ siren_dir + '.')
                msgbox.setIcon(QtWidgets.QMessageBox.Question)
                msgbox.setStandardButtons(QtWidgets.QMessageBox.Yes | QtWidgets.QMessageBox.No)
                reply = msgbox.exec_()
                if reply == QtWidgets.QMessageBox.Yes:
                    if os.path.exists(siren_dir + 'SIREN.ini'):
                        if os.path.exists(siren_dir + 'SIREN.ini~'):
                            os.remove(siren_dir + 'SIREN.ini~')
                        os.rename(siren_dir + 'SIREN.ini', siren_dir + 'SIREN.ini~')
                    copy(mydir + fldr_div +  'SIREN.ini', siren_dir + 'SIREN.ini')
        return siren_dir + ini_file

    siren_dir = ''
    if len(args) > 0:
        ini_file = args[0]
        if os.path.exists('siren_models_location.txt'):
            mf = open('siren_models_location.txt', 'r')
            siren_dir = mf.readline()
            mf.close()
            siren_dir = siren_dir.strip('\n')
            siren_dir = siren_dir + ini_file
            return siren_dir
        else:
            return ini_file
    ini_file = ''
    if os.path.exists('siren_models_location.txt'):
        mf = open('siren_models_location.txt', 'r')
        siren_dir = mf.readline()
        mf.close()
        siren_dir = siren_dir.strip('\n')
        if not(os.path.exists(siren_dir)):
            app = QtWidgets.QApplication.instance()
            if app is None:
                app = QtWidgets.QApplication(sys.argv)
            msgbox = QtWidgets.QMessageBox()
            msgbox.setWindowTitle("SIREN - Can't find Models")
            msgbox.setWindowIcon(QtGui.QIcon('sen_icon32.ico'))
            msgbox.setText('SIREN Models folder missing')
            msgbox.setInformativeText('Do you want to reset the Models location (Y)?')
            msgbox.setDetailedText("Can't find " + siren_dir + '. ' + \
                "If you reply 'Y'es you can choose a new location for the Models.")
            msgbox.setIcon(QtWidgets.QMessageBox.Question)
            msgbox.setStandardButtons(QtWidgets.QMessageBox.Yes | QtWidgets.QMessageBox.No)
            reply = msgbox.exec_()
            if reply == QtWidgets.QMessageBox.Yes:
                siren_dir = set_models_locn(ini_file)
            else:
                sys.exit(8)
        siren_dir = siren_dir + ini_file
    else:
        siren_dir = set_models_locn(ini_file)
    return siren_dir class Grid:
    def decode2(self, proces, substation=False):
        operands = ['+', '-', '*', '/', 'x']
        proces = proces.replace(' ', '')
        proces = proces.split('=')
        numbr = [0., 0.]
        for p in range(len(proces)):
            stack = []
            i = 0
            j = 0
            while j < len(proces[p]):
                if proces[p][j] in operands:
                    stack.append(proces[p][i:j])
                    stack.append(proces[p][j])
                    i = j + 1
                j += 1
            if j > i:
                stack.append(proces[p][i:])
            if p == 1:
                i = 0
                while i < len(stack):
                    try:
                        if substation:
                            stack[i] = self.substation_costs[stack[i]]
                        else:
                            stack[i] = self.line_costs[stack[i]]
                    except:
                        pass
                    i += 2
            numbr[p] = float(stack[0])
            for i in range(1, len(stack), 2):
                if stack[i] == '*' or stack[i] == 'x':
                    numbr[p] = numbr[p] * float(stack[i + 1])
                elif stack[i] == '/':
                    numbr[p] = numbr[p] / float(stack[i + 1])
                if stack[i] == '+':
                    numbr[p] = numbr[p] + float(stack[i + 1])
                if stack[i] == '-':
                    numbr[p] = numbr[p] - float(stack[i + 1])
        return [numbr[0], numbr[1], proces[1]]

    def decode(self, s_lines, the_lines=None):
        if the_lines is None:
            new_lines = []
        else:
            new_lines = the_lines
        l = 0
        while True:
            l = s_lines.find('(', l)
            if l < 0:
                break
            r = s_lines.find(')', l)
            l2 = s_lines.find('(', l + 1)
            if l2 > r:   # normal
                proces = s_lines[l + 1:r]
                l += 1
            else:
                b = 1
                while l2 < r:
                    b += 1
                    l2 = s_lines.find('(', l2 + 1)
                    if l2 < 0:
                        break
                while b > 0:
                    b -= 1
                    r = s_lines.find(')', r + 1)
                proces = s_lines[l + 1:r]
                l = r
            if proces[:4] == 'for(':
                i = 4
                while proces[i] != '=':
                    i += 1
                subs = proces[4:i]
                j = i + 1
                while proces[j] != ',':
                    j += 1
                frm = int(proces[i + 1:j])
                i = j + 1
                while proces[i] != ',':
                    i += 1
                too = int(proces[j + 1:i])
                for j in range(frm, too + 1):
                    this_bit = proces[i + 1:-1].replace(subs, str(j))
                    if j == 1:
                        this_bit = this_bit.replace('1*', '')
                        this_bit = this_bit.replace('*1', '')
                    self.decode(this_bit, new_lines)
            else:
                new_lines.append(self.decode2(proces))
        return new_lines

    def get_config(self):
        config = configparser.RawConfigParser()
        if len(sys.argv) > 1:
            config_file = sys.argv[1]
        else:
            config_file = getModelFile('SIREN.ini')
        config.read(config_file)
        try:
            self.base_year = config.get('Base', 'year')
        except:
            self.base_year = '2012'
        parents = []
        try:
            parents = getParents(aparents = config.items('Parents'))
        except:
            pass
        try:
            self.kml_file = config.get('Files', 'grid_network')
            for key, value in parents:
                self.kml_file = self.kml_file.replace(key, value)
            self.kml_file = self.kml_file.replace('$USER$', getUser())
            self.kml_file = self.kml_file.replace('$YEAR$', self.base_year)
        except:
            self.kml_file = ''
        try:
            self.kml_file2 = config.get('Files', 'grid_network2')
            for key, value in parents:
                self.kml_file2 = self.kml_file2.replace(key, value)
            self.kml_file2 = self.kml_file2.replace('$USER$', getUser())
            self.kml_file2 = self.kml_file2.replace('$YEAR$', self.base_year)
        except:
            self.kml_file2 = ''
        try:
            mapc = config.get('Map', 'map_choice')
        except:
            mapc = ''
        self.colors = {}
        self.colors['grid_boundary'] = 'blue'
        self.grid2_colors = {'500': '#%02x%02x%02x' % (255, 237, 0),
                             '400': '#%02x%02x%02x' % (242, 170, 45),
                             '330': '#%02x%02x%02x' % (242, 170, 45),
                             '275': '#%02x%02x%02x' % (249, 13, 227),
                             '220': '#%02x%02x%02x' % (0, 0, 255),
                             '132': '#%02x%02x%02x' % (228, 2, 45),
                             '110': '#%02x%02x%02x' % (228, 2, 45),
                             '88': '#%02x%02x%02x' % (119, 65, 16),
                             '66': '#%02x%02x%02x' % (119, 65, 16),
                             'dc': '#%02x%02x%02x' % (86, 154, 105)}
        try:
            colours = config.items('Colors')
            for item, colour in colours:
                if item in ['grid_boundary', 'grid_trace']:
                    self.colors[item] = colour
                elif item[:5] == 'grid_':
                    self.colors[item[5:]] = colour
                elif item[:6] == 'grid2_':
                    self.grid2_colors[item[6:]] = colour
        except:
            pass
        if mapc != '':
            try:
                colours = config.items('Colors' + mapc)
                for item, colour in colours:
                    if item in ['grid_boundary', 'grid_trace']:
                        self.colors[item] = colour
                    elif item[:5] == 'grid_':
                        self.colors[item[5:]] = colour
                    elif item[:6] == 'grid2_':
                        self.grid2_colors[item[6:]] = colour
            except:
                pass
        upper_left = [0., 0.]
        lower_right = [-90., 180.]
        try:
             upper_left = config.get('Map', 'upper_left' + mapc).split(',')
             upper_left[0] = float(upper_left[0].strip())
             upper_left[1] = float(upper_left[1].strip())
             lower_right = config.get('Map', 'lower_right' + mapc).split(',')
             lower_right[0] = float(lower_right[0].strip())
             lower_right[1] = float(lower_right[1].strip())
        except:
             try:
                 lower_left = config.get('Map', 'lower_left' + mapc).split(',')
                 upper_right = config.get('Map', 'upper_right' + mapc).split(',')
                 upper_left[0] = float(upper_right[0].strip())
                 upper_left[1] = float(lower_left[1].strip())
                 lower_right[0] = float(lower_left[0].strip())
                 lower_right[1] = float(upper_right[1].strip())
             except:
                 pass
        self.map_polygon = [upper_left, [upper_left[0], lower_right[1]], lower_right,
           [lower_right[0], upper_left[1]], upper_left]
        self.default_length = -1
        try:
            trace_existing = config.get('View', 'trace_existing')
            if trace_existing.lower() in ['true', 'yes', 'on']:
                self.default_length = 0
        except:
            try:
                trace_existing = config.get('Grid', 'trace_existing')
                if trace_existing.lower() in ['true', 'yes', 'on']:
                   self.default_length = 0
            except:
                pass
        self.line_costs = {}
        try:
            line_costs = config.get('Grid', 'line_costs')
            line_costs = line_costs.replace('(', '')
            line_costs = line_costs.replace(')', '')
            line_costs = line_costs.split(',')
            for i in range(len(line_costs)):
                line = line_costs[i].split('=')
                if line[1].strip()[-1] == 'K':
                    self.line_costs[line[0].strip()] = float(line[1].strip()[:-1]) * pow(10, 3)
                elif line[1].strip()[-1] == 'M':
                    self.line_costs[line[0].strip()] = float(line[1].strip()[:-1]) * pow(10, 6)
                else:
                    self.line_costs[line[0].strip()] = float(line[1].strip()[:-1])
        except:
            pass
        self.substation_costs = {}
        try:
            substation_costs = config.get('Grid', 'substation_costs')
            substation_costs = substation_costs.replace('(', '')
            substation_costs = substation_costs.replace(')', '')
            substation_costs = substation_costs.split(',')
            for i in range(len(substation_costs)):
                line = substation_costs[i].split('=')
                if line[1].strip()[-1] == 'K':
                    self.substation_costs[line[0].strip()] = float(line[1].strip()[:-1]) * pow(10, 3)
                elif line[1].strip()[-1] == 'M':
                    self.substation_costs[line[0].strip()] = float(line[1].strip()[:-1]) * pow(10, 6)
                else:
                    self.substation_costs[line[0].strip()] = float(line[1].strip()[:-1])
        except:
            pass
        try:
            s_lines = config.get('Grid', 's_lines')
            s_lines = s_lines.strip()
            self.s_line_table = self.decode(s_lines)
        except:
            self.s_line_table = []
        try:
            d_lines = config.get('Grid', 'd_lines')
            d_lines = d_lines.strip()
            self.d_line_table = self.decode(d_lines)
        except:
            self.d_line_table = []
        self.dummy_fix = False
        try:
            dummy_fix = config.get('Grid', 'dummy_fix')
            if dummy_fix.lower() in ['true', 'yes', 'on']:
                self.dummy_fix = True
        except:
            pass

    def __init__(self, grid2=False):
        self.get_config()
        self.lines = []
        if grid2:
            kml_file = self.kml_file2
        else:
            kml_file = self.kml_file
        if not os.path.exists(kml_file):
            if grid2:
                self.kml_file = ''
            else:
                self.kml_file2 = ''
            return
        style = {}
        styl = ''
        zipped = False
        if kml_file[-4:] == '.kmz': # zipped file?
            zipped = True
            zf = zipfile.ZipFile(kml_file, 'r')
            inner_file = ''
            for name in zf.namelist():
                if name[-4:] == '.kml':
                    inner_file = name
                    break
            if inner_file == '':
                return
            memory_file = io.BytesIO()
            memory_file.write(zf.open(inner_file).read())
            root = ElementTree(fromstring(memory_file.getvalue()))
        else:
            kml_data = open(kml_file, 'rb')
            root = ElementTree(fromstring(kml_data.read()))
         # Create an iterator
        if sys.version_info[1] < 9: # before python 3.9
            iterat = root.getiterator()
        else:
            iterat = root.iter()
        placemark_id = ''
        line_names = []
        stylm = ''
        for element in iterat:
            elem = element.tag[element.tag.find('}') + 1:]
            if elem == 'Style':
                for name, value in list(element.items()):
                    if name == 'id':
                        styl = value
            elif elem == 'StyleMap':
                for name, value in list(element.items()):
                    if name == 'id':
                        stylm = value
            elif elem == 'color':
                if styl in self.colors:
                    style[styl] = self.colors[styl]
                else:
                    style[styl] = '#' + element.text[-2:] + element.text[-4:-2] + element.text[-6:-4]
                if stylm != '':
                    if stylm in self.colors:
                        style[stylm] = self.colors[stylm]
                    else:
                        style[stylm] = '#' + element.text[-2:] + element.text[-4:-2] + element.text[-6:-4]
            elif elem == 'name':
                line_name = element.text
                if placemark_id != '':
                    line_name += placemark_id
                    placemark_id = ''
            elif elem == 'Placemark' and grid2:
                for key, value in list(element.items()):
                    if key == 'id':
                        if value[:4] == 'kml_':
                            placemark_id = value[3:]
                        else:
                            placemark_id = value
            elif elem == 'SimpleData' and grid2:
                for key, value in list(element.items()):
                    if key == 'name' and (value == 'CAPACITY_kV' or value == 'CAPACITYKV'):
                        try:
                            styl = self.grid2_colors[element.text]
                        except:
                            styl = self.grid2_colors['66']
            elif elem == 'styleUrl':
                styl = element.text[1:]
            elif elem == 'coordinates':
                coords = []
                coordinates = ' '.join(element.text.split()).split()
                for i in range(len(coordinates)):
                    coords.append([float(coordinates[i].split(',')[1]), float(coordinates[i].split(',')[0])])
                inmap = False
                for coord in coords:
                    if within_map(coord[0], coord[1], self.map_polygon):
                        inmap = True
                        break
                if inmap:
                    if self.default_length >= 0:
                        grid_len = 0.
                        for j in range(1, len(coords)):
                            grid_len += self.actualDistance(coords[j - 1][0], coords[j - 1][1],
                                        coords[j][0], coords[j][1])
                    else:
                        grid_len = self.default_length
                    if line_name in line_names:
                        i = 2
                        while line_name + '#' + str(i) in line_names:
                            i += 1
                        line_name += '#' + str(i)
                    line_names.append(line_name)
                    if grid2:
                        self.lines.append(Line(line_name, styl, coords, length=grid_len))
                    else:
                        try:
                            self.lines.append(Line(line_name, style[styl], coords, length=grid_len))
                        except:
                            style[styl] = '#FFFFFF'
                            self.lines.append(Line(line_name, style[styl], coords, length=grid_len))
        if zipped:
            memory_file.close()
            zf.close()
        else:
            kml_data.close()
     # connect together
     # if load_centres connect closest end to closest load centre
     #    for i in range(len(self.lines)):
     #        connect = []
     #        con = -1
     #        connect.append(self.gridConnect(self.lines[i].coordinates[0][0], self.lines[i].coordinates[0][1], \
     #                       ignore=[i]))
     #        connect.append(self.gridConnect(self.lines[i].coordinates[-1][0], self.lines[i].coordinates[-1][1], \
     #                       ignore=[i]))
     #        if connect[0][0] > 0:
     #            if connect[1][0] < connect[0][0]:
     #                con = connect[1][2]
     #                self.lines[i].coordinates.append([connect[1][1], connect[1][2]])
     #            else:
     #                con = connect[0][2]
     #                self.lines[i].coordinates.insert(0, [connect[0][1], connect[0][2]])
     #        self.lines[i].connector = con

    def gridConnect(self, lat, lon, ignore=[]):
        shortest = [99999, -1., -1., -1]
        for l in range(len(self.lines)):
            if l in ignore:
                continue
            for i in range(len(self.lines[l].coordinates) - 1):
                if self.dummy_fix:
                    dist = dust(lat, lon, self.lines[l].coordinates[i][0], self.lines[l].coordinates[i][1],
                           self.lines[l].coordinates[i + 1][0], self.lines[l].coordinates[i + 1][1])
                else:
                    dist = self.DistancePointLine(lat, lon, self.lines[l].coordinates[i][0], self.lines[l].coordinates[i][1],
                           self.lines[l].coordinates[i + 1][0], self.lines[l].coordinates[i + 1][1])
                if dist[0] >= 0 and dist[0] < shortest[0]:
                    shortest = dist[:]
                    shortest.append(l)
        if shortest[0] == 99999:
             shortest[0] = -1
        return shortest   # length, lat, lon, line#

    def DistancePointLine(self, pyd, pxd, y1d, x1d, y2d, x2d):
# px,py is the point to test.
# x1,y1,x2,y2 is the line to check distance.
# Returns distance from the line, or if the intersecting point on the line nearest
# the point tested is outside the endpoints of the line, the distance to the
# nearest endpoint.
# Returns -1 on 0 denominator conditions.
        px = radians(pxd)
        py = radians(pyd)
        x1 = radians(x1d)
        y1 = radians(y1d)
        x2 = radians(x2d)
        y2 = radians(y2d)
        b13 = self.Bearing(y1, x1, py, px)
        b12 = self.Bearing(y1, x1, y2, x2)
        d13 = self.Distance(y1, x1, py, px)
        d23 = self.Distance(y2, x2, py, px)
        dxt = asin(sin(d13) * sin(b13 - b12))
        dat = acos(cos(d13) / cos(dxt))
        iy = asin(sin(y1) * cos(dat) + cos(y1) * sin(dat) * cos(b12))
        ix = x1 + atan2(sin(b12) * sin(dat) * cos(y1), cos(dat) - sin(y1) * sin(iy))
        if abs(ix - x1) > abs(x1 - x2) or abs(iy - y1) > abs(y1 - y2):
            dst = d13
            ix = x1
            iy = y1
        else:
            dst = self.Distance(iy, ix, py, px)
            if d13 < dst:   # must be another way but this'll do for now
                dst = d13
                ix = x1
                iy = y1
        if d23 < dst:   # must be another way but this'll do for now
            dst = d23
            ix = x2
            iy = y2
        return [round(abs(dst) * RADIUS, 2), round(degrees(iy), 6), round(degrees(ix), 6)]

    def Bearing(self, y1, x1, y2, x2):
# find the bearing between the coordinates
        return atan2(sin(x2 - x1) * cos(y2), cos(y1) * sin(y2) - sin(y1) * cos(y2) * cos(x2 - x1))

    def Distance(self, y1, x1, y2, x2):
# find the differences between the coordinates
        dy = y2 - y1
        dx = x2 - x1
        ra13 = pow(sin(dy / 2.), 2) + cos(y1) * cos(y2) * pow(sin(dx / 2.), 2)
        return 2 * asin(min(1, sqrt(ra13)))

    def actualDistance(self, y1d, x1d, y2d, x2d):
        x1 = radians(x1d)
        y1 = radians(y1d)
        x2 = radians(x2d)
        y2 = radians(y2d)
        dst = self.Distance(y1, x1, y2, x2)
        return round(abs(dst) * RADIUS, 2)

    def Line_Cost(self, peak_load, peak_dispatchable):
        if peak_load is None:
            peak_load = 0
        if peak_dispatchable is None or peak_dispatchable == 0:
            s = 0
            while self.s_line_table[s][0] < peak_load and (s + 1) < len(self.s_line_table):
                s += 1
            return self.s_line_table[s][1], self.s_line_table[s][2]
        elif peak_load == peak_dispatchable:
            p = 0
            while self.d_line_table[p][0] < peak_load and (p + 1) < len(self.d_line_table):
                p += 1
            return self.d_line_table[p][1], self.d_line_table[p][2]
        else:
            if peak_load > peak_dispatchable + peak_dispatchable:
                p = 0
                while self.d_line_table[p][0] < peak_load and (p + 1) < len(self.d_line_table):
                    p += 1
                return self.d_line_table[p][1], self.d_line_table[p][2]
            p = 0
            while self.d_line_table[p][0] < peak_dispatchable and (p + 1) < len(self.d_line_table):
                p += 1
            s = 0
            while self.s_line_table[s][0] < peak_load - peak_dispatchable and (s + 1) < len(self.s_line_table):
                s += 1
            if self.d_line_table[p][0] >= self.s_line_table[s][0]:
                return self.d_line_table[p][1], self.d_line_table[p][2]
            else:
                return self.d_line_table[p][1], self.d_line_table[p][2]

    def Substation_Cost(self, line):
        if line in self.substation_costs:
            return self.substation_costs[line]
        return 0. class Grid_Boundary:
    def get_config(self):
        config = configparser.RawConfigParser()
        if len(sys.argv) > 1:
            config_file = sys.argv[1]
        else:
            config_file = getModelFile('SIREN.ini')
        config.read(config_file)
        try:
            self.base_year = config.get('Base', 'year')
        except:
            self.base_year = '2012'
        parents = []
        try:
            parents = getParents(config.items('Parents'))
        except:
            pass
        try:
            self.kml_file = config.get('Files', 'grid_boundary')
            for key, value in parents:
                self.kml_file = self.kml_file.replace(key, value)
            self.kml_file = self.kml_file.replace('$USER$', getUser())
            self.kml_file = self.kml_file.replace('$YEAR$', self.base_year)
        except:
            self.kml_file = ''
        self.colour = '#0000FF'
        try:
            mapc = config.get('Map', 'map_choice')
        except:
            mapc = ''
        try:
            self.colour = config.get('Colors', 'grid_boundary')
        except:
            pass
        if mapc != '':
            try:
                self.colour = config.get('Colors' + mapc, 'grid_boundary')
            except:
                pass
        upper_left = [0., 0.]
        lower_right = [-90., 180.]
        try:
             upper_left = config.get('Map', 'upper_left' + mapc).split(',')
             upper_left[0] = float(upper_left[0].strip())
             upper_left[1] = float(upper_left[1].strip())
             lower_right = config.get('Map', 'lower_right' + mapc).split(',')
             lower_right[0] = float(lower_right[0].strip())
             lower_right[1] = float(lower_right[1].strip())
        except:
             try:
                 lower_left = config.get('Map', 'lower_left' + mapc).split(',')
                 upper_right = config.get('Map', 'upper_right' + mapc).split(',')
                 upper_left[0] = float(upper_right[0].strip())
                 upper_left[1] = float(lower_left[1].strip())
                 lower_right[0] = float(lower_left[0].strip())
                 lower_right[1] = float(upper_right[1].strip())
             except:
                 pass
        self.map_polygon = [upper_left, [upper_left[0], lower_right[1]], lower_right,
           [lower_right[0], upper_left[1]], upper_left]

    def __init__(self):
        self.get_config()
        self.lines = []
        if not os.path.exists(self.kml_file):
            return
        style = {}
        styl = ''
        zipped = False
        if self.kml_file[-4:] == '.kmz': # zipped file?
            zipped = True
            zf = zipfile.ZipFile(kml_file, 'r')
            inner_file = ''
            for name in zf.namelist():
                if name[-4:] == '.kml':
                    inner_file = name
                    break
            if inner_file == '':
                return
            memory_file = io.BytesIO()
            memory_file.write(zf.open(inner_file).read())
            root = ElementTree(fromstring(memory_file.getvalue()))
        else:
            kml_data = open(self.kml_file, 'rb')
            root = ElementTree(fromstring(kml_data.read()))
         # Create an iterator
        if sys.version_info[1] < 9: # before python 3.9
            iterat = root.getiterator()
        else:
            iterat = root.iter()
        for element in iterat:
            elem = element.tag[element.tag.find('}') + 1:]
            if elem == 'Style':
                for name, value in list(element.items()):
                    if name == 'id':
                        styl = value
            elif elem == 'color':
                style[styl] = self.colour
            elif elem == 'name':
                line_name = element.text
            elif elem == 'styleUrl':
                styl = element.text[1:]
            elif elem == 'coordinates':
                coords = []
                coordinates = ' '.join(element.text.split()).split()
                for i in range(len(coordinates)):
                    coords.append([round(float(coordinates[i].split(',')[1]), 6),
                      round(float(coordinates[i].split(',')[0]), 6)])
                i = int(len(coords) / 2)
                if within_map(coords[0][0], coords[0][1], self.map_polygon) and \
                   within_map(coords[i][0], coords[i][1], self.map_polygon):
                    try:
                        self.lines.append(Line(line_name, style[styl], coords))
                    except:
                        self.lines.append(Line(line_name, self.colour, coords))
        if zipped:
            memory_file.close()
            zf.close()
        else:
            kml_data.close() class Line:
    def __init__(self, name, style, coordinates, length=0., connector=-1, dispatchable=None, line_cost=None, peak_load=None,
                 peak_dispatchable=None, peak_loss=None, line_table=None, substation_cost=None):
        self.name = name
        self.style = style
        self.coordinates = []
        for i in range(len(coordinates)):
            self.coordinates.append([])
            for j in range(len(coordinates[i])):
                self.coordinates[-1].append(round(coordinates[i][j], 6))
        self.connector = connector
        self.length = length
        self.dispatchable = dispatchable
        self.line_cost = line_cost
        self.peak_load = peak_load
        self.peak_dispatchable = peak_dispatchable
        self.peak_loss = peak_loss
        self.line_table = line_table
        self.substation_cost = substation_cost def getParents(aparents):
    parents = []
    for key, value in aparents:
        for key2, value2 in aparents:
            if key2 == key:
                continue
            value = value.replace(key2, value2)
        for key2, value2 in parents:
            if key2 == key:
                continue
            value = value.replace(key2, value2)
        parents.append((key, value))
    return parents def getUser():
    if sys.platform == 'win32' or sys.platform == 'cygwin':   # windows
        return os.environ.get("USERNAME")

    elif sys.platform == 'darwin':   # osx64
        return os.environ.get("USERNAME")
    elif sys.platform == 'linux' or sys.platform == 'linux2':   # linux
        return pwd.getpwuid(os.geteuid()).pw_name
    else:
        return os.environ.get("USERNAME") def techClean(tech, full=False):
    cleantech = tech.replace('_', ' ').title()
    cleantech = cleantech.replace('Bm', 'BM')
    cleantech = cleantech.replace('Ccgt', 'CCGT')
    cleantech = cleantech.replace('Ccg', 'CCG')
    cleantech = cleantech.replace('Cst', 'CST')
    cleantech = cleantech.replace('Lng', 'LNG')
    cleantech = cleantech.replace('Ocgt', 'OCGT')
    cleantech = cleantech.replace('Ocg', 'OCG')
    cleantech = cleantech.replace('Phs', 'PHS')
    cleantech = cleantech.replace('Pv', 'PV')
    if full:
        alll = [['Cf', 'CF'], ['hi ', 'HI '], ['Lcoe', 'LCOE'], ['Mw', 'MW'],
                ['ni ', 'NI '], ['Npv', 'NPV'], ['Re', 'RE'], ['Tco2E', 'tCO2e']]
        for each in alll:
            cleantech = cleantech.replace(each[0], each[1])
    return cleantech class Station:
    def __init__(self, name, technology, lat, lon, capacity, turbine, rotor, no_turbines, area, scenario, generation=None,
                 power_file=None, grid_line=None, grid_len=None, grid_path_len=None, direction=None, tilt=None, storage_hours=None):
        self.name = name
        self.technology = technology
        self.lat = lat
        self.lon = lon
        self.capacity = capacity
        self.turbine = turbine
        self.rotor = rotor
        self.no_turbines = no_turbines
        self.area = area
        self.scenario = scenario
        self.generation = generation
        self.power_file = power_file
        self.grid_line = grid_line
        self.grid_len = grid_len
        self.grid_path_len = grid_path_len
        self.direction = direction
        self.storage_hours = storage_hours
        if tilt is not None:
            self.tilt = tilt class Stations:
    def get_config(self):
        config = configparser.RawConfigParser()
        if __name__ == '__main__':
            for i in range(1, len(sys.argv)):
                if sys.argv[i][-4:] == '.ini':
                    config_file = sys.argv[i]
                    break
        else:
            if len(sys.argv) > 1:
                config_file = sys.argv[1]
            else:
                config_file = getModelFile('SIREN.ini')
        config.read(config_file)
        try:
            self.base_year = config.get('Base', 'year')
        except:
            self.base_year = '2012'
        parents = []
        try:
            parents = getParents(config.items('Parents'))
        except:
            pass
        try:
            self.sam_file = config.get('Files', 'sam_turbines')
            for key, value in parents:
                self.sam_file = self.sam_file.replace(key, value)
            self.sam_file = self.sam_file.replace('$USER$', getUser())
            self.sam_file = self.sam_file.replace('$YEAR$', self.base_year)
        except:
            self.sam_file = ''
        try:
            self.pow_dir = config.get('Files', 'pow_files')
            for key, value in parents:
                self.pow_dir = self.pow_dir.replace(key, value)
            self.pow_dir = self.pow_dir.replace('$USER$', getUser())
            self.pow_dir = self.pow_dir.replace('$YEAR$', self.base_year)
        except:
            self.pow_dir = ''
        self.fac_files = []
        try:
            fac_file = config.get('Files', 'grid_stations')
            for key, value in parents:
                fac_file = fac_file.replace(key, value)
            fac_file = fac_file.replace('$USER$', getUser())
            fac_file = fac_file.replace('$YEAR$', self.base_year)
            self.fac_files.append(fac_file)
        except:
            pass
        if self.stations2:
            try:
                fac_file = config.get('Files', 'grid_stations2')
                for key, value in parents:
                    fac_file = fac_file.replace(key, value)
                fac_file = fac_file.replace('$USER$', getUser())
                fac_file = fac_file.replace('$YEAR$', self.base_year)
                self.fac_files.append(fac_file)
            except:
                pass
        self.ignore_deleted = True
        try:
            if config.get('Grid', 'ignore_deleted_existing').lower() in ['false', 'off', 'no']:
                self.ignore_deleted = False
        except:
            pass
        self.technologies = ['']
        self.areas = {}
        try:
            technologies = config.get('Power', 'technologies')
            for item in technologies.split():
                itm = techClean(item)
                self.technologies.append(itm)
                try:
                    self.areas[itm] = float(config.get(itm, 'area'))
                except:
                    self.areas[itm] = 0.
        except:
            pass
        try:
            technologies = config.get('Power', 'fossil_technologies')
            technologies = technologies.split()
            for item in technologies:
                itm = techClean(item)
                try:
                    self.areas[itm] = float(config.get(itm, 'area'))
                except:
                    self.areas[itm] = 0.
        except:
            pass
        try:
            mapc = config.get('Map', 'map_choice')
        except:
            mapc = ''
        upper_left = [0., 0.]
        lower_right = [-90., 180.]
        try:
             upper_left = config.get('Map', 'upper_left' + mapc).split(',')
             upper_left[0] = float(upper_left[0].strip())
             upper_left[1] = float(upper_left[1].strip())
             lower_right = config.get('Map', 'lower_right' + mapc).split(',')
             lower_right[0] = float(lower_right[0].strip())
             lower_right[1] = float(lower_right[1].strip())
        except:
             try:
                 lower_left = config.get('Map', 'lower_left' + mapc).split(',')
                 upper_right = config.get('Map', 'upper_right' + mapc).split(',')
                 upper_left[0] = float(upper_right[0].strip())
                 upper_left[1] = float(lower_left[1].strip())
                 lower_right[0] = float(lower_left[0].strip())
                 lower_right[1] = float(upper_right[1].strip())
             except:
                 pass
        self.map_polygon = [upper_left, [upper_left[0], lower_right[1]], lower_right,
           [lower_right[0], upper_left[1]], upper_left]

    def haversine(self, lat1, lon1, lat2, lon2):
        """
        Calculate the great circle distance between two points
        on the earth (specified in decimal degrees)
        """
     # convert decimal degrees to radians
        lon1, lat1, lon2, lat2 = list(map(radians, [lon1, lat1, lon2, lat2]))

     # haversine formula
        dlon = lon2 - lon1
        dlat = lat2 - lat1
        a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
        c = 2 * asin(sqrt(a))

     # 6367 km is the radius of the Earth
        km = 6367 * c
        return km

    def __init__(self, code=False, fossil=False, existing=True, stations2=True):
        self.stations2 = stations2
        self.get_config()
        try:
            self.stations in locals()
        except:
            self.stations = []
        if not existing:
            return
        if not os.path.exists(self.sam_file):
            return
        sam = open(self.sam_file)
        sam_turbines = csv.DictReader(sam)
        for fac_file in self.fac_files:
            if os.path.exists(fac_file):
                if fac_file[-4:] == '.csv':
                    facile = open(fac_file)
                    facilities = csv.DictReader(facile)
                    for facility in facilities:
                        if self.ignore_deleted and facility['Balancing Status'] == 'Deleted':
                            continue
                        if facility['Longitude'] != '':
                            if not within_map(float(facility['Latitude']),
                              float(facility['Longitude']), self.map_polygon):
                                continue
                            if 'Facility Code' in facilities.fieldnames:   # IMO format
                                bit = facility['Facility Code'].split('_')
                                rotor = 0.
                                turbine = ''
                                no_turbines = 0
                                if bit[-1][:2] == 'WW' or bit[-1][:2] == 'WF':
                                    tech = 'Wind'
                                    turbine = facility['Turbine']
                                    if turbine[:7] == 'Enercon':
                                        bit = turbine[9:].split(' _')
                                        if len(bit) == 1:
                                            rotor = bit[0].split('_')[0]
                                        else:
                                            bit = bit[0].split('_')
                                            if len(bit) == 1:
                                                rotor = bit[0]
                                            else:
                                                if bit[1][-1] == 'm':
                                                    rotor = bit[1][:-1]
                                                else:
                                                    rotor = bit[0]
                                    else:
                                        turb = turbine.split(';')
                                        if len(turb) > 1:
                                            turbine = turb[1]
                                        else:
                                            turbine = turb[0]
                                        sam.seek(0)
                                        for turb in sam_turbines:
                                            if turb['Name'] == turbine:
                                                rotor = turb['Rotor Diameter']
                                                break
                                        else: # try and find .pow file
                                            pow_file = self.pow_dir + '/' + turbine + '.pow'
                                            if os.path.exists(pow_file):
                                                tf = open(pow_file, 'r')
                                                lines = tf.readlines()
                                                tf.close()
                                                rotor = float(lines[1].strip('" \t\n\r'))
                                                del lines
                                    no_turbines = int(facility['No. turbines'])
                                    try:
                                        rotor = float(rotor)
                                    except:
                                        pass
                                    area = self.areas[tech] * float(no_turbines) * pow((rotor * .001), 2)
                                elif bit[0] == 'MERSOLAR':
                                    tech = 'Tracking PV'
                                    area = self.areas[tech] * float(facility['Maximum Capacity (MW)'])
                                elif bit[-1][:2] == 'PV':
                                    tech = 'Fixed PV'
                                    area = self.areas[tech] * float(facility['Maximum Capacity (MW)'])
                                elif bit[-1] == 'PLANT' and bit[-2] == 'BIOMASS':
                                    tech = 'Biomass'
                                    area = self.areas[tech] * float(facility['Maximum Capacity (MW)'])
                                elif facility['Participant Code'] in ['CTE', 'LNDFLLGP', 'PERTHNRGY', 'WGRES']:
                                    tech = 'Biomass'
                                    area = self.areas[tech] * float(facility['Maximum Capacity (MW)'])
                                else:
                                    tech = 'Fossil '
                                    if 'Fossil' in facilities.fieldnames:
                                        tech += facility['Fossil']
                                    else:
                                        if bit[-1][:2] == 'GT' or bit[-1][0] == 'U':
                                            tech += 'OCGT'
                                        elif bit[-1][:2] == 'CC':
                                            tech += 'CCGT'
                                        elif bit[-1][:3] == 'COG':
                                            tech += 'Cogen'
                                        elif bit[0] == 'TESLA' or bit[-1] == 'WGP':
                                            tech += 'Distillate'
                                        elif bit[-2] == 'WGP':
                                            tech += 'Mixed'
                                        elif bit[-1][0] == 'G' and (bit[-1][1] >= '1' and bit[-1][1] <= '9'):
                                            tech += 'Coal'
                                        else:
                                            tech += 'Mixed'
                                           #  tech += 'Distillate'
                                    try:
                                        area = self.areas[tech] * float(facility['Maximum Capacity (MW)'])
                                    except:
                                        area = 0
                                if code:
                                    nice_name = facility['Facility Code']
                                else:
                                    nice_name = facility['Facility Name']
                                    if nice_name == '':
                                        name_split = facility['Facility Code'].split('_')
                                        if len(name_split) > 1:
                                            nice_name = ''
                                            for i in range(len(name_split) - 1):
                                                nice_name += name_split[i].title() + '_'
                                            nice_name += name_split[-1]
                                        else:
                                            nice_name = facility['Facility Code']
                                stn = self.Get_Station(nice_name)
                                if stn is None:   # new station?
                                    self.stations.append(Station(nice_name, tech,
                                        float(facility['Latitude']), float(facility['Longitude']),
                                        float(facility['Maximum Capacity (MW)']), turbine, rotor, no_turbines, area, 'Existing'))
                                    if tech == 'Fixed PV':
                                        try:
                                            if facility['Tilt'] != '':
                                                self.stations[-1].tilt = float(facility['Tilt'])
                                        except:
                                            pass
                                else:   # additional generator in existing station
                                    if stn.technology != tech:
                                        if stn.technology[:6] == 'Fossil' and tech[:6] == 'Fossil':
                                            stn.technology = 'Fossil Mixed'
                                    stn.capacity = stn.capacity + float(facility['Maximum Capacity (MW)'])
                                    stn.area += area
                                    stn.no_turbines = stn.no_turbines + no_turbines
                            else:   # SIREN format
                                try:
                                    turbs = int(facility['No. turbines'])
                                except:
                                    turbs = 0
                                self.stations.append(Station(facility['Station Name'],
                                             facility['Technology'],
                                             float(facility['Latitude']),
                                             float(facility['Longitude']),
                                             float(facility['Maximum Capacity (MW)']),
                                             facility['Turbine'],
                                             0.,
                                             turbs,
                                             float(facility['Area']),
                                             'Existing'))
                                if 'Wind' in self.stations[-1].technology:
                                    if self.stations[-1].rotor == 0 or self.stations[-1].rotor == '':
                                        rotor = 0
                                        if self.stations[-1].turbine[:7] == 'Enercon':
                                            bit = self.stations[-1].turbine[9:].split(' _')
                                            if len(bit) == 1:
                                                rotor = bit[0].split('_')[0]
                                            else:
                                                bit = bit[0].split('_')
                                                if len(bit) == 1:
                                                    rotor = bit[0]
                                                else:
                                                    if bit[1][-1] == 'm':
                                                        rotor = bit[1][:-1]
                                                    else:
                                                        rotor = bit[0]
                                        else:
                                            turb = self.stations[-1].turbine.split(';')
                                            if len(turb) > 1:
                                                sam.seek(0)
                                                turbine = turb[1]
                                                for turb in sam_turbines:
                                                    if turb['Name'] == turbine:
                                                        rotor = turb['Rotor Diameter']
                                    try:
                                        self.stations[-1].rotor = float(rotor)
                                    except:
                                        pass
                                if self.stations[-1].area == 0 or self.stations[-1].area == '':
                                    if 'Wind' in self.stations[-1].technology:
                                        self.stations[-1].area = self.areas[self.stations[-1].technology] * \
                                                                 float(self.stations[-1].no_turbines) * \
                                                                 pow((self.stations[-1].rotor * .001), 2)
                                    else:
                                        self.stations[-1].area = self.areas[self.stations[-1].technology] * \
                                                                 float(self.stations[-1].capacity)
                                try:
                                    if facility['Power File'] != '':
                                        self.stations[-1].power_file = facility['Power File']
                                except:
                                    pass
                                try:
                                    if facility['Grid Line'] != '':
                                        self.stations[-1].grid_line = facility['Grid Line']
                                except:
                                    pass
                                if 'PV' in self.stations[-1].technology:
                                    try:
                                        if facility['Direction'] != '':
                                            self.stations[-1].direction = facility['Direction']
                                    except:
                                        pass
                                    try:
                                        if facility['Tilt'] != '':
                                            self.stations[-1].tilt = float(facility['Tilt'])
                                    except:
                                        pass
                                if self.stations[-1].technology in ['CST', 'Solar Thermal']:
                                    try:
                                        if facility['Storage Hours'] != '':
                                            self.stations[-1].storage_hours = float(facility['Storage Hours'])
                                    except:
                                        pass
                    facile.close()
                else:   # assume excel and in our format
                    var = {}
                    workbook = xlrd.open_workbook(fac_file)
                    worksheet = workbook.sheet_by_index(0)
                    num_rows = worksheet.nrows - 1
                    num_cols = worksheet.ncols - 1
                    if worksheet.cell_value(0, 0) == 'Description:' or worksheet.cell_value(0, 0) == 'Comment:':
                        curr_row = 1
                        self.description = worksheet.cell_value(0, 1)
                    else:
                        curr_row = 0
                        self.description = ''
#                   get column names
                    curr_col = -1
                    while curr_col < num_cols:
                        curr_col += 1
                        var[worksheet.cell_value(curr_row, curr_col)] = curr_col
                    while curr_row < num_rows:
                        curr_row += 1
                        if not within_map(worksheet.cell_value(curr_row, var['Latitude']),
                          worksheet.cell_value(curr_row, var['Longitude']), self.map_polygon):
                            continue
                        self.stations.append(Station(str(worksheet.cell_value(curr_row, var['Station Name'])),
                                             str(worksheet.cell_value(curr_row, var['Technology'])),
                                             worksheet.cell_value(curr_row, var['Latitude']),
                                             worksheet.cell_value(curr_row, var['Longitude']),
                                             worksheet.cell_value(curr_row, var['Maximum Capacity (MW)']),
                                             str(worksheet.cell_value(curr_row, var['Turbine'])),
                                             worksheet.cell_value(curr_row, var['Rotor Diam']),
                                             worksheet.cell_value(curr_row, var['No. turbines']),
                                             worksheet.cell_value(curr_row, var['Area']),
                                             'Existing'))
                        if 'Wind' in self.stations[-1].technology:
                            if self.stations[-1].rotor == 0 or self.stations[-1].rotor == '':
                                rotor = 0
                                if self.stations[-1].turbine[:7] == 'Enercon':
                                    bit = self.stations[-1].turbine[9:].split(' _')
                                    if len(bit) == 1:
                                        rotor = bit[0].split('_')[0]
                                    else:
                                        bit = bit[0].split('_')
                                        if len(bit) == 1:
                                            rotor = bit[0]
                                        else:
                                            if bit[1][-1] == 'm':
                                                rotor = bit[1][:-1]
                                            else:
                                                rotor = bit[0]
                                else:
                                    turb = self.stations[-1].turbine.split(';')
                                    if len(turb) > 1:
                                        sam.seek(0)
                                        turbine = turb[1]
                                        for turb in sam_turbines:
                                            if turb['Name'] == turbine:
                                                rotor = turb['Rotor Diameter']
                                try:
                                    self.stations[-1].rotor = float(rotor)
                                except:
                                    pass
                        try:
                            if self.stations[-1].area == 0 or self.stations[-1].area == '':
                                if 'Wind' in self.stations[-1].technology:
                                    self.stations[-1].area = self.areas[self.stations[-1].technology] * \
                                                             float(self.stations[-1].no_turbines) * \
                                                             pow((self.stations[-1].rotor * .001), 2)
                                else:
                                    self.stations[-1].area = self.areas[self.stations[-1].technology] * \
                                                             float(self.stations[-1].capacity)
                        except:
                            self.stations[-1].area = 0.
                        try:
                            power_file = worksheet.cell_value(curr_row, var['Power File'])
                            if power_file != '':
                                self.stations[-1].power_file = power_file
                        except:
                            pass
                        try:
                            grid_line = worksheet.cell_value(curr_row, var['Grid Line'])
                            if grid_line != '':
                                self.stations[-1].grid_line = grid_line
                        except:
                            pass
                        if 'PV' in self.stations[-1].technology:
                            try:
                                direction = worksheet.cell_value(curr_row, var['Direction'])
                                if direction != '':
                                    self.stations[-1].direction = direction
                            except:
                                pass
                            try:
                                tilt = worksheet.cell_value(curr_row, var['Tilt'])
                                if tilt != '':
                                    self.stations[-1].tilt = tilt
                            except:
                                pass
                        if self.stations[-1].technology in ['CST', 'Solar Thermal']:
                            try:
                                storage_hours = worksheet.cell_value(curr_row, var['Storage Hours'])
                                if storage_hours != '':
                                    self.stations[-1].storage_hours = storage_hours
                            except:
                                pass
        sam.close()

    def Nearest(self, lat, lon, distance=False, fossil=False, ignore=None):
        hdr = ''
        distnce = 999999
        for station in self.stations:
            if station.technology[:6] == 'Fossil' and not fossil:
                continue
            dist = self.haversine(lat, lon, station.lat, station.lon)
            if dist < distnce:
                if ignore is not None and ignore == station.name:
                    continue
                hdr = station.name
                distnce = dist
        for station in self.stations:
            if station.name == hdr:
                if distance:
                    return station, distnce
                else:
                    return station
        return None

    def Stn_Location(self, name):
        for station in self.stations:
            if station.name == name:
                return str(station.lat) + ' ' + str(station.lon)
        return ''

    def Stn_Turbine(self, name):
        for station in self.stations:
            if station.name == name:
                return station.turbine
        return ''

    def Get_Station(self, name):
        for station in self.stations:
            if station.name == name:
                return station
        return None

    def Description(self):
        return self.description class Shortest:
    def isBetween(self, a, b, c):
        crossproduct = (c[1] - a[1]) * (b[0] - a[0]) - (c[0] - a[0]) * (b[1] - a[1])
        if abs(crossproduct) > 0.0001:  # sys.float_info.epsilon:
            return False    # (or != 0 if using integers)
        dotproduct = (c[0] - a[0]) * (b[0] - a[0]) + (c[1] - a[1]) * (b[1] - a[1])
        if dotproduct < 0:
            return False
        squaredlengthba = (b[0] - a[0]) * (b[0] - a[0]) + (b[1] - a[1]) * (b[1] - a[1])
        if dotproduct > squaredlengthba:
            return False
        return True

    def Distance(self, y1, x1, y2, x2):
# find the differences between the coordinates
        dy = y2 - y1
        dx = x2 - x1
        ra13 = pow(sin(dy / 2.), 2) + cos(y1) * cos(y2) * pow(sin(dx / 2.), 2)
        return 2 * asin(min(1, sqrt(ra13)))

    def actualDistance(self, y1d, x1d, y2d, x2d):
        x1 = radians(x1d)
        y1 = radians(y1d)
        x2 = radians(x2d)
        y2 = radians(y2d)
        dst = self.Distance(y1, x1, y2, x2)
        return round(abs(dst) * RADIUS, 2)

    def __init__(self, lines, source, target, grid):
        self.source = source
        self.target = target
        self.lines = lines
        self.grid = grid # existing grid lines count
        self.edges = {}
        self.g = Graph()
        for li in range(len(self.lines)):
            vert1 = str(self.lines[li].coordinates[0])
            for pt in range(1, len(self.lines[li].coordinates)):
                vert2 = str(self.lines[li].coordinates[pt])
                dist = self.actualDistance(self.lines[li].coordinates[pt][0], self.lines[li].coordinates[pt][1],
                       self.lines[li].coordinates[pt - 1][0], self.lines[li].coordinates[pt - 1][1])
                self.g.add_edge(vert1, vert2, dist)
                self.edges[vert1 + vert2] = li
                vert1 = vert2
        for li in range(self.grid, len(self.lines)):
            for l2 in range(len(self.lines)):
                if li == l2:
                    continue
                for pt in range(1, len(self.lines[l2].coordinates)):
                    if self.lines[l2].coordinates[pt - 1] == self.lines[li].coordinates[-1] or \
                      self.lines[l2].coordinates[pt] == self.lines[li].coordinates[-1]:
                        continue
                    if self.isBetween(self.lines[l2].coordinates[pt - 1], self.lines[l2].coordinates[pt],
                      self.lines[li].coordinates[-1]):
                        dist = self.actualDistance(self.lines[l2].coordinates[pt - 1][0], self.lines[l2].coordinates[pt - 1][1],
                               self.lines[li].coordinates[-1][0], self.lines[li].coordinates[-1][1])
                        vert1 = str(self.lines[l2].coordinates[pt - 1])
                        vert2 = str(self.lines[li].coordinates[-1])
                        self.g.add_edge(vert1, vert2, dist)
                        self.edges[vert1 + vert2] = l2
                        dist = self.actualDistance(self.lines[l2].coordinates[pt][0], self.lines[l2].coordinates[pt][1],
                               self.lines[li].coordinates[-1][0], self.lines[li].coordinates[-1][1])
                        vert1 = str(self.lines[l2].coordinates[pt])
                        vert2 = str(self.lines[li].coordinates[-1])
                        self.g.add_edge(vert1, vert2, dist)
                        self.edges[vert1 + vert2] = l2
                        break
     #    print 'Graph data:'
      #   print 'source', self.source, 'target', self.target
       #  for v in self.g:
        #     for w in v.get_connections():
         #        vid = v.get_id()
          #       wid = w.get_id()
           #      print '(%s ,%s, %s)'  % (vid, wid, str(v.get_weight(w)))
        if self.g.get_vertex(str(self.source)) is None:
            self.g.add_vertex(str(self.source))
        dijkstra(self.g, self.g.get_vertex(str(self.source)))
        target = self.g.get_vertex(str(self.target))
        try:
            self.path = [target.get_id()]
            shortest(target, self.path)
        except:
            self.path = []

    def getPath(self):
        the_path = []
        for i in range(len(self.path)):
            it = self.path[i].replace('[', '')
            it = it.replace(']', '')
            it = it.replace(' ', '')
            bit = it.split(',')
            the_path.append([float(bit[0]), float(bit[1])])
        return the_path

    def getLines(self):
        the_lines = []
        for i in range(1, len(self.path)):
            try:
                the_lines.append(self.edges[self.path[i - 1] + self.path[i]])
            except:
                try:
                    the_lines.append(self.edges[self.path[i] + self.path[i - 1]])
                except:
                    pass
        the_lines = list(set(the_lines))
        return the_lines
import csv
import datetime
import os
import sys
import xlrd
import configparser   # decode .ini file
from math import sin, cos, pi, sqrt, degrees, radians, asin, atan2
from PyQt5 import QtCore, QtGui, QtWidgets
from mpl_toolkits.basemap import pyproj as pyproj
from towns import Towns
from getmodels import getModelFile
from grid import Grid, Grid_Boundary, Line
from parents import getParents
from senuser import getUser, techClean
from station import Station, Stations
from dijkstra_4 import Shortest
            self._fcapacityGroup.addToGroup(el)
        else:
            self._capacityGroup.addToGroup(el)
        self.addLine(st)
        txt = QtWidgets.QGraphicsSimpleTextItem(st.name)
        new_font = txt.font()
        new_font.setPointSizeF(self.width() * self.txt_ratio)
        txt.setFont(new_font)
        txt.setPos(p + QtCore.QPointF(1.5, -0.5))
        if st.technology[:6] == 'Fossil':
            txt.setBrush(QtGui.QColor(self.colors['fossil_name']))
        else:
            txt.setBrush(QtGui.QColor(self.colors['station_name']))
        txt.setZValue(2)
        self.addItem(txt)
        self._stationGroups[st.name].append(txt)
        self._stationLabels.append(txt)
        if st.technology[:6] == 'Fossil':
            self._fnameGroup.addToGroup(txt)
        else:
            self._nameGroup.addToGroup(txt)
        return

    def addGeneration(self, st):
        if (st.technology[:6] != 'Fossil' or self.show_fossil) and st.generation > 0:
            p = self.mapFromLonLat(QtCore.QPointF(st.lon, st.lat))
            size = sqrt(float(st.generation) * (self.capacity_area / 2000) / pi)
            east = self.destinationxy(st.lon, st.lat, 90., size)
            pe = self.mapFromLonLat(QtCore.QPointF(east[0], st.lat))
            north = self.destinationxy(st.lon, st.lat, 0., size)
            pn = self.mapFromLonLat(QtCore.QPointF(st.lon, north[1]))
            x_d = p.x() - pe.x()
            y_d = pn.y() - p.y()
            el = QtWidgets.QGraphicsEllipseItem(p.x() - x_d / 2, p.y() - y_d / 2, x_d, y_d)
            if self.show_capacity_fill:
                brush = QtGui.QBrush()
                brush.setColor(QtGui.QColor(self.colors[st.technology]))
            #     brush.setStyle(QtCore.Qt.Dense1Pattern)
                brush.setStyle(QtCore.Qt.SolidPattern)
                el.setBrush(brush)
                el.setOpacity(self.capacity_opacity)
            if self.colors['border'] != '':
                el.setPen(QtGui.QColor(self.colors['border']))
            else:
                el.setPen(QtGui.QColor(self.colors[st.technology]))
            el.setZValue(1)
            self.addItem(el)
            self._stationGroups[st.name].append(el)
            self._stationCircles[st.technology].append(el)
            self._generationGroup.addToGroup(el)

    def addLine(self, st):
        if st.technology == 'Rooftop PV':
            return
        if self.new_grid and (st.scenario != 'Existing' or self.trace_existing):
            dims = [[st.lat, st.lon]]
            grid_len = 0.
            grid_path_len = 0.
            gridn = ''
            if st.grid_line is not None:
                gridn = st.grid_line.replace('(', '')
                gridn = gridn.replace('))', '')
                gridn = gridn.replace('),', ',')
                gridn = gridn.replace(')', '')
                gridn = gridn.split(',')
                for j in range(2, len(gridn), 2):
                    dims.append([float(gridn[j]), float(gridn[j + 1])])
                    grid_len += self.lines.actualDistance(dims[-2][0], dims[-2][1], dims[-1][0], dims[-1][1])
            grid_point = self.lines.gridConnect(dims[-1][0], dims[-1][1])
            if grid_point[0] == -1:
                if self.load_centre is not None:   # start with a load centre
                    nearest = 99999
                    j = -1
                    for i in range(len(self.load_centre)):
                        thisone = self.lines.actualDistance(self.load_centre[i][1], self.load_centre[i][2],
                                  dims[-1][0], dims[-1][1])
                        if thisone < nearest:
                            nearest = thisone
                            j = i
                    grid_point[0] = nearest
                    grid_point[1] = self.load_centre[j][1]
                    grid_point[2] = self.load_centre[j][2]
                    grid_point[3] = -1
                else:
                    try:
                        nearest, dist = self._stations.Nearest(dims[-1][0], dims[-1][1], distance=True, fossil=True,
                                        ignore=st.name)
                    except:
                        return
                    grid_point[0] = dist
                    grid_point[1] = nearest.lat
                    grid_point[2] = nearest.lon
                    grid_point[3] = -1
            if grid_point[0] > 0 or len(gridn) > 1:
                grid_len += grid_point[0]
                dims.append([grid_point[1], grid_point[2]])
                if st.technology in self.dispatchable:
                    dispatchable = 'Y'
                else:
                    dispatchable = ''
                l = Line(st.name, self.colors['new_grid'], dims, round(grid_len, 1), grid_point[3], dispatchable, 0,
                    0.)
            #     l.peak_loss = round(st.capacity * self.line_loss * l.length + self.subs_loss * st.capacity, 3)
              #   subs_cost = self.subs_cost * st.capacity
                if l.dispatchable == 'Y':
                    cost, l.line_table = self.lines.Line_Cost(st.capacity, st.capacity)
                else:
                    cost, l.line_table = self.lines.Line_Cost(st.capacity, 0.)
                l.line_cost = round(l.length * cost, 0)
                next = grid_point[3]
                self.lines.lines.append(l)
                if self.trace_grid and self.load_centre is not None:
                    nearest = 99999
                    j = -1
                    for i in range(len(self.load_centre)):
                        thisone = self.lines.actualDistance(self.load_centre[i][1], self.load_centre[i][2],
                                  dims[0][0], dims[0][1])
                        if thisone < nearest:
                            nearest = thisone
                            j = i

                    path = Shortest(self.lines.lines, dims[0], [self.load_centre[j][1],