prep.py

#
# Title: Data Preparation for LINK Configs and Transmissions
# Author: Fabian Kovac <fabian.kovac@fhstp.ac.at>
# Team: University of Applied Sciences St. Pölten
# Version: 4.1
# Last changed: 2022-01-21
#

import sys
import gzip
import json
import pathlib
import argparse
import datetime

import numpy as np
import pandas as pd

from shapely.geometry import Point, Polygon


def parse_arguments() -> argparse.Namespace:
    """Parses provided commandline arguments for LINK
	
	Returns:
		args (argparse.Namespace): object with paths to provided files
	"""
    
    # create argument parser with description
    desc = '# Data Preparation for LINK\n'
    desc += '-'*64 + '\n'
    desc += 'Script outputs the same files with a "_clean" suffix.\n'
    desc += 'Existing clean versions are automatically overwritten!'
    
    parser = argparse.ArgumentParser(
        prog = 'prep.py',
        usage = 'python %(prog)s -c <config_file> -t <transmissions_file> -i <inca_dir>',
        description = desc,
        formatter_class = argparse.RawTextHelpFormatter
    )
    
    # add required argument group
    # add config parameter to parser
    # add transmissions parameter to parser
    required_args = parser.add_argument_group('required arguments')
    required_args.add_argument('-c', '--config', type = str, required = True, help = 'Path to Config-File')
    required_args.add_argument('-t', '--transmissions', type = str, required = True, help = 'Path to Transmissions-File')
    required_args.add_argument('-i', '--inca', type = str, required = True, help = 'Path to Inca-Dir')
    
    # parse arguments
    args = parser.parse_args()
    
    return args


def _log(msg: str) -> None:
    """Logs message
    
    Parameters:
        msg (str): Message to log to console
	"""
    
    # add marker to log message
    marker = '%'
    if msg[:1] == '\n':
        msg = f'\n{marker} {msg[1:]}'
    else:
        msg = f'{marker} {msg}'
    
    # print message
    print(msg)


def load_areas(file: pathlib.Path = 'areas.json') -> dict:
    """Load json file with areas and corresponding polygons
    
    Parameters:
       file (pathlib.Path): Path to areas json file
    
    Returns:
        areas_poly (dict): Dict containing areas and corresponding polygons
    """
    
    # area dictionary
    areas_poly = {}

    # load areas geojson file
    with open(file) as f:   
        for area, points in json.load(f).items():        
            transformed = Polygon([tuple(p) for p in points])
            areas_poly.update({area : transformed})
    
    return areas_poly


def load_inca_file(file_inca: pathlib.Path) -> np.array:
    """Loads gzipped INCA data of a given file
    
    Parameters:
       file_inca (pathlib.Path): Path to INCA file
    
    Returns:
        x (np.array): Matrix with shape (401, 701)
    """
    
    # open zipped file and bring data to right shape (resolution x,y: 701x401 km2)
    with gzip.open(file_inca, 'rb') as file:
        x = file.read()
        x = np.fromstring(x, dtype = np.float32, sep = ' ')
        x = np.reshape(x, (401, 701))

        return x


def load_inca_data(dir_inca: pathlib.Path) -> np.array:
    """Loads inca files of a given inca dir
    
    Parameters:
        dir_inca (pathlib.Path): Directory for INCA files
    
    Returns:
        inca_Data (np.array): Tensor with shape (96, 401, 701)
    """
    
    # initialize tensor with 96 15min intervals
    inca_data = np.zeros((96, 401, 701))
    
    # load inca dates from inca dir
    for i, file_inca in enumerate(sorted([file for file in dir_inca.iterdir() if file.is_file()])):        
        # load zipped ascii data
        data = load_inca_file(file_inca)
        
        # update inca tensor
        inca_data[i] = data
    
    return inca_data


