pysoleng/tests/test_utils.py

#******************************************************************************
#******************************************************************************

# import libraries

import math

from src import pysoleng as pse

# test using "!python -m pytest -s --cov --cov-report term-missing"

#******************************************************************************
#******************************************************************************

def tests():
    
    example_jrc_tmy()
    
    example_jrc_tmy(assume_year=2012)

#******************************************************************************
#******************************************************************************

def example_jrc_tmy(assume_year: int = None):
    
    filename = 'tests/TMY/tmy_54.027_9.803_2005_2020.csv'

    (npdata, 
     names, 
     latitude, 
     longitude, 
     elevation) = pse.utils.read_jrc_tmy_csv(
         filename=filename,
         assume_year=assume_year
         )

    # try to use the data

    site = pse.timeplace.Location(
        coordinates=(latitude, longitude), 
        crs_string='epsg:4326'
        )
    
    #**************************************************************************
    #**************************************************************************
    
    # tracking surface

    site_ground_reflectance = 0.6

    list_I_tilted_total = [] 
    list_I_tilted_beam = [] 
    list_I_tilted_sky_diffuse = [] 
    list_I_tilted_albedo = []

    for i, local_time in enumerate(npdata['timeUTC']):
        
        site_horizontal_total = npdata['Gh'][i]
        site_horizontal_diffuse = npdata['Gdh'][i]
        
        # find solar altitude and solar azimuth
        
        day_year = pse.solar.day_of_the_year(local_time, False)
        
        declination_angle = pse.solar.declination_angle_spencer1971(day_year)
        
        ast = pse.solar.lst_to_ast(local_time, site)
        
        solar_ha = pse.solar.hour_angle(ast)
        
        cos_zen_angle = pse.solar.cosine_zenith_angle(
            latitude*math.pi/180, 
            declination_angle*math.pi/180, 
            solar_ha)
        
        solar_az = pse.solar.solar_azimuth_angle(latitude*math.pi/180, 
                                                 declination_angle*math.pi/180, 
                                                 math.acos(cos_zen_angle), 
                                                 solar_ha)
        
        solar_alt = pse.solar.solar_altitude_angle(latitude*math.pi/180, 
                                                   declination_angle*math.pi/180, 
                                                   solar_ha)
        
        # define the new surface
        
        surface = pse.timeplace.FlatSurface(
            azimuth=solar_az, 
            slope=math.pi/2-solar_alt
            )
        
        # compute the radiation on the tilted plane
        
        (I_tilted_total, 
         I_tilted_beam, 
         I_tilted_sky_diffuse, 
         I_tilted_albedo) = pse.models.tilted_plane_irradiance_isotropic(
             site, 
             surface.slope,
             surface.azimuth,
             local_time, 
             site_ground_reflectance, 
             site_horizontal_total, 
             site_horizontal_diffuse)
             
        list_I_tilted_total.append(I_tilted_total)
        
        list_I_tilted_beam.append(I_tilted_beam)
     
        list_I_tilted_sky_diffuse.append(I_tilted_sky_diffuse)
     
        list_I_tilted_albedo.append(I_tilted_albedo) 
             
    # compare the results

    dni_pysoleng = sum(list_I_tilted_beam)*3600*2.7777778E-07
    dni_jrcpvgis = sum(npdata['Gbn'])*3600*2.7777778E-07

    # print('Results for DNI:')
    # print('pysoleng (isotropic): '+str(dni_pysoleng)+' kWh/m2/y')
    # print('jrc-pvgis (?): '+str(dni_jrcpvgis)+' kWh/m2/y')
    # print('The error between them amounts to ' +
    #       str(abs(dni_pysoleng-dni_jrcpvgis))+' kWh/m2/y.')
    
    assert math.isclose(dni_pysoleng,
                        dni_jrcpvgis,
                        abs_tol=30)
    
    #**************************************************************************
    #**************************************************************************
    
    # fixed plane
        
    surface = pse.timeplace.FlatSurface(
        azimuth=0*math.pi/180, 
        slope=45*math.pi/180
        )
    
    site_ground_reflectance = 0.6
    
    list_I_tilted_total = [] 
    list_I_tilted_beam = [] 
    list_I_tilted_sky_diffuse = [] 
    list_I_tilted_albedo = []
    
    for i, local_time in enumerate(npdata['timeUTC']):
        
        site_horizontal_total = npdata['Gh'][i]
        site_horizontal_diffuse = npdata['Gdh'][i]
        
        #local_time = datetime.fromisoformat(time_interval)
        
        (I_tilted_total, 
          I_tilted_beam, 
          I_tilted_sky_diffuse, 
          I_tilted_albedo) = pse.models.tilted_plane_irradiance_isotropic(
              site, 
              surface.slope,
              surface.azimuth,
              local_time, 
              site_ground_reflectance, 
              site_horizontal_total, 
              site_horizontal_diffuse)
             
        list_I_tilted_total.append(I_tilted_total)
        
        list_I_tilted_beam.append(I_tilted_beam)
     
        list_I_tilted_sky_diffuse.append(I_tilted_sky_diffuse)
     
        list_I_tilted_albedo.append(I_tilted_albedo) 
             
    # sum the components
    
    true_I_tilted_total = 1136.397502028676
    
    abs_tol = 10
    
    assert math.isclose(sum(list_I_tilted_total)*3600*2.7777778E-07, 
                        true_I_tilted_total,
                        abs_tol=abs_tol)
    
    #**************************************************************************
    #**************************************************************************

