diff --git a/tests.py b/tests.py index d7502bd..2e76e91 100644 --- a/tests.py +++ b/tests.py @@ -1,8 +1,9 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -import unittest import datetime +import unittest + from tmgpm import Tmgpm @@ -12,7 +13,7 @@ def setUp(self): self.tide = Tmgpm() def test_default_station(self): - # check default station name is BREST + # check default station name is BREST self.assertEqual(self.tide.station_name, 'BREST') def test_default_date(self): @@ -27,7 +28,7 @@ def setUp(self): self.tide = Tmgpm('CONCARNEAU', 2014, 1, 1) def test_default_station(self): - # check station name is properly set + # check station name is properly set self.assertEqual(self.tide.station_name, 'CONCARNEAU') def test_default_date(self): @@ -50,24 +51,33 @@ def test_date_after_2100(self): # check default date is today self.assertRaises(ValueError, self.tide.set_date, 2100, 3, 1) + class TestTmgpmGetStationList(unittest.TestCase): def test_returns_a_list(self): # check we get a list station_list = Tmgpm.get_station_list() - self.assertTrue(type(station_list) is list) - + self.assertTrue(isinstance(station_list, list)) + def test_list_is_not_empty(self): # check the list is not empty station_list = Tmgpm.get_station_list() self.assertTrue(len(station_list) > 0) - + def test_list_contains_brest(self): # check the list contains at least BREST station station_list = Tmgpm.get_station_list() self.assertTrue('BREST' in station_list) +class TestTmgpmHarmonicConstituents(unittest.TestCase): + + def test_concarneau_1982(self): + # Check that tide heigh in CONCARNEAU on the 1st January 1982 + # at 00:00 is 2302 mm + tide = Tmgpm('CONCARNEAU', 1982, 1, 1) + self.assertEqual(int(tide.height(0)), 2302) + + if __name__ == '__main__': unittest.main() - diff --git a/tmgpm.py b/tmgpm.py index d028ebc..795aa96 100755 --- a/tmgpm.py +++ b/tmgpm.py @@ -9,74 +9,92 @@ @contact: hleroy@hleroy.com ''' -import os import csv import datetime from math import cos, acos, sin, asin, degrees, radians, floor, sqrt, pow, copysign +import os + + +def sign(x): + return copysign(1, x) -sign = lambda x: copysign(1, x) class Tmgpm(object): - + # Load stations data from csv file stations_data = {} - filename = os.path.join(os.path.abspath(os.path.dirname(__file__)), 'tmgpm.csv') + filename = os.path.join( + os.path.abspath(os.path.dirname(__file__)), 'tmgpm.csv') with open(filename, newline='') as csvfile: - csvreader = csv.DictReader(filter(lambda row: row[0]!='#', csvfile)) + csvreader = csv.DictReader(filter(lambda row: row[0] != '#', csvfile)) for row in csvreader: name = row.pop('NAME') stations_data[name] = row # Constants initialization - n1 = [ [0,0,0,0,0,0,0,0,0,0,0], [0,-2,-3,0,-2,0,0,0,0,0,0], [-2,-3,0,-4,-4,-3,-1,0,0,-2,0], [0,0,0,0,0,0,0,0,0,0,0], [-5,-4,-2,0,0,0,0,0,0,0,0] ] - n2 = [ [1,0,0,0,0,0,0,0,0,0,0], [1,1,1,-1,1,1,0,0,0,0,0], [2,2,0,2,4,4,2,2,-1,2,2], [0,0,0,0,0,0,0,0,0,0,0], [4,4,2,0,0,0,0,0,0,0,0] ] - n3 = [ [0,0,0,0,0,0,0,0,0,0,0], [0,0,1,0,0,0,0,0,0,0,0], [0,1,0,2,0,-1,-1,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0], [1,0,0,0,0,0,0,0,0,0,0] ] - n4 = [ [0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,-1,1,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,-1,1], [0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0] ] - n5 = [ [0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0] ] - n6 = [ [0,0,0,0,0,0,0,0,0,0,0], [1,-1,-1,1,-1,1,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0] ] - + n1 = [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, -2, -3, 0, -2, 0, 0, 0, 0, 0, 0], [-2, -3, 0, -4, -4, -3, -1, 0, 0, -2, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [-5, -4, -2, 0, 0, 0, 0, 0, 0, 0, 0]] + n2 = [[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 1, 1, -1, 1, 1, 0, 0, 0, 0, 0], [2, 2, 0, 2, 4, 4, 2, 2, -1, 2, 2], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0]] + n3 = [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 2, 0, -1, -1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]] + n4 = [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]] + n5 = [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]] + n6 = [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, -1, -1, 1, -1, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]] def __init__(self, station_name=None, year=None, month=None, day=None): """Initialise a Tmpgpm object - + If the station name is omitted, it defaults to 'BREST'. The date must be comprised between 1/3/1900 and 28/02/2100. If the date is omitted, it defaults to today. - + Args: station_name (str, optional)-- station name, defaults to 'BREST' year (int, optional) -- year (4 digits, e.g. 1980), defaults to current year month (int, optional) -- month (1 to 12), defaults to current month day (int, optional) -- day (1 to 31), defaults to current day """ - + + # Initialise instance attributes to None to make PyLint happy + self.station_data = None + self.station_name = None + self.year = None + self.month = None + self.day = None + self.LO = None + self.LA = None + self.UTC = None + self.Z0 = None + self.A = None + self.G = None + self.R0 = None + self.R24 = None + self.PHI0 = None + self.PHI24 = None + # Set station name or defaults to BREST self.set_station(station_name if (station_name is not None) else "BREST") # Set date or defaults to datetime.now() - if not any([year, month, day]): # if year, month and date are all None - now = datetime.datetime.now() # use datetime.now() + if not any([year, month, day]): # if year, month and date are all None + now = datetime.datetime.now() # use datetime.now() year = now.year month = now.month day = now.day self.set_date(year, month, day) - @staticmethod def get_station_list(): """Return the list of stations - + Returns: A list of strings with the station names """ return list(Tmgpm.stations_data.keys()) - def set_station(self, station_name): """Set station name and initialize harmonic data for this station - + Args: station_name (str): station name - + Raises: ValueError: if the station name is not recognized """ @@ -86,18 +104,17 @@ def set_station(self, station_name): self._init_harmonic_data() else: raise ValueError("Station name not recognized") - - + def set_date(self, year, month, day): """Set date and compute data for this date - + The date must be comprised between 1/3/1900 and 28/02/2100. - + Args: year (int) -- Year (4 digits, e.g. 1980) month (int) -- Month (1 to 12) day (int) -- Day (1 to 31) - + Raises: ValueError: if the date is not a valid or out of bounds """ @@ -106,7 +123,7 @@ def set_date(self, year, month, day): d = datetime.date(year, month, day) except: raise ValueError("The date is not valid") - + # Check if the date is comprised between 1/3/1900 and 28/02/2100 if datetime.date(1900, 3, 1) < d < datetime.date(2100, 2, 28): self.year = year @@ -115,8 +132,7 @@ def set_date(self, year, month, day): self._init_precalc() else: raise ValueError("The date is out of bounds. It must be comprised between 1/3/1900 and 28/02/2100") - - + def height(self, t): """Calculate tide height at the given time @@ -127,8 +143,9 @@ def height(self, t): Tide height in mmillimeters """ # Initialize tide height at the station Z0 value - # Z0 is stored in cm and the tide calculations are done in millimeters, hence the * 10 - height = self.Z0 * 10; + # Z0 is stored in cm and the tide calculations are done in millimeters, + # hence the * 10 + height = self.Z0 * 10 # Calculate tide height for j in range(5): @@ -139,35 +156,34 @@ def height(self, t): DELTAj += 360 if DELTAj > 180: DELTAj -= 360 - PHIj = self.PHI0[j] + t / 24.0 * (j * 360 + DELTAj); - height += Rj * cos(radians(PHIj)); + PHIj = self.PHI0[j] + t / 24.0 * (j * 360 + DELTAj) + height += Rj * cos(radians(PHIj)) - # and return the value - return height; + # and return the value + return height - def _init_harmonic_data(self): - self.A = [[0. for x in range(11)] for x in range(5)] - self.G = [[0. for x in range(11)] for x in range(5)] + self.A = [[0. for _ in range(11)] for _ in range(5)] + self.G = [[0. for _ in range(11)] for _ in range(5)] self.Z0 = float(self.station_data['Z0']) self.A[0][0] = float(self.station_data['ASa']) self.A[1][0] = float(self.station_data['AK1']) self.A[1][1] = float(self.station_data['AO1']) self.A[1][2] = float(self.station_data['AQ1']) - self.A[1][3] = -1/3.0 * float(self.station_data['AK1']) - self.A[1][4] = 1/5.3 * float(self.station_data['AO1']) - self.A[1][5] = 1/7.4 * float(self.station_data['AK1']) + self.A[1][3] = -1 / 3.0 * float(self.station_data['AK1']) + self.A[1][4] = 1 / 5.3 * float(self.station_data['AO1']) + self.A[1][5] = 1 / 7.4 * float(self.station_data['AK1']) self.A[2][0] = float(self.station_data['AM2']) self.A[2][1] = float(self.station_data['AN2']) self.A[2][2] = float(self.station_data['AS2']) - self.A[2][3] = 1/7.6 * float(self.station_data['AN2']) - self.A[2][4] = 1/6.3 * float(self.station_data['AN2']) - self.A[2][5] = 1/5.3 * float(self.station_data['AN2']) - self.A[2][6] = -1/35.0 * float(self.station_data['AM2']) - self.A[2][7] = 1/3.7 * float(self.station_data['AS2']) - self.A[2][8] = 1/17.0 * float(self.station_data['AS2']) - self.A[2][9] = -1/27.0 * float(self.station_data['AM2']) - self.A[2][10] = 1/12.0 * float(self.station_data['AS2']) + self.A[2][3] = 1 / 7.6 * float(self.station_data['AN2']) + self.A[2][4] = 1 / 6.3 * float(self.station_data['AN2']) + self.A[2][5] = 1 / 5.3 * float(self.station_data['AN2']) + self.A[2][6] = -1 / 35.0 * float(self.station_data['AM2']) + self.A[2][7] = 1 / 3.7 * float(self.station_data['AS2']) + self.A[2][8] = 1 / 17.0 * float(self.station_data['AS2']) + self.A[2][9] = -1 / 27.0 * float(self.station_data['AM2']) + self.A[2][10] = 1 / 12.0 * float(self.station_data['AS2']) self.A[4][0] = float(self.station_data['AMN4']) self.A[4][1] = float(self.station_data['AM4']) self.A[4][2] = float(self.station_data['AMS4']) @@ -192,61 +208,58 @@ def _init_harmonic_data(self): self.G[4][0] = float(self.station_data['GMN4']) self.G[4][1] = float(self.station_data['GM4']) self.G[4][2] = float(self.station_data['GMS4']) - self.LA = self.station_data['LA'] - self.LO = self.station_data['LO'] - self.UTC = self.station_data['UTC'] - #TODO: add control mechanism - + self.LA = self.station_data['LA'] + self.LO = self.station_data['LO'] + self.UTC = self.station_data['UTC'] + # TODO: add control mechanism def _init_precalc(self): def compute_R_and_PHI_at(t): R = [0., 0., 0., 0., 0.] PHI = [0., 0., 0., 0., 0.] - - T = floor(30.6001 * (1 + self.month + 12 * floor(1 / (self.month + 1.0) + 0.7))) + floor(365.25 * (self.year - floor(1 / (self.month + 1.0) + 0.7))) + self.day + t / 24 - 723258; - h = 279.82 + 0.98564734 * T; - s = 78.16 + 13.17639673 * T; - p = 349.5 + 0.11140408 * T; - N = 208.1 + 0.05295392 * T; - p1 = 282.6 + 0.000047069 * T; - D = 90; - + + T = floor(30.6001 * (1 + self.month + 12 * floor(1 / (self.month + 1.0) + 0.7))) + floor( + 365.25 * (self.year - floor(1 / (self.month + 1.0) + 0.7))) + self.day + t / 24 - 723258 + h = 279.82 + 0.98564734 * T + s = 78.16 + 13.17639673 * T + p = 349.5 + 0.11140408 * T + N = 208.1 + 0.05295392 * T + p1 = 282.6 + 0.000047069 * T + D = 90 + for j in range(5): if j != 3: - Xj = 0; - Yj = 0; + Xj = 0 + Yj = 0 for i in range(11): - Vij = 15 * j * t + Tmgpm.n1[j][i] * s + Tmgpm.n2[j][i] * h + Tmgpm.n3[j][i] * p + Tmgpm.n4[j][i] * N + Tmgpm.n5[j][i] * p1 + Tmgpm.n6[j][i] * D; - Xj += self.A[j][i] * cos(radians(Vij - self.G[j][i])); - Yj += self.A[j][i] * sin(radians(Vij - self.G[j][i])); - - R[j] = sqrt(pow(Xj, 2) + pow(Yj, 2)); + Vij = 15 * j * t + Tmgpm.n1[j][i] * s + Tmgpm.n2[j][i] * h + Tmgpm.n3[j][i] * p + \ + Tmgpm.n4[j][i] * N + Tmgpm.n5[j][i] * p1 + Tmgpm.n6[j][i] * D + Xj += self.A[j][i] * cos(radians(Vij - self.G[j][i])) + Yj += self.A[j][i] * sin(radians(Vij - self.G[j][i])) + + R[j] = sqrt(pow(Xj, 2) + pow(Yj, 2)) if R[j] == 0: - PHI[j] = 90; + PHI[j] = 90 else: - PHI[j] = degrees(acos(Xj / R[j])) * sign(degrees(asin(Yj / R[j]))); - + PHI[j] = degrees( + acos(Xj / R[j])) * sign(degrees(asin(Yj / R[j]))) + return R, PHI self.R0, self.PHI0 = compute_R_and_PHI_at(0) self.R24, self.PHI24 = compute_R_and_PHI_at(24) - def __repr__(self): """Return a string containing a printable representation of a Tmgpm object""" return """'{}.{}("{}", {}, {}, {})'""".format(self.__class__.__module__, - self.__class__.__qualname__, + self.__class__.__name__, self.station_name, self.year, self.month, self.day) - def __str__(self): """Return a string containing a nicely printable representation of a Tmgpm""" return """Tide object for {} on {:%d, %b %Y}""".format(self.station_name, - datetime.date(self.year, - self.month, - self.day)) - + datetime.date(self.year, self.month, self.day))