def get_midpoint(lon_a: np.array, lat_a: np.array, lon_b: np.array, lat_b: np.array) -> np.array:
    """Calculcates the midpoint between two wgs 84 coordinate pairs
    
    Parameters:
        lon_a (np.array): Longitudes of point A
        lat_a (np.array): Latitudes of point A
        lon_b (np.array): Longitudes of point B
        lat_b (np.array): Latitudes of point B
    
    Returns:
        lon, lat (np.array): Vectors with longitudes and latitudes of the midpoints
    """
    
    # convert degrees to radians
    lon_a, lat_a = np.radians(lon_a), np.radians(lat_a)
    lon_b, lat_b = np.radians(lon_b), np.radians(lat_b)
    
    # calculate midpoint
    Bx = np.cos(lat_b) * np.cos(lon_b - lon_a)
    By = np.cos(lat_b) * np.sin(lon_b - lon_a)
    
    lon_mid = lon_a + np.arctan2(By, np.cos(lat_a) + Bx)
    lat_mid = np.arctan2(
        np.sin(lat_a) + np.sin(lat_b),
        np.sqrt(np.square(np.cos(lat_a) + Bx) + np.square(By))
    )
    
    return np.degrees(lon_mid), np.degrees(lat_mid)


def get_distance(lon_a: np.array, lat_a: np.array, alt_a: np.array, lon_b: np.array, lat_b: np.array, alt_b: np.array) -> np.array:
    """Calculcates distance between two coordinates in km using a WGS84 rotation-ellipsoid
    
    Parameters:
        lon_a (np.array): Longitudes of point A
        lat_a (np.array): Latitudes of point A
        alt_a (np.array): Altitudes of point A
        lon_b (np.array): Longitudes of point B
        lat_b (np.array): Latitudes of point B
        alt_b (np.array): Altitudes of point B
    
    Returns:
        length (np.array): Vector with distances in km
    """
    
    # constants (equator radius and pole radius in km)
    # r_equator is the 'semi-major axis' and r_pole the 'semi-minor axis' on a WGS84 ellipsoid
    r_equator = 6378.137
    r_pole = 6356.7523142
    
    # calculate points on rotation-ellipsoid
    # point A
    za = np.sin(np.radians(lat_a)) * r_pole
    ra = np.cos(np.radians(lat_a)) * r_equator
    xa = np.sin(np.radians(lon_a)) * ra
    ya = np.cos(np.radians(lon_a)) * ra
    
    # point B
    zb = np.sin(np.radians(lat_b)) * r_pole
    rb = np.cos(np.radians(lat_b)) * r_equator
    xb = np.sin(np.radians(lon_b)) * rb
    yb = np.cos(np.radians(lon_b)) * rb
    
    # calculate differences between A and B
    dx = xa - xb
    dy = ya - yb
    dz = za - zb
    
    # calculate distances between points using the square root of the sum of squares
    distance = np.sqrt(np.square(dx) + np.square(dy) + np.square(dz))
    
    # incorporate altitudes into the distance
    distance = np.sqrt(np.square(distance) + np.square((alt_a/1e3) - (alt_b/1e3)))
    	
    return distance


def get_link_area(lon: float, lat: float) -> str:
    """Get geographical area depending on the start of the link
    
    Parameters:
        lon (float): Float containing wgs 84 longitudes
        lat (float): Float containing wgs 84 latitudes
    
    Returns:
        area (str): Geographical area the link belongs to
	"""    
    
    # define point with start-coordinates of link
    point = Point(lon, lat)
    
    # get area depending on point
    for area, poly in areas_poly.items():
        if point.within(poly):
            return area
    
    return 'ukn'