#******************************************************************************
#******************************************************************************

# def test_read_tmy():
    
#     filename = 'tests/TMY/tmy_54.027_9.803_2005_2020.csv'
    
#     (npdata, 
#      names, 
#      latitude, 
#      longitude, 
#      elevation) = pse.utils.read_jrc_tmy_csv(
#          filename=filename)
    
#     # try to use the data
    
#     site = pse.timeplace.Location(
#         coordinates=(latitude, longitude), 
#         crs_string='epsg:4326')
    
#     surface = pse.timeplace.Plane(
#         azimuth_zero_south_positive_west=0*math.pi/180, 
#         slope_zero_horizontal_positive_south=45*math.pi/180)
    
#     site_ground_reflectance = 0.6
    
#     list_I_tilted_total = [] 
#     list_I_tilted_beam = [] 
#     list_I_tilted_sky_diffuse = [] 
#     list_I_tilted_albedo = []
    
#     for i, local_time in enumerate(npdata['timeUTC']):
        
#         site_horizontal_total = npdata['Gh'][i]
#         site_horizontal_diffuse = npdata['Gdh'][i]
        
#         #local_time = datetime.fromisoformat(time_interval)
        
#         (I_tilted_total, 
#          I_tilted_beam, 
#          I_tilted_sky_diffuse, 
#          I_tilted_albedo) = pse.models.tilted_plane_irradiance_isotropic(
#              site, 
#              surface,
#              local_time, 
#              site_ground_reflectance, 
#              site_horizontal_total, 
#              site_horizontal_diffuse)
             
#         list_I_tilted_total.append(I_tilted_total)
        
#         list_I_tilted_beam.append(I_tilted_beam)
     
#         list_I_tilted_sky_diffuse.append(I_tilted_sky_diffuse)
     
#         list_I_tilted_albedo.append(I_tilted_albedo) 
             
#     # sum the components
    
#     true_I_tilted_total = 1136.397502028676
    
#     abs_tol = 10
    
#     assert math.isclose(sum(list_I_tilted_total)*3600*2.7777778E-07, 
#                         true_I_tilted_total,
#                         abs_tol=abs_tol)

#******************************************************************************
#******************************************************************************

# def test_read_tmy_year_set():
    
#     filename = 'tests/TMY/tmy_54.027_9.803_2005_2020.csv'
    
#     sample_year = 2012
    
#     (npdata, 
#      names, 
#      latitude, 
#      longitude, 
#      elevation) = pse.utils.read_jrc_tmy_csv(
#          filename=filename,
#          assume_year=sample_year)
    
#     # try to use the data
    
#     site = pse.timeplace.Location(
#         coordinates=(latitude, longitude), 
#         crs_string='epsg:4326')
    
#     surface = pse.timeplace.Plane(
#         azimuth_zero_south_positive_west=0*math.pi/180, 
#         slope_zero_horizontal_positive_south=45*math.pi/180)
    
#     site_ground_reflectance = 0.6
    
#     list_I_tilted_total = [] 
#     list_I_tilted_beam = [] 
#     list_I_tilted_sky_diffuse = [] 
#     list_I_tilted_albedo = []
    
#     for i, local_time in enumerate(npdata['timeUTC']):
        
#         site_horizontal_total = npdata['Gh'][i]
#         site_horizontal_diffuse = npdata['Gdh'][i]
        
