Develop

.travis.yml

@@ -14,9 +14,10 @@ install:
   - conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION numpy matplotlib 
   - source activate test-environment
   - pip install pytest pytest-cov coveralls rasterio
+  - pip install pytest pytest-cov coveralls opencv-python-headless
 
 script:
-  pytest --cov-report term-missing:skip-covered --cov=coral tests/
+  pytest --cov-report term-missing:skip-covered --cov=coral tests/test_coral.py
 
 after_success:
   - coveralls

Calc_CR_response.py

@@ -1,55 +1,106 @@
 """
-Ths is an example python script for running CoRAL analysis
+Ths is the main python script for running CoRAL analysis
 
-Example command line usage: python3 Calc_CR_response.py
+Example command line usage: python3 Calc_CR_response_config.py /path/to/config_file
+
+An example for the config_file is here: ./data/coral_serf.conf
+An example for the input file list is here: ./data/mli_desc.list
+An example for the file containing target coordinates is here: ./data/site_lat_lon_hgt_date1_date2_az_rg_initial.txt
+Test run: python3 Calc_CR_response.py ./data/coral_serf.conf
+
+required packages to be installed (on Gadi):
+# pip install --user joblib
+# pip install --user rasterio
 """
-import glob
-import numpy as np
 from datetime import datetime
+from joblib import Parallel, delayed
+import multiprocessing
+import sys, os.path
+from coral import config as cf
+from coral.dataio import read_input_files, write_radar_coords, plot_results
 from coral.corner_reflector import loop
-from coral.plot import *
-
-# start/end date for plots
-start = datetime(2018, 7, 1)
-end = datetime(2018, 11, 1)
-
-# Get list of image files in current directory
-files = []
-for file in glob.glob("data/*.tif"):
-    files.append(file)
-
-files.sort()
-
-# final target positions manual entry
-sites = {
-#   site name            rg   az
-    'SERF' : np.array([[ 87, 110]]),
-}
-
-#cropped image size
-sub_im = 51
-
-# define target_window size
-targ_win_sz = 5
-
-# define clutter window size
-clt_win_sz = 9
-
-# loop through all corner reflector sites
-for name, cr in sites.items():
-    print(name,'is desc',cr[0])
-
-    avgI, rcs, scr, Avg_clt, t = loop(files, sub_im, cr[0], targ_win_sz, clt_win_sz)
-
-    cr_pos = np.array([sub_im, sub_im])
-
-    # Plot mean intensity image
-    plot_mean_intensity(avgI, cr_pos, targ_win_sz, clt_win_sz, name)
-
-    # Plot RCS_SCR time series
-    plot_rcs_scr(t, rcs, scr, start, end, name)
-
-    # Plot average clutter time series
-    plot_clutter(t, Avg_clt, start, end, name)
 
 
+# check if config-file has been given as an argument of the main script call
+print('')
+if len(sys.argv) != 2:
+    print('Exiting: Provide path to config-file as command line argument')
+    print('')
+    print('Usage: python3 Calc_CR_response.py <config-file>')
+    exit()
+else:
+    config_file = sys.argv[1]
+    print(f'Looking for CoRAL input data defined in {config_file}')
+
+# read config-file and convert parameters to required data type
+params = cf.get_config_params(config_file)
+
+# General screen output
+print(f'Results will be saved into {params[cf.OUT_DIR]}')
+# Start the processing
+print(' ')
+print('Running CoRAL...')
+print(' ')
+# used to calculate runtime
+run_start = datetime.now()
+
+# check and read paths to input data
+files_a, files_d, sites = read_input_files(params)
+
+
+#######################################
+# function for parallel loop processing
+sitenames = sites.keys()
+# note that dictionaries are unsorted and the variable name is hence not ordered
+def processInput(name):
+
+    cr = sites[name]
+
+    if files_a:
+        avgI_a, rcs_a, scr_a, clt_a, t_a, cr_new_a, cr_pos_a = loop(files_a, params[cf.SUB_IM], cr[0], params[cf.TARG_WIN_SZ], params[cf.CLT_WIN_SZ])
+    if files_d:
+        avgI_d, rcs_d, scr_d, clt_d, t_d, cr_new_d, cr_pos_d = loop(files_d, params[cf.SUB_IM], cr[1], params[cf.TARG_WIN_SZ], params[cf.CLT_WIN_SZ])
+
+    if files_a and files_d:
+        return name, cr_pos_a, avgI_a, rcs_a, scr_a, clt_a, t_a, cr_new_a, \
+               cr_pos_d, avgI_d, rcs_d, scr_d, clt_d, t_d, cr_new_d
+    elif files_a and not files_d:
+        return name, cr_pos_a, avgI_a, rcs_a, scr_a, clt_a, t_a, cr_new_a
+    elif files_d and not files_a:
+        return name, cr_pos_d, avgI_d, rcs_d, scr_d, clt_d, t_d, cr_new_d
+
+
+# parallel processing of MLI read and RCS, SCR calculation
+num_cores = multiprocessing.cpu_count()
+# num_cores results in 32 on the NCI which in turn results in an error
+# hence the number of 16 cores is hard-coded here
+results = Parallel(n_jobs=params[cf.N_JOBS])(delayed(processInput)(name) for name in sitenames)
+#######################################
+
+
+# extract results and plot images
+print(' ') 
+print('Creating output data...')
+print(' ') 
+
+# create output dir if it doesn't exist
+if not os.path.exists(params[cf.OUT_DIR]):
+    os.makedirs(params[cf.OUT_DIR])
+
+plot_results(sites, results, params)
+
+
+# write updated radar coordinates to a new file
+print(' ')
+if files_a:
+    # new file with updated radar coordinates of CRs
+    write_radar_coords(params[cf.ASC_CR_FILE_ORIG], params[cf.ASC_CR_FILE_NEW], sites, "asc")
+if files_d:
+    # new file with updated radar coordinates of CRs
+    write_radar_coords(params[cf.DESC_CR_FILE_ORIG], params[cf.DESC_CR_FILE_NEW], sites, "desc")
+
+
+# print out runtime
+runtime = datetime.now() - run_start
+print(' ')
+print('Runtime: %s' % runtime)