def wgs84_to_lambert(lon: np.array, lat: np.array) -> tuple:
    """Convert UTM coordinates to an Austrian Lambert Conic Conformal Projection (lcc)
    EPSG code: 31287, see https://epsg.io/31287
    
    Parameters:
        lon (np.array): Vector containing wgs 84 longitudes
        lat (np.array): Vector containing wgs 84 latitudes
    
    Returns:
        x, y (np.array): Vectors containing lcc coordinates
	"""
    
    # convert wgs 84 coordinates to radians
    lon = np.radians(lon)
    lat = np.radians(lat)

    # define standard parallels according to EPSG:31287 - MGI/Austria Lambert
    lat1 = np.radians(49)
    lat2 = np.radians(46)
    
    # point of reference:
    #   lon: 13°20'E
    #   lat: 47°30'N
    lon0 = np.radians(13.33333333333333)
    lat0 = np.radians(47.5)
    
    # false easting/northing
    # INCA grid: 701x401 km
    # compensate for point of reference (middle of the grid)
    # move point of reference (0, 0) to following x and y coordinates (x0, y0) to e.g. get strictly positive coordinates
    # false easting: 3505000m
    # false northing: 200500m
    x0 = 350500
    y0 = 200500
    
    # volumetric mean radius of the earth in m
    R = 6371000

    # lambert conformal conic projection as in https://pubs.usgs.gov/pp/1395/report.pdf
    # reference:
    # >>    J. P. Snyder, Map projections--a working manual. Washington: U.S. G.P.O.,
    # >>    For sale by the Supt. of Docs., 1987., pp. 104-110
    n = np.log(np.cos(lat1) * (1 / np.cos(lat2))) / np.log(np.tan(np.pi/4 + lat2/2) * (1 / np.tan(np.pi/4 + lat1/2)))
    F = (np.cos(lat1) * np.tan(np.pi/4 + lat1/2)**n) / n
    p = R * F * (1 / np.tan(np.pi/4 + lat/2))**n
    p0 = R * F * (1 / np.tan(np.pi/4 + lat0/2))**n

    # calculate x and y in lcc
    x = p * np.sin(n * (lon - lon0)) + x0
    y = p0 - p * np.cos(n * (lon - lon0)) + y0

    return x, y


def lambert_to_inca_idx(x: np.array, y: np.array) -> tuple:
    """Convert x and y of Lambert Conic Conformal Projection to an INCA index
    
    Parameters:
        x (np.array): Vector containing x values in meter of Lambert Conformal Conic Projection
        y (np.array): Vector containing y values in meter of in Lambert Conformal Conic Projection
    
    Returns:
        ix, iy (np.array): indices for INCA data
	"""
    
    # convert from m to km
    x /= 1000
    y /= 1000
    
    # convert coordinates to INCA index (nearest coordinates as int)
    ix = np.round(x, decimals = 0).astype(int)
    iy = np.round(y, decimals = 0).astype(int)
    
    return ix, iy