#         #local_time = datetime.fromisoformat(time_interval)
        
#         (I_tilted_total, 
#          I_tilted_beam, 
#          I_tilted_sky_diffuse, 
#          I_tilted_albedo) = pse.models.tilted_plane_irradiance_isotropic(
#              site, 
#              surface,
#              local_time, 
#              site_ground_reflectance, 
#              site_horizontal_total, 
#              site_horizontal_diffuse)
             
#         list_I_tilted_total.append(I_tilted_total)
        
#         list_I_tilted_beam.append(I_tilted_beam)
     
#         list_I_tilted_sky_diffuse.append(I_tilted_sky_diffuse)
     
#         list_I_tilted_albedo.append(I_tilted_albedo) 
             
#     # sum the components
    
#     true_I_tilted_total = 1136.397502028676
    
#     abs_tol = 10
    
#     assert math.isclose(sum(list_I_tilted_total)*3600*2.7777778E-07, 
#                         true_I_tilted_total,
#                         abs_tol=abs_tol)

#******************************************************************************
#******************************************************************************

def test_angle_format_conversion():
    
    format_table = {
        0: (True, True),
        1: (True, False),
        2: (False, True),
        3: (False, False)
        }
    
    truth_table = {
        0: [0,    45,  90,  135, 180,  225, 270, 315, 360],
        1: [360, 315, 270,  225, 180,  135,  90,  45,   0],
        2: [0,    45,  90,  135, 180, -135, -90, -45,  0],
        3: [0,   -45, -90, -135, 180,  135,  90,  45,  0]
        }
    
    number_points = len(truth_table[0])
    
    for key in truth_table:
        
        assert len(truth_table[key]) == number_points
    
    # for each format
    
    for in_key, in_value in format_table.items():
    
    # for every other format
    
        for out_key, out_value in format_table.items():
            
            if out_key == in_key:
                
                continue
    
            # for each angle
            
            for point_index in range(number_points):
                
                # convert from the first into the second formats
                
                angle_degrees = truth_table[in_key][point_index]
                
                new_angle_degrees = pse.utils.convert_angle(
                    angle_degrees, 
                    in_angle_positive_only=in_value[0], 
                    in_angle_positive_anticlockwise=in_value[1], 
                    out_angle_positive_only=out_value[0], 
                    out_angle_positive_anticlockwise=out_value[1]
                    )
            
                # assert that it adds up
                
                assert math.isclose(
                    math.cos(math.radians(new_angle_degrees)), 
                    math.cos(math.radians(truth_table[out_key][point_index]))
                    )
                
                assert math.isclose(
                    pse.utils.simplify_angle(
                        new_angle_degrees,
                        True), 
                    pse.utils.simplify_angle(
                        truth_table[out_key][point_index],
                        True)
                    )
                
#******************************************************************************
#******************************************************************************

def test_angle_conversion_errors():
    
    # test negative angle
    
    try:
        _ = pse.utils.convert_angle(
            -45, 
            in_angle_positive_only=True, 
            in_angle_positive_anticlockwise=True, 
            out_angle_positive_only=False, 
            out_angle_positive_anticlockwise=False
            )
    except pse.utils.AngleFormatError:
        
        assert True
        
    # test wider than 180 angle
    
    try:
        _ = pse.utils.convert_angle(
            270, 
            in_angle_positive_only=False, 
            in_angle_positive_anticlockwise=True, 
            out_angle_positive_only=False, 
            out_angle_positive_anticlockwise=False
            )
    except pse.utils.AngleFormatError:
        
        assert True
        
#******************************************************************************
#******************************************************************************

def test_angle_simplification():
    
    tests = [
        # degrees
        (90+360, True, 90),
        (3*360, True, 0),
        (270+360, True, 270),
        (-90+360, True, 270),
        (-3*360, True, 0),
        (-270+360, True, 90),
        # radians
        (math.pi/2+2*math.pi, False, math.pi/2),
        (3*2*math.pi, False, 0),
        (math.pi*3/2+2*math.pi, False, math.pi*3/2),
        (-math.pi/2+2*math.pi, False, math.pi*3/2),
        (-3*2*math.pi, False, 0),
        (-math.pi*3/2+2*math.pi, False, math.pi/2),
        ]
    
    for test in tests:
        
        new_angle = pse.utils.simplify_angle(
            test[0],
            test[1])
        
        assert math.isclose(new_angle, test[2])
            
#******************************************************************************
#******************************************************************************