[i] calc. INCA RR data for each km of path length

FHSTP/prep.py

@@ -2,8 +2,8 @@
 # Title: Data Preparation for LINK Configs and Transmissions
 # Author: Fabian Kovac <ds191008@fhstp.ac.at>
 # Team: University of Applied Sciences St. Pölten
-# Version: 1.3
-# Last changed: 2021-06-23
+# Version: 1.4
+# Last changed: 2021-06-24
 #
 
 import sys
@@ -31,7 +31,7 @@ def parse_arguments() -> argparse.Namespace:
     
     parser = argparse.ArgumentParser(
         prog = 'prep.py',
-        usage = 'python %(prog)s -c <config_file> -t <transmissions_file>',
+        usage = 'python %(prog)s -c <config_file> -t <transmissions_file> -i <inca_dir>',
         description = desc,
         formatter_class = argparse.RawTextHelpFormatter
     )
@@ -82,7 +82,7 @@ def load_inca_file(file_inca: pathlib.Path) -> np.array:
     # 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, sep = ' ')
+        x = np.fromstring(x, dtype = np.float32, sep = ' ')
         x = np.reshape(x, (401, 701))
 
         return x
@@ -246,25 +246,24 @@ def lambert_to_inca_idx(x: np.array, y: np.array) -> tuple:
     return np.round(x/1000, decimals = 0).astype(int), np.round(y/1000, decimals = 0).astype(int)
 
 
-def get_inca_data(inca_data: np.array, datetimes: np.array, x: np.array, y: np.array) -> np.array:
-    """Get INCA RR data based on Lamber Conic Conformal Coordinates
+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:
-        inca_data (np.array): Tensor containing INCA RR data
         datetimes (np.array): Vector containing datetimes of transmissions
-        x (np.array): Vector containing x values of LINK
-        y (np.array): Vector containing y values of LINK
+        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 based on datetimes and lambert coordinates
+        inca_RR (np.array): Vector containing INCA RR data for each transmission
 	"""
     
-    # convert utm coordinates to lambert conic conformal projection
-    lccX, lccY = utm_to_lambert(x, y)
-    
-    # convert lambert coordinates to INCA indices
-    idx, idy = lambert_to_inca_idx(lccX, lccY)
-    
+    # 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
@@ -275,16 +274,38 @@ def get_inca_data(inca_data: np.array, datetimes: np.array, x: np.array, y: np.a
     
     # 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 utm coordinates to lambert conic conformal projection
+    xa, ya = utm_to_lambert(lon_a, lat_a)
+    xb, yb = utm_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
+    n = np.round(lengths, decimals = 0).astype(int)
+    n = np.where(n < 3, 3, n)
     
-    # get INCA RR data based on time indices and lambert coordinates
-    inca_RR = inca_data[idx_times, idy, idx]
+    # set list for INCA RR data
+    inca_RR = []
     
-    # TODO: discuss retrieving INCA RR data
-    # currently data is based on the exact point of the link
-    # maybe get mean RR of 5x5 km grid around point?
-    # maybe get RR data based on conv2d with a 5x5 gaussian filter around point?
-    # --> let's discuss!
+    # 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
 
@@ -476,10 +497,10 @@ def prep() -> None:
         'DIFFLEVEL' : 'DiffLevel',
         'LONGITUDE_A' : 'XStart',
         'LATITUDE_A' : 'YStart',
-        'LONGITUDE_B' : 'XEnd',
-        'LATITUDE_B' : 'YEnd',
         'LONGITUDE_MID' : 'XMid',
         'LATITUDE_MID' : 'YMid',
+        'LONGITUDE_B' : 'XEnd',
+        'LATITUDE_B' : 'YEnd',
         'LENGTH' : 'PathLength',
         'FREQUENCY' : 'Frequency',
     }
@@ -497,15 +518,9 @@ def prep() -> None:
     
     _log('\n******************************** MERGE INCA ********************************')
     
-    # 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}')
-    
-    # set INCA RR data based on datetime and coordinates
-    df_link['RRStart'] = get_inca_data(inca_data, df_link['DateTime'], df_link['XStart'], df_link['YStart'])
-    df_link['RREnd'] = get_inca_data(inca_data, df_link['DateTime'], df_link['XEnd'], df_link['YEnd'])
-    df_link['RRMid'] = get_inca_data(inca_data, df_link['DateTime'], df_link['XMid'], df_link['YMid'])
-    _log('Merged INCA RR data to LINK dataframe')
+    # 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 ********************************')