def get_inca_data(datetimes: np.array, lon_a: np.array, lat_a: np.array, lon_b: np.array, lat_b: np.array, lengths: np.array) -> np.array:
    """Get INCA RR data based for each km of link path
    
    Parameters:
        datetimes (np.array): Vector containing datetimes of transmissions
        lon_a (np.array): Vector containing longitude values of LINK site a
        lat_a (np.array): Vector containing latitude values of LINK site a
        lon_b (np.array): Vector containing longitude values of LINK site b
        lat_b (np.array): Vector containing latiotude values of LINK site b
        length (np.array): Vector containing distance between sites in km
    
    Returns:
        inca_RR (np.array): Vector containing INCA RR data for each transmission
	"""
    
    # load inca data
    inca_data = load_inca_data(dir_inca)
    _log(f'Loaded INCA data from {str(dir_inca).split("/")[-1]} with shape {inca_data.shape}')
    
    # generate times of day in 15min (window) intervals
    window = 15
    inca_times = sorted([str(i * datetime.timedelta(seconds = window))[2:] for i in range(24*60//window)])
    
    # generate times of LINK data (last 4 chars of datetime vector)
    time_func = np.vectorize(lambda x: f'{x[-4:-2]}:{x[-2:]}')
    link_times = time_func(datetimes)
    
    # get INCA indices of LINK times
    idx_times = np.searchsorted(inca_times, link_times)
    _log(f'Created search indices vector for each transmission for INCA RR tensor')
    
    
    # convert wgs 84 coordinates to lambert conic conformal projection
    xa, ya = wgs84_to_lambert(lon_a, lat_a)
    xb, yb = wgs84_to_lambert(lon_b, lat_b)
    
    # create lambert coordinate matrices for start and end points    
    start = np.array([xa, ya]).T
    end = np.array([xb, yb]).T
    
    # define n INCA RR points for each km of link length
    # --> min. points: 3 (start, mid and end) if length < 3km
    # --> get odd amount of RR points to always include the middle of the link
    n = np.round(lengths, decimals = 0).astype(int)
    n = np.where(n < 3, 3, n)
    n = np.where(n%2 == 0, n+1, n)
    
    # set list for INCA RR data
    inca_RR = []
    
    # get n INCA RR data points for each transmission
    # TODO: find vectorized solution without loop
    # --> n is dynamic (diffent lengths of links) and np.linspace() only supports n as a scalar
    for i in range(len(link_times)):        
        # create n points for each km of link
        # --> lambert projection is in meters, therefore linspace equally divides the link into n points (including start, end and mid)
        x, y = np.linspace(start[i], end[i], n[i], axis = 1)
        
        # get nearest inca coordinates for each point on link
        idx, idy = lambert_to_inca_idx(x, y)
        
        # get INCA RR data for each point on link
        inca_RR.append(inca_data[idx_times[i], idy, idx])
    
    return inca_RR


def prep() -> None:
    """Data preparation for LINK config and transmissions"""
    
    _log('\n******************************** READ FILES ********************************')
    
    # read config and transmission files
    df_config = pd.read_csv(file_config, sep = ';')
    _log(f'Read config file with shape {df_config.shape}')
    df_trans = pd.read_csv(file_trans, sep = ';')
    _log(f'Read transmissions file with shape {df_trans.shape}')
    
    
    _log('\n******************************** BASIC PREP ********************************')
    
    # remove test-link with link id 1
    df_config = df_config[df_config['LINKID'] != 1]
    df_trans = df_trans[df_trans['RADIOLINKID'] != 1]
    _log('Removed all entries of test-link with linkid 1')
    
    
    # drop links that are officially not in use ('na' in CAPACITYINTERFACE and/or FREQUENCY)
    # --> see Q&A Phillip Scheffknecht (05 Feb 2021)
    df_config = df_config.dropna(axis = 0, subset = ['CAPACITYINTERFACE', 'FREQUENCY'])
    
    # check if RXFREQUENCY and TXFREQUENCY exists (only present after 2021-05)
    if 'TXFREQUENCY' in df_config.columns and 'RXFREQUENCY' in df_config.columns:
        df_config = df_config.dropna(axis = 0, subset = ['RXFREQUENCY', 'TXFREQUENCY'])
    
    _log('Dropped configs with NA in CAPACITYINTERFACE and/or FREQUENCY (links officially not in use)')
    
    
    # delete rows with unused link ids
    # get link ids of config and transmissions
    config_ids = set(df_config['LINKID'].unique().tolist())
    trans_ids = set(df_trans['RADIOLINKID'].unique().tolist())
    
    # delete link ids in transmissions without config
    unused_trans_ids = trans_ids - config_ids
    df_trans = df_trans[~df_trans['RADIOLINKID'].isin(list(unused_trans_ids))]
    _log('Removed all links in transmissions where no config is present')
    
    # delete link ids in config without transmissions
    unused_config_ids = config_ids - trans_ids
    df_config = df_config[~df_config['LINKID'].isin(list(unused_config_ids))]
    _log('Removed all links in config where no transmission is present')
    
    
    # delete duplicates in config (same values, different link ids), where corresponding link ids are not used in transmissions
    # gather duplicated rows in config file
    col_subset = ['LINKTYPE', 'SITEID_A', 'LATITUDE_A', 'LONGITUDE_A', 'SITEID_B', 'LATITUDE_B', 'LONGITUDE_B', 'CAPACITYINTERFACE', 'FREQUENCY']
    duplicated_config_ids = set(df_config[df_config.duplicated(subset = col_subset)]['LINKID'].unique().tolist())
    
    # gather duplicated link ids of config file in transmissions file
    found_trans_ids = set(df_trans[df_trans['RADIOLINKID'].isin(duplicated_config_ids)]['RADIOLINKID'].unique().tolist())
    
    # calculate unused duplicated ids in config file
    duplicated_used_ids = duplicated_config_ids - found_trans_ids
    
    # delete rows with unused duplicated link ids in config file
    df_config = df_config[~df_config['LINKID'].isin(list(duplicated_used_ids))]
    _log('Removed duplicated links which are not in use')
    
    
    # check (RADIOLINKID - LOCATION) pairs with (LINKID - SITEID_A) or (LINKID - SITEID_B) pairs
    # gather unique combinations of link ids and site ids of config
    # (new temporary column is much faster then pandas agg function)
    df_config['TEMP_LOC_TUPLE'] = df_config['LINKID'].astype(str) + '--' + df_config['SITEID_A']
    config_loc_tuples = set(df_config['TEMP_LOC_TUPLE'].unique().tolist())
    df_config['TEMP_LOC_TUPLE'] = df_config['LINKID'].astype(str) + '--' + df_config['SITEID_B']
    config_loc_tuples.update(set(df_config['TEMP_LOC_TUPLE'].unique().tolist()))
    
    # gather unique combinations of link ids and locations in transmissions
    # (new temporary column is much faster then pandas agg function)
    df_trans['TEMP_LOC_TUPLE'] = df_trans['RADIOLINKID'].astype(str) + '--' + df_trans['LOCATION']
    trans_loc_tuples = set(df_trans['TEMP_LOC_TUPLE'].unique().tolist())
    
    # calculate link id - location tuples which are not in config
    invalid_loc_tuples = trans_loc_tuples - config_loc_tuples
    
    # remove invalid tuples from transmissions
    df_trans = df_trans[~df_trans['TEMP_LOC_TUPLE'].isin(list(invalid_loc_tuples))]
    _log('Removed all transmissions with invalid RADIOLINKID - LOCATION tuples (not present in config)')
    
    # remove temp columns
    df_config = df_config.drop(['TEMP_LOC_TUPLE'], axis = 1)
    df_trans = df_trans.drop(['TEMP_LOC_TUPLE'], axis = 1)
    
    
    # check if ANT_HEIGHT_ALTITUDE_A and ANT_HEIGHT_ALTITUDE_B exists (only present after 2021-09)
    # otherwise get altitudes of config file from "VIW_ZAMG_MW_CONFIG_DY_20210928_SH.csv"
    # --> see mail Jürgen Köhler (juergen.koehler@drei.com) at 2021-09-27 16:15
    if 'ANT_HEIGHT_ALTITUDE_A' not in df_config.columns and 'ANT_HEIGHT_ALTITUDE_B' not in df_config.columns:
        df_altitudes = pd.read_csv('files_from_H3A/CONFIG_ALTITUDES.csv', sep = ';')
        df_altitudes = df_altitudes[['LINKID', 'ANT_HEIGHT_ALTITUDE_A', 'ANT_HEIGHT_ALTITUDE_B']]
        df_config = pd.merge(df_config, df_altitudes, how = 'inner', left_on = 'LINKID', right_on = 'LINKID')
        del df_altitudes
        _log('Added altitude of links')
    
    
    # calculate midpoint of links
    df_config['LONGITUDE_MID'], df_config['LATITUDE_MID'] = get_midpoint(
        df_config['LONGITUDE_A'], df_config['LATITUDE_A'],
        df_config['LONGITUDE_B'], df_config['LATITUDE_B']
    )
    df_config['ANT_HEIGHT_ALTITUDE_MID'] = (df_config['ANT_HEIGHT_ALTITUDE_A'] + df_config['ANT_HEIGHT_ALTITUDE_B']) // 2
    df_config['ANT_HEIGHT_ALTITUDE_MID'] = df_config['ANT_HEIGHT_ALTITUDE_MID'].astype('int')
    _log('Calculated midpoint of links')
    
    # calculate LENGTH in km between links
    df_config['LENGTH'] = get_distance(
        df_config['LONGITUDE_A'],
        df_config['LATITUDE_A'],
        df_config['ANT_HEIGHT_ALTITUDE_A'],
        df_config['LONGITUDE_B'],
        df_config['LATITUDE_B'],
        df_config['ANT_HEIGHT_ALTITUDE_B']
    )
    _log('Calculated distances between sites using a WGS84 ellipsoid')
    
    
    # calculate geographical area of link
    df_config['AREA'] = df_config.apply(lambda row: get_link_area(row['LONGITUDE_A'], row['LATITUDE_A']), axis = 1)
    
    
    # convert CAPACITYINTERFACE to int
    df_config['CAPACITYINTERFACE'] = df_config['CAPACITYINTERFACE'].map(lambda x: str(x)[:-5]).astype('int')
    _log('Converted CAPACITYINTERFACE to int')
    
    # convert FREQUENCY to float
    df_config['FREQUENCY'] = df_config['FREQUENCY'].map(lambda x: str(x)[:-3]).astype('float')
    _log('Converted FREQUENCY to float')
    
    # convert RXFREQUENCY and TXFREQUENCY to float and from MHz to GHz
    # check if columns exists (only present with 2021-05)
    if 'TXFREQUENCY' in df_config.columns and 'RXFREQUENCY' in df_config.columns:
        df_config['RXFREQUENCY'] = df_config['RXFREQUENCY'].astype('float')
        df_config['RXFREQUENCY'] = df_config['RXFREQUENCY']/1e3
        df_config['TXFREQUENCY'] = df_config['TXFREQUENCY'].astype('float')
        df_config['TXFREQUENCY'] = df_config['TXFREQUENCY']/1e3
        _log('Converted RXFREQUENCY and TXFREQUENCY to float and GHz')
    
    
    # drop transmissions with (operational) status unequal 1
    df_trans = df_trans[df_trans['STATUS'] == 1]
    df_trans = df_trans[df_trans['OPERATIONALSTATUS'] == 1]
    _log('Removed transmissions with STATUS and/or OPERATIONALSTATUS unequal 1')
    
    
    _log('\n******************************** BUILD LINK DF *****************************')
    
    # copy transmissions dataframe to link dataframe
    df_link = df_trans.copy()
    _log('Copy transmissions dataframe to link dataframe')
    
    
    # convert begintime to utc
    df_link['BEGINTIME'] = pd.to_datetime(df_link['BEGINTIME'], format = '%Y-%m-%d %H:%M:%S')
    df_link['BEGINTIME'] = df_link['BEGINTIME'].dt.tz_localize('Europe/Vienna').dt.tz_convert('UTC').dt.tz_localize(None)
    _log('Converted BEGINTIME to UTC')
    
    # only use transmissions with begintime matching INCA date
    date_inca = np.datetime64(f'{str(dir_inca.stem)[0:4]}-{str(dir_inca.stem)[-4:-2]}-{str(dir_inca.stem)[-2:]}')
    df_link = df_link[df_link['BEGINTIME'].dt.date == date_inca]
    _log('Filtered transmissions to match BEGINTIME with INCA date')
    
    
    # copy REMOTERXLEVEL to PMIN and PMAX (for aggregation in 15min time window)
    df_link['PMIN'] = df_link['REMOTERXLEVEL']
    df_link['PMAX'] = df_link['REMOTERXLEVEL']
    _log('Created PMIN, PMAX and mean REMOTERXLEVEL')

    # convert 3min windows to 15min windows to get min, max and mean power levels
    group_cols = [df_link['BEGINTIME'].dt.round('15Min'), 'RADIOLINKID']
    agg_cols = {
        'PMIN' : 'min',
        'PMAX' : 'max',
        'REMOTERXLEVEL' : 'mean',
        'TXLEVEL' : 'mean',
        'SPEED' : 'median',
        'CURRRXBITRATE' : 'median',
        'CURRTXBITRATE' : 'median',
        'CURRRXPROFILE' : 'median',
        'CURRTXPROFILE' : 'median'
    }
    df_link = df_link.groupby(group_cols).agg(agg_cols).reset_index()
    _log('Converted 3min windows to 15min windows to get min, max and mean power levels')
    
    
    # convert BEGINTIME to format "%Y%m%d%H%M"
    df_link['BEGINTIME'] = df_link['BEGINTIME'].dt.strftime('%Y%m%d%H%M')
    _log('Converted BEGINTIME to format "%Y%m%d%H%M"')
    
    
    # merge config and link dataframe
    drop_cols = ['RADIOLINKID', 'LINKTYPE', 'SITEID_A', 'SITEID_B', 'CAPACITYINTERFACE']
    df_link = pd.merge(df_link, df_config, how = 'inner', left_on = 'RADIOLINKID', right_on = 'LINKID').drop(drop_cols, axis = 1)
    _log('Merged config data to link dataframe')
    
    
    _log('\n******************************** CALC POWER LEVELS *************************')
    
    # load dry LINK attenuation statistics
    df_stats = pd.read_csv('files_from_H3A/00_Link_dry_attnstats.csv', sep = ';')
    
    # merge mean and std attenuation to transmission dataframe
    df_link = pd.merge(df_link, df_stats, how = 'left', left_on = 'LINKID', right_on = 'LINKID')
    df_link = df_link.dropna(axis = 0, subset = ['ATTENUATIONMEANDRY', 'ATTENUATIONSTDDRY'])
    del df_stats
    _log('Merged mean and std of link attenuations to link dataframe')
    
    
    
    # calculate attenuation of transmitted and received signal strengths per transmission
    df_link['ATTENUATION'] = df_link['REMOTERXLEVEL'] - df_link['TXLEVEL']
    _log('Calculated attenuation per transmission')

    # calculate differences from attenuations and mean attenuations
    df_link['ATTENUATIONDIFF'] = df_link['ATTENUATION'] - df_link['ATTENUATIONMEANDRY']
    
    df_std = df_link.groupby(['LINKID']).agg({'ATTENUATION' : 'std'}).reset_index()
    df_std = df_std.rename(columns = {'LINKID' : 'LINKID', 'ATTENUATION' : 'ATTENUATIONSTD'})
    
    df_link = pd.merge(df_link, df_std, how = 'left', left_on = 'LINKID', right_on = 'LINKID')
    df_link['ATTENUATIONSTDDIFF'] = np.sqrt(df_link['ATTENUATIONSTD']**2 + df_link['ATTENUATIONSTDDRY']**2)
    del(df_std)
    _log('Calculated attenuation differences of dry means')
    
    # calculate extinction coefficients for links based on attenuation and lengths of the links based on beer-lambert law
    # Coeff = -ln(10) * (Attenuation/PathLength)
    df_link['EXTINCTIONCOEFF'] = -np.log(10) * (df_link['ATTENUATION'] / (df_link['LENGTH']))
    _log('Calculated extinction coefficients based on beer-lambert law')
    
    
    # rename and reorder columns
    name_cols = {
        'LINKID' : 'ID',
        'BEGINTIME' : 'DateTime',
        'REMOTERXLEVEL' : 'RxLevel',
        'PMIN' : 'RxLevelMin',
        'PMAX' : 'RxLevelMax',
        'TXLEVEL' : 'TxLevel',
        'ATTENUATION' : 'Attn',
        'ATTENUATIONSTD' : 'AttnStd',
        'ATTENUATIONMEANDRY' : 'AttnMeanDry',
        'ATTENUATIONSTDDRY' : 'AttnStdDry',
        'ATTENUATIONDIFF' : 'AttnDiff',
        'ATTENUATIONSTDDIFF' : 'AttnStdDiff',
        'EXTINCTIONCOEFF' : 'ExtinctionCoeff',
        'SPEED' : 'Speed',
        'CURRRXBITRATE' : 'RxBitrate',
        'CURRTXBITRATE' : 'TxBitrate',
        'CURRRXPROFILE' : 'RxProfile',
        'CURRTXPROFILE' : 'TxProfile',
        'LONGITUDE_A' : 'XStart',
        'LATITUDE_A' : 'YStart',
        'ANT_HEIGHT_ALTITUDE_A' : 'AltStart',
        'LONGITUDE_MID' : 'XMid',
        'LATITUDE_MID' : 'YMid',
        'ANT_HEIGHT_ALTITUDE_MID' : 'AltMid',
        'LONGITUDE_B' : 'XEnd',
        'LATITUDE_B' : 'YEnd',
        'ANT_HEIGHT_ALTITUDE_B' : 'AltEnd',
        'LENGTH' : 'PathLength',
        'FREQUENCY' : 'Frequency',
        'AREA' : 'Area'
    }
    
    # check if RXFREQUENCY and TXFREQUENCY exists (only present after 2021-05)
    if 'TXFREQUENCY' in df_link.columns and 'RXFREQUENCY' in df_link.columns:
        name_cols.update({
            'RXFREQUENCY' : 'RxFrequency',
            'TXFREQUENCY' : 'TxFrequency'
        })
    
    df_link = df_link.rename(columns = name_cols).reindex(columns = list(name_cols.values()))
    _log('Renamed and reordered link dataframe columns')
    
    
    _log('\n******************************** MERGE INCA ********************************')
    
    # set INCA RR data for each km of link path
    df_link['RRpath'] = get_inca_data(
        df_link['DateTime'],
        df_link['XStart'],
        df_link['YStart'],
        df_link['XEnd'],
        df_link['YEnd'],
        df_link['PathLength']
     )
    _log('Merged INCA RR data to transmissions')
    
    
    _log('\n******************************** SAVE FILES ********************************')
    
    # build path for clean config and transmissions destination files
    dest_config = file_config.parent.parent / 'clean' / f'{file_config.stem}{file_config.suffix}'
    dest_trans = file_trans.parent.parent / 'clean' / f'{file_trans.stem}{file_trans.suffix}'
    
    # build path for clean link destination file (same folder, date and extension as transmissions file)
    date = str(file_trans.stem).split('_')[-1]
    dest_link = pathlib.Path(dest_trans.parents[0], f'LINK_{date}{file_trans.suffix}')
    
    
    # save clean files
    df_config.to_csv(dest_config, sep = ';', header = True, index = False)
    _log(f'Saved clean config file with shape {df_config.shape} to "{str(dest_config)}"')
    df_trans.to_csv(dest_trans, sep = ';', header = True, index = False)
    _log(f'Saved clean transmissions file with shape {df_trans.shape} to "{str(dest_trans)}"')
    df_link.to_csv(dest_link, sep = ';', header = True, index = False)
    _log(f'Saved clean link file with shape {df_link.shape} to "{str(dest_link)}"')


if __name__ == '__main__':
    # flag for data prep
    start_prep = True
    
    # get config and transmissions file from arguments
    args = parse_arguments()
    
    # convert config and transmissions arguments to paths
    file_config = pathlib.Path(args.config)
    file_trans = pathlib.Path(args.transmissions)
    dir_inca = pathlib.Path(args.inca)
    
    # check if config files exists
    if not file_config.exists():
        _log('Invalid path for config file!')
        start_prep = False
        
    # check if transmissions file exists
    if not file_trans.exists():
        _log('Invalid path for transmissions file!')
        start_prep = False
    
    # chec if inca dir exists
    if not dir_inca.exists():
        _log('Invalid path for inca directory!')
        start_prep = False
    
    # load areas json file
    areas_poly = load_areas(pathlib.Path('areas.json'))
    
    # start prep if flag is True, otherwise exit with code 2
    if start_prep:
        prep()
    else:
        sys.exit(2)