README.md

@@ -43,7 +43,7 @@ A Jupyter notebook is included in the package that demonstrates an example CoRAL
 
 To run unit tests from within CoRAL repository directory: 
 
-```python -m unittest discover tests/```
+```python3 -m unittest discover tests/```
 
 ### License
 

coral/config.py

@@ -0,0 +1,134 @@
+"""
+This Python module contains utilities to parse CoRAL configuration files.
+Most of this code is copied from PyRate (see https://github.com/GeoscienceAustralia/PyRate)
+"""
+from typing import Dict
+import os
+
+
+# constants for lookups
+#: STR; Name of input interferogram list file
+ASC_LIST = 'backscatter_data_asc'
+DESC_LIST = 'backscatter_data_desc'
+ASC_CR_FILE_ORIG = 'cr_file_asc'
+ASC_CR_FILE_NEW = 'cr_file_asc_new'
+DESC_CR_FILE_ORIG = 'cr_file_desc'
+DESC_CR_FILE_NEW = 'cr_file_desc_new'
+OUT_DIR = 'path_out'
+TARG_WIN_SZ = 'target_window_size'
+CLT_WIN_SZ = 'clutter_window_size'
+SUB_IM = 'sub_image_size'
+N_JOBS = 'n_jobs'
+YMAX_RCS = 'ymax_rcs'
+YMAX_SCR = 'ymax_scr'
+YMIN_CLUTTER = 'ymin_clutter'
+YMAX_CLUTTER = 'ymax_clutter'
+
+
+# Lookup to help convert args to correct type/defaults
+# format is	key : (conversion, default value)
+PARAM_CONVERSION = {
+    TARG_WIN_SZ: (int, 3),
+    CLT_WIN_SZ: (int, 7),
+    SUB_IM: (int, 51),
+    N_JOBS: (int, 16),
+    YMAX_RCS: (float, 35),
+    YMAX_SCR: (float, 30),
+    YMIN_CLUTTER: (float, -16),
+    YMAX_CLUTTER: (float, -2),
+}
+
+# path variables
+PATHS = [
+    ASC_LIST,
+    DESC_LIST,
+    ASC_CR_FILE_ORIG,
+    ASC_CR_FILE_NEW,
+    DESC_CR_FILE_ORIG,
+    DESC_CR_FILE_NEW,
+    OUT_DIR,
+]
+
+
+def get_config_params(path: str) -> Dict:
+    """
+    Reads the parameters file provided by the user and converts it into a dictionary.
+
+    :param str path: path to config file
+    :return: dict params: config parameters
+    """
+    txt = ''
+    with open(path, 'r') as inputFile:
+        for line in inputFile:
+            if any(x in line for x in PATHS):
+                pos = line.find('~')
+                if pos != -1:
+                    # create expanded line
+                    line = line[:pos] + os.environ['HOME'] + line[(pos + 1):]
+            txt += line
+    params = _parse_conf_file(txt)
+
+    return params
+
+
+def _parse_conf_file(content) -> Dict:
+    """
+    Converts the parameters from their text form into a dictionary.
+
+    :param str content: Parameters as text
+    :return: dict params: config parameters
+    """
+
+    def _is_valid(line):
+        """
+        Check if line is not empty or has % or #
+        """
+        return line != "" and line[0] not in "%#"
+
+    lines = [ln.split() for ln in content.split('\n') if _is_valid(ln)]
+
+    # convert "field:   value" lines to [field, value]
+    kvpair = [(e[0].rstrip(":"), e[1]) for e in lines if len(e) == 2] + \
+             [(e[0].rstrip(":"), None) for e in lines if len(e) == 1]
+    parameters = dict(kvpair)
+    for p in PATHS:
+        if p not in parameters:
+            parameters[p] = None
+
+    if not parameters:
+        raise ConfigException('Cannot parse any parameters from config file')
+
+    return _parse_pars(parameters)
+
+
+# todo: check why conversion of parameters to int is not working properly
+def _parse_pars(pars) -> Dict:
+    """
+    Takes dictionary of parameters, converting values to required type
+    and providing defaults for missing values.
+
+    :param dict pars: config parameters (as strings)
+    :return: dict params: converted config parameters (according to PARAM_CONVERSION lookup table)
+    """
+    # Fallback to default for missing values and perform conversion.
+    for k in PARAM_CONVERSION:
+        if pars.get(k) is None:
+            pars[k] = PARAM_CONVERSION[k][1]
+            # _logger.warning(f"No value found for parameter '{k}'. Using "f"default value {pars[k]}.")
+        else:
+            conversion_func = PARAM_CONVERSION[k][0]
+            if conversion_func:
+                try:
+                    pars[k] = conversion_func(pars[k])
+                except ValueError as e:
+                    _logger.error(
+                        f"Unable to convert '{k}': {pars[k]} to " f"expected type {conversion_func.__name__}.")
+                    raise e
+
+    return pars
+
+
+class ConfigException(Exception):
+    """
+    Default exception class for configuration errors.
+    """

coral/corner_reflector.py

@@ -30,19 +30,38 @@ def loop(files, sub_im, cr, targ_win_sz, clt_win_sz):
 
     cr_pos = np.array([sub_im, sub_im])
 
+
+    # find location of pixel with highest intensity inside target window
+    # calculate target window bounds
+    xmin_t = int(np.ceil(cr_pos[0] - targ_win_sz / 2))
+    xmax_t = int(np.floor(cr_pos[0] + targ_win_sz / 2 + 1))
+    ymin_t = int(np.ceil(cr_pos[1] - targ_win_sz / 2))
+    ymax_t = int(np.floor(cr_pos[1] + targ_win_sz / 2 + 1))
+    # crop mean intensity matrix to target window size
+    avgI_t = avgI[ymin_t:ymax_t, xmin_t:xmax_t]
+    # find matrix position of maximum mean intesity
+    max_ix = np.unravel_index(np.argmax(avgI_t, axis=None), avgI_t.shape)
+    # calculate shift w.r.t. central pixel
+    centre_ix =  (targ_win_sz - 1) / 2
+    shift = np.array([max_ix[1] - centre_ix, max_ix[0] - centre_ix], dtype=int)
+    # updated cr position
+    cr_new = cr + shift
+    # also update cr_pos accordingly
+    cr_pos = cr_pos + shift
+
+
     # calculate target energy
     En, Ncr = calc_integrated_energy(d, cr_pos, targ_win_sz)
     # calculate clutter energy for window centred in same spot as target window
     E1, N1 = calc_integrated_energy(d, cr_pos, clt_win_sz)
     # calculate average clutter
     Avg_clt, Eclt, Nclt = calc_clutter_intensity(En, E1, Ncr, N1)
-    #
+    # Calculate total energy, signal-to-clutter ratio and radar cross section
     Ecr = calc_total_energy(Ncr, Nclt, Eclt, En)
     scr = calc_scr(Ecr, Eclt, Nclt)
-
     rcs = calc_rcs(Ecr, rho_r, rho_a, theta)
 
-    return avgI, rcs, scr, Avg_clt, t
+    return avgI, rcs, scr, Avg_clt, t, cr_new, cr_pos
 
 
 def get_win_bounds(pos, winsz):