diff --git a/.travis.yml b/.travis.yml index 1493ef2..cd417b3 100644 --- a/.travis.yml +++ b/.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 diff --git a/Calc_CR_response.py b/Calc_CR_response.py index 196a76c..c659f9d 100644 --- a/Calc_CR_response.py +++ b/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 ') + 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) diff --git a/README.md b/README.md index f5511d2..f011592 100644 --- a/README.md +++ b/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 diff --git a/coral/config.py b/coral/config.py new file mode 100644 index 0000000..370ff94 --- /dev/null +++ b/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. + """ diff --git a/coral/corner_reflector.py b/coral/corner_reflector.py index 412f71b..d589a78 100644 --- a/coral/corner_reflector.py +++ b/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): diff --git a/coral/dataio.py b/coral/dataio.py index 0611255..0801808 100644 --- a/coral/dataio.py +++ b/coral/dataio.py @@ -3,9 +3,89 @@ Radar images """ import numpy as np -import rasterio +import rasterio as rio from rasterio.windows import Window -import os +import sys, os +from coral.plot import * +from coral.plot2 import * + + +def read_input_files(params): + """ + Reads input files and performs checks on existence + + :param dict params: config parameters + :return: list files_a: containing paths to asc backscatter files (none if not given) + :return: list files_d: containing paths to desc backscatter files (none if not given) + :return: dict sites: containing radar target sites and array with range and azimuth coordinates + (-999 is used for non-existent geometry if applicable) + """ + # if no datapath is given in the config file, this geometry is not used + # read asc pass data files + if params[cf.ASC_LIST]: + if not os.path.exists(params[cf.ASC_LIST]): + raise Exception(f'{params[cf.ASC_LIST]} does not exist') + else: + print(f'Reading data from file list {params[cf.ASC_LIST]}') + # read the paths to ascending-pass backscatter images + with open(params[cf.ASC_LIST]) as f_in: + files_a = [line.rstrip() for line in f_in] + files_a.sort() + # read the ascending-pass radar coordinates of targets + filename = params[cf.ASC_CR_FILE_ORIG] + if not os.path.exists(params[cf.ASC_CR_FILE_ORIG]): + raise Exception(f'{params[cf.ASC_CR_FILE_ORIG]} does not exist') + else: + print('') + sites_a, az_a, rg_a = read_radar_coords(filename) + # set files variable to None if no asc file list supplied + else: + print('No ascending-pass file list supplied') + files_a = None + # read desc pass data files + if params[cf.DESC_LIST]: + if not os.path.exists(params[cf.DESC_LIST]): + raise Exception(f'{params[cf.DESC_LIST]} does not exist') + else: + print(f'Reading data from file list {params[cf.DESC_LIST]}') + # read the paths to descending-pass backscatter images + with open(params[cf.DESC_LIST]) as f_in: + files_d = [line.rstrip() for line in f_in] + files_d.sort() + # read the descending-pass radar coordinates of targets + filename = params[cf.DESC_CR_FILE_ORIG] + if not os.path.exists(params[cf.DESC_CR_FILE_ORIG]): + raise Exception(f'{params[cf.DESC_CR_FILE_ORIG]} does not exist') + else: + print('') + sites_d, az_d, rg_d = read_radar_coords(filename) + # set files variable to None if no desc file list supplied + else: + print('No descending-pass file list supplied') + files_d = None + + # check if CR files have the same length + if files_a and files_d: + if sites_a == sites_d: + print("CR files for asc and desc tracks contain same set of sites.") + keys = sites_a + values = (np.array([[rg_a[i], az_a[i]], [rg_d[i], az_d[i]]], dtype=int) \ + for i in range(0, len(sites_a))) + else: + print("ERROR: different CR sites given for asc and desc track.") + sys.exit() + print(' ') + elif files_a and not files_d: + keys = sites_a + values = (np.array([[rg_a[i], az_a[i]], [-999, -999]], dtype=int) \ + for i in range(0, len(sites_a))) + elif files_d and not files_a: + keys = sites_d + values = (np.array([[-999, -999], [rg_d[i], az_d[i]]], dtype=int) \ + for i in range(0, len(sites_d))) + sites = dict(zip(keys, values)) + + return files_a, files_d, sites def readfile(file, sub_im, cr): @@ -16,7 +96,12 @@ def readfile(file, sub_im, cr): if ext == '.tif': print('Reading tiff image:', file) par = readpar(root + '.mli.par') - data = readtiff(file, sub_im, cr) + data = readimg(file, sub_im, cr) + + elif ext == '.img': + print('Reading ENVI image:', file) + #par = readhdr(root + '.hdr') + data = readimg(file, sub_im, cr) else: # must be GAMMA flat binary float format print('Reading flat binary image', file) @@ -60,12 +145,190 @@ def readmli(datafile, par, sub_im, cr): return d[cr[1]-sub_im:cr[1]+sub_im,cr[0]-sub_im:cr[0]+sub_im] -def readtiff(datafile, sub_im, cr): - """Function to read a tiff and provide a subsetted image""" +def readimg(datafile, sub_im, cr): + """Function to read a file and provide a subsetted image""" - with rasterio.open(datafile) as src: + with rio.open(datafile) as src: d = src.read(1, window=Window(cr[0]-sub_im, cr[1]-sub_im, sub_im*2, sub_im*2)) #print("Number of elements and size of the array is",d.size, d.shape) #d[d==0]= np.nan # convert zeros to nan return d + + +def read_radar_coords(filename): + print("Reading textfile with CR positions...") + site = [] + az = [] + rg = [] + if os.path.isfile(filename) and os.access(filename, os.R_OK): + print("File", filename, "exists and is readable.") + f = open(filename) + lines = f.readlines() + idx = 0 + for line in lines: + # get site name + line = line.strip("\n") + site.append(line.split("\t")[0]) + az.append(line.split("\t")[-2]) + rg.append(line.split("\t")[-1]) + idx = idx + 1 + print("Radar coordinates at %d sites read" % (idx)) + f.close() + else: + print("ERROR: Can't read file", filename) + sys.exit() + print() + return site, az, rg + + +def plot_results(sites, results, params): + """ + Function to plot results of radar target response analysis, also outputs RCS values to file + + :param dict sites: containts site name and radar coordinates + :param: dict results: containing paths to asc backscatter files (none if not given) + :param dict params: config parameters + """ + sitenames = sites.keys() + for i in range(0, len(sitenames)): + + # ascending and descending CRs in one plot + if params[cf.ASC_LIST] and params[cf.DESC_LIST]: + # read result arrays of parallel function, both geometries + name = results[i][0] + cr_pos_a = results[i][1] + avgI_a = results[i][2] + rcs_a = results[i][3] + scr_a = results[i][4] + clt_a = results[i][5] + t_a = results[i][6] + cr_new_a = results[i][7] + cr_pos_d = results[i][8] + avgI_d = results[i][9] + rcs_d = results[i][10] + scr_d = results[i][11] + clt_d = results[i][12] + t_d = results[i][13] + cr_new_d = results[i][14] + + # add new coords to sites dictionary + cr = np.array([cr_new_a, cr_new_d]) + sites[name] = cr + + # Visualisation + print('Site %s' % name) + # Plot mean intensity image + plot_mean_intensity2(avgI_a, avgI_d, cr_pos_a, cr_pos_d, name, params) + + # extract start and end date + start_time = min(t_a[0], t_d[0]) + end_time = max(t_a[-1], t_d[-1]) + margin = (end_time - start_time) / 50 + start = start_time - margin + end = end_time + margin + + # Plot average clutter time series + plot_clutter2(t_a, t_d, clt_a, clt_d, start, end, name, params) + # Plot scr time series + plot_scr2(t_a, t_d, scr_a, scr_d, start, end, name, params) + # Plot rcs time series + plot_rcs2(t_a, t_d, rcs_a, rcs_d, start, end, name, params) + + # write RCS data to file + filename_out = params[cf.OUT_DIR] + "/rcs_values_" + name + "_" + "Ascending.txt" + fout = open(filename_out, 'w') + for time, value in zip(t_a, rcs_a): + timestr = time.strftime("%Y%m%d") + fout.write("%s %f\n" % (timestr, value)) + fout.close() + filename_out = params[cf.OUT_DIR] + "/rcs_values_" + name + "_" + "Descending.txt" + fout = open(filename_out, 'w') + for time, value in zip(t_d, rcs_d): + timestr = time.strftime("%Y%m%d") + fout.write("%s %f\n" % (timestr, value)) + fout.close() + + # one geometry (ascending or descending) + else: + # read result arrays of parallel function, one geometry + name = results[i][0] + cr_pos = results[i][1] + avgI = results[i][2] + rcs = results[i][3] + scr = results[i][4] + clt = results[i][5] + t = results[i][6] + cr_new = results[i][7] + + # add new coords to sites dictionary + if params[cf.ASC_LIST]: + geom = 'Ascending' + cr = np.array([cr_new, [-999, -999]]) + if params[cf.DESC_LIST]: + geom = 'Descending' + cr = np.array([[-999, -999], cr_new]) + sites[name] = cr + + # Visualisation + print('Site %s' % name) + # Plot mean intensity image + plot_mean_intensity(avgI, cr_pos, name, params) + + # extract start and end date + start_time = t[0] + end_time = t[-1] + margin = (end_time - start_time) / 50 + start = start_time - margin + end = end_time + margin + + # Plot average clutter time series + plot_clutter(t, clt, start, end, name, geom, params) + # Plot scr time series + plot_scr(t, scr, start, end, name, geom, params) + # Plot rcs time series + plot_rcs(t, rcs, start, end, name, geom, params) + + # write RCS data to file + filename_out = params[cf.OUT_DIR] + "/rcs_values_" + name + "_" + geom + ".txt" + fout = open(filename_out, 'w') + for time, value in zip(t, rcs): + timestr = time.strftime("%Y%m%d") + fout.write("%s %f\n" % (timestr, value)) + fout.close() + + +def write_radar_coords(filename_init, filename, sites, geom): + print("Writing new CR positions to textfile...") + fout = open(filename,'w') + if os.path.isfile(filename_init) and os.access(filename_init, os.R_OK): + print("File", filename_init, "exists and is readable.") + f = open(filename_init) + lines = f.readlines() + for line in lines: + # get site name + name = line.split("\t")[0] + cr = sites[name] + temp = line.split("\t")[:-2] + temp2 = "\t".join(temp) + if geom == "asc": + out_line = temp2 + "\t" + str(cr[0][1]) + \ + "\t" + str(cr[0][0]) + "\n" + if geom == "desc": + out_line = temp2 + "\t" + str(cr[1][1]) + \ + "\t" + str(cr[1][0]) + "\n" + fout.write(out_line) + f.close() + fout.close() + else: + print("ERROR: Can't read file", filename_init) + sys.exit() + print() + print("Azimuth and Range number of CR pixels has been written to " + \ + filename + ".") + print() + return + + +class Exception(Exception): + """IO generic exception class""" \ No newline at end of file diff --git a/coral/plot.py b/coral/plot.py index e040ac8..1954104 100644 --- a/coral/plot.py +++ b/coral/plot.py @@ -2,14 +2,14 @@ This python module contains functions for plotting CoRAL output """ import matplotlib -matplotlib.use('Agg') +matplotlib.use('agg') from matplotlib import pyplot as plt from matplotlib.patches import RegularPolygon -from matplotlib.colors import LogNorm import numpy as np +from coral import config as cf -def plot_mean_intensity(avgI, cr_pos, targ_win_sz, clt_win_sz, name): +def plot_mean_intensity(avgI, cr_pos, name, params): '''Plot image of mean SAR intensity''' # set black/white colormap for plots cmap = plt.set_cmap('gist_gray') @@ -22,11 +22,11 @@ def plot_mean_intensity(avgI, cr_pos, targ_win_sz, clt_win_sz, name): #cax = ax2.matshow(avgI_d, vmin=-20, vmax=10, cmap=cmap) # define target window - p1 = RegularPolygon(cr_pos, 4, (targ_win_sz/2)+1, \ + p1 = RegularPolygon(cr_pos, 4, (params[cf.TARG_WIN_SZ]/2)+1, \ orientation = np.pi / 4, linewidth=1, \ edgecolor='r',facecolor='none') # define clutter window - p2 = RegularPolygon(cr_pos, 4, (clt_win_sz/2)+2, \ + p2 = RegularPolygon(cr_pos, 4, (params[cf.CLT_WIN_SZ]/2)+2, \ orientation = np.pi / 4, linewidth=1, \ edgecolor='y',facecolor='none') @@ -56,44 +56,90 @@ def plot_mean_intensity(avgI, cr_pos, targ_win_sz, clt_win_sz, name): #fig.set_size_inches(w=6,h=4) # save PNG file - fig.savefig('mean_intensity_%s.png' % name, dpi=300, bbox_inches='tight') - + filename = params[cf.OUT_DIR] + "/mean_intensity_" + name + ".png" + fig.savefig(filename, dpi=300, bbox_inches='tight') + + # avoid "RuntimeWarning: More than 20 figures have been opened" + # by closing all open figures + plt.close('all') + return -def plot_rcs_scr(t, rcs, scr, start, end, name): - '''Plot RCS and SCR time series''' +def plot_clutter(t, clt, start, end, name, geom, params): + '''Plot average clutter time series''' fig = plt.figure() ax = fig.add_subplot(1,1,1) - plt.plot(t, rcs, 'ro-', label='RCS') - plt.plot(t, scr, 'bo-', label='SCR') + plt.plot(t, clt, 'bo-', label=geom) plt.xlim(start, end) - plt.ylim(0, 40) + plt.ylim(params[cf.YMIN_CLUTTER], params[cf.YMAX_CLUTTER]) plt.xlabel('Date') - plt.ylabel('RCS / SCR (dB)') - plt.legend(loc=4) + plt.ylabel('Average Clutter (dB)') + plt.legend(loc=1) plt.grid(True) - plt.title('Corner Reflector response at %s' % name) + plt.title('Average Clutter at %s' % name) for label in ax.get_xticklabels(): label.set_rotation(90) - fig.savefig('rcs_scr_%s.png' % name, dpi=300, bbox_inches='tight') + + # save PNG file + filename = params[cf.OUT_DIR] + "/clutter_" + name + ".png" + fig.savefig(filename, dpi=300, bbox_inches='tight') + + # avoid "RuntimeWarning: More than 20 figures have been opened" + # by closing all open figures + plt.close('all') return -def plot_clutter(t, clt, start, end, name): - '''Plot average clutter time series''' + +def plot_scr(t, scr, start, end, name, geom, params): + '''Plot RCS time series''' fig = plt.figure() ax = fig.add_subplot(1,1,1) - plt.plot(t, clt, 'bo-', label='Clutter') + plt.plot(t, scr, 'bo-', label=geom) plt.xlim(start, end) - plt.ylim(-16, -2) + plt.ylim(0, params[cf.YMAX_SCR]) plt.xlabel('Date') - plt.ylabel('Average Clutter (dB)') - plt.legend(loc=1) + plt.ylabel('SCR (dB)') + plt.legend(loc=4) plt.grid(True) - plt.title('Average Clutter at %s' % name) + plt.title('Signal to Clutter Ratio at site %s' % name) for label in ax.get_xticklabels(): label.set_rotation(90) - fig.savefig('clutter_%s.png' % name, dpi=300, bbox_inches='tight') + + # save PNG file + filename = params[cf.OUT_DIR] + "/scr_" + name + ".png" + fig.savefig(filename, dpi=300, bbox_inches='tight') + + # avoid "RuntimeWarning: More than 20 figures have been opened" + # by closing all open figures + plt.close('all') + + return + +def plot_rcs(t, rcs, start, end, name, geom, params): + '''Plot RCS time series''' + fig = plt.figure() + ax = fig.add_subplot(1,1,1) + plt.plot(t, rcs, 'bo-', label=geom) + plt.xlim(start, end) + plt.ylim(0, params[cf.YMAX_RCS]) + plt.xlabel('Date') + plt.ylabel('RCS (dB$\mathregular{m^2}$)') + plt.legend(loc=4) + plt.grid(True) + plt.title('Radar Cross Section at site %s' % name) + for label in ax.get_xticklabels(): + label.set_rotation(90) + + # save PNG file + filename = params[cf.OUT_DIR] + "/rcs_" + name + ".png" + fig.savefig(filename, dpi=300, bbox_inches='tight') + + # avoid "RuntimeWarning: More than 20 figures have been opened" + # by closing all open figures + plt.close('all') + return + diff --git a/coral/plot2.py b/coral/plot2.py new file mode 100644 index 0000000..77ff65a --- /dev/null +++ b/coral/plot2.py @@ -0,0 +1,153 @@ +""" +This python module contains functions for plotting CoRAL output +""" +import matplotlib +matplotlib.use('agg') +from matplotlib import pyplot as plt +from matplotlib.patches import RegularPolygon +import numpy as np +from coral import config as cf + + +def plot_mean_intensity2(avgI_a, avgI_d, cr_pos_a, cr_pos_d, name, params): + '''Plot image of mean SAR intensity''' + # set black/white colormap for plots + cmap = plt.set_cmap('gist_gray') + + # draw new plot + fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(14,5)) + axlist = [ax1, ax2] + im1 = ax1.matshow(avgI_d, vmin=-20, vmax=10, cmap=cmap) + im2 = ax2.matshow(avgI_a, vmin=-20, vmax=10, cmap=cmap) + + # define target window + p1 = RegularPolygon(cr_pos_d, 4, (params[cf.TARG_WIN_SZ]/2)+2, \ + orientation = np.pi / 4, linewidth=1, \ + edgecolor='r',facecolor='none') + # define clutter window + p2 = RegularPolygon(cr_pos_d, 4, (params[cf.CLT_WIN_SZ]/2)+3, \ + orientation = np.pi / 4, linewidth=1, \ + edgecolor='y',facecolor='none') + # define target window + p3 = RegularPolygon(cr_pos_a, 4, (params[cf.TARG_WIN_SZ]/2)+2, \ + orientation = np.pi / 4, linewidth=1, \ + edgecolor='r',facecolor='none') + # define clutter window + p4 = RegularPolygon(cr_pos_a, 4, (params[cf.CLT_WIN_SZ]/2)+3, \ + orientation = np.pi / 4, linewidth=1, \ + edgecolor='y',facecolor='none') + # add windows to plot + ax1.add_patch(p1) + ax1.add_patch(p2) + ax2.add_patch(p3) + ax2.add_patch(p4) + + # add text labels + ax1.text(45, 42, name, color='w', fontsize=10) + ax2.text(45, 42, name, color='w', fontsize=10) + # plot labels + ax1.set_xlabel('Range') + ax2.set_xlabel('Range') + ax1.set_ylabel('Azimuth') + # add colorbar + cbar = fig.colorbar(im1, ax=axlist) + cbar.set_label('dB') + # add title + #fig.set_title('Mean intensity at site %s' % name) + ax1.set_title('Descending') + ax2.set_title('Ascending') + # x-axis labels at bottom + ax1.xaxis.set_tick_params(labeltop='off', labelbottom='on') + ax2.xaxis.set_tick_params(labeltop='off', labelbottom='on') + + # save PNG file + filename = params[cf.OUT_DIR] + "/mean_intensity_" + name + ".png" + fig.savefig(filename, dpi=300, bbox_inches='tight') + + # avoid "RuntimeWarning: More than 20 figures have been opened" + # by closing all open figures + plt.close('all') + + return + + +def plot_clutter2(t_a, t_d, clt_a, clt_d, start, end, name, params): + '''Plot average clutter time series''' + fig = plt.figure() + ax = fig.add_subplot(1,1,1) + plt.plot(t_a, clt_a, 'ro-', label='Asc') + plt.plot(t_d, clt_d, 'bo-', label='Desc') + plt.xlim(start, end) + plt.ylim(params[cf.YMIN_CLUTTER], params[cf.YMAX_CLUTTER]) + plt.xlabel('Date') + plt.ylabel('Average Clutter (dB)') + plt.legend(loc=1) + plt.grid(True) + plt.title('Average Clutter at site %s' % name) + for label in ax.get_xticklabels(): + label.set_rotation(90) + + # save PNG file + filename = params[cf.OUT_DIR] + "/clutter_" + name + ".png" + fig.savefig(filename, dpi=300, bbox_inches='tight') + + # avoid "RuntimeWarning: More than 20 figures have been opened" + # by closing all open figures + plt.close('all') + + return + + +def plot_scr2(t_a, t_d, scr_a, scr_d, start, end, name, params): + '''Plot RCS time series''' + fig = plt.figure() + ax = fig.add_subplot(1,1,1) + plt.plot(t_a, scr_a, 'ro-', label='Ascending') + plt.plot(t_d, scr_d, 'bo-', label='Descending') + plt.xlim(start, end) + plt.ylim(0, params[cf.YMAX_SCR]) + plt.xlabel('Date') + plt.ylabel('SCR (dB)') + plt.legend(loc=4) + plt.grid(True) + plt.title('Signal to Clutter Ratio at site %s' % name) + for label in ax.get_xticklabels(): + label.set_rotation(90) + + # save PNG file + filename = params[cf.OUT_DIR] + "/scr_" + name + ".png" + fig.savefig(filename, dpi=300, bbox_inches='tight') + + # avoid "RuntimeWarning: More than 20 figures have been opened" + # by closing all open figures + plt.close('all') + + return + + +def plot_rcs2(t_a, t_d, rcs_a, rcs_d, start, end, name, params): + '''Plot RCS time series''' + fig = plt.figure() + ax = fig.add_subplot(1,1,1) + plt.plot(t_a, rcs_a, 'ro-', label='Ascending') + plt.plot(t_d, rcs_d, 'bo-', label='Descending') + plt.xlim(start, end) + plt.ylim(0, params[cf.YMAX_RCS]) + plt.xlabel('Date') + plt.ylabel('RCS (dB$\mathregular{m^2}$)') + plt.legend(loc=4) + plt.grid(True) + plt.title('Radar Cross Section at site %s' % name) + for label in ax.get_xticklabels(): + label.set_rotation(90) + + # save PNG file + filename = params[cf.OUT_DIR] + "/rcs_" + name + ".png" + fig.savefig(filename, dpi=300, bbox_inches='tight') + + # avoid "RuntimeWarning: More than 20 figures have been opened" + # by closing all open figures + plt.close('all') + + return + diff --git a/data/clutter_SERF.png b/data/clutter_SERF.png deleted file mode 100644 index 20f6092..0000000 Binary files a/data/clutter_SERF.png and /dev/null differ diff --git a/data/coral_serf.conf b/data/coral_serf.conf new file mode 100644 index 0000000..290f8d1 --- /dev/null +++ b/data/coral_serf.conf @@ -0,0 +1,41 @@ +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +# CoRAL configuration file for Corner Reflector site SERF +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +# Input/Output file locations: +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +# file containing the list of asc and/or desc SAR intensity images +# leave blank, if one image geometry is not used +backscatter_data_asc: +backscatter_data_desc: data/mli_desc.list + +# files containing initial and updated radar coordinates of CRs +# leave blank, if one image geometry is not used +cr_file_asc: +cr_file_asc_new: +cr_file_desc: data/site_lon_lat_hgt_date1_date2_az_rg_initial.txt +cr_file_desc_new: data/out/site_lon_lat_hgt_date1_date2_az_rg.txt + +# output directory +path_out: data/out + +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +# Parameters: +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +# processing parameters +# target window size (pixels, default for multilooked data: 3) +target_window_size: 3 +# clutter window size (pixels, default for multilooked data: 11) +clutter_window_size: 11 +# cropped image size used for plotting (pixels, default for multilooked data: 51) +sub_image_size: 51 +# number of jobs for parallel processing (default: 16) +n_jobs: 16 +# plotting parameters (y-axis limits in dB) +ymax_rcs: 35 +ymax_scr: 30 +ymin_clutter: -16 +ymax_clutter: -2 \ No newline at end of file diff --git a/data/coral_serf2.conf b/data/coral_serf2.conf new file mode 100644 index 0000000..5f63477 --- /dev/null +++ b/data/coral_serf2.conf @@ -0,0 +1,41 @@ +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +# CoRAL configuration file for Corner Reflector site SERF +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +# Input/Output file locations: +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +# file containing the list of asc and/or desc SAR intensity images +# leave blank, if one image geometry is not used +backscatter_data_asc: data/mli_asc_mockup.list +backscatter_data_desc: data/mli_desc.list + +# files containing initial and updated radar coordinates of CRs +# leave blank, if one image geometry is not used +cr_file_asc: data/site_lon_lat_hgt_date1_date2_az_rg_initial_Asc.txt +cr_file_asc_new: data/out2/site_lon_lat_hgt_date1_date2_az_rg_Asc.txt +cr_file_desc: data/site_lon_lat_hgt_date1_date2_az_rg_initial.txt +cr_file_desc_new: data/out2/site_lon_lat_hgt_date1_date2_az_rg.txt + +# output directory +path_out: data/out2 + +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +# Parameters: +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +# processing parameters +# target window size (pixels, default for multilooked data: 3) +target_window_size: 3 +# clutter window size (pixels, default for multilooked data: 11) +clutter_window_size: 11 +# cropped image size used for plotting (pixels, default for multilooked data: 51) +sub_image_size: 51 +# number of jobs for parallel processing (default: 16) +n_jobs: 16 +# plotting parameters (y-axis limits in dB) +ymax_rcs: 35 +ymax_scr: 30 +ymin_clutter: -16 +ymax_clutter: -2 \ No newline at end of file diff --git a/data/mean_intensity_SERF.png b/data/mean_intensity_SERF.png deleted file mode 100644 index 60cef9e..0000000 Binary files a/data/mean_intensity_SERF.png and /dev/null differ diff --git a/data/mli_asc_mockup.list b/data/mli_asc_mockup.list new file mode 100644 index 0000000..2372343 --- /dev/null +++ b/data/mli_asc_mockup.list @@ -0,0 +1,9 @@ +./data/20180726_VV.tif +./data/20180807_VV.tif +./data/20180819_VV.tif +./data/20180831_VV.tif +./data/20180912_VV.tif +./data/20180924_VV.tif +./data/20181006_VV.tif +./data/20181018_VV.tif +./data/20181030_VV.tif diff --git a/data/mli_desc.list b/data/mli_desc.list new file mode 100644 index 0000000..c681b5a --- /dev/null +++ b/data/mli_desc.list @@ -0,0 +1,9 @@ +./data/20180726_VV.mli +./data/20180807_VV.mli +./data/20180819_VV.mli +./data/20180831_VV.mli +./data/20180912_VV.mli +./data/20180924_VV.mli +./data/20181006_VV.mli +./data/20181018_VV.mli +./data/20181030_VV.mli diff --git a/data/rcs_scr_SERF.png b/data/rcs_scr_SERF.png deleted file mode 100644 index 11c695c..0000000 Binary files a/data/rcs_scr_SERF.png and /dev/null differ diff --git a/data/site_az_rg.txt b/data/site_az_rg.txt new file mode 100644 index 0000000..79fa383 --- /dev/null +++ b/data/site_az_rg.txt @@ -0,0 +1 @@ +SERF 110 87 \ No newline at end of file diff --git a/data/site_lon_lat_hgt_date1_date2_az_rg.txt b/data/site_lon_lat_hgt_date1_date2_az_rg.txt new file mode 100644 index 0000000..aba4d6b --- /dev/null +++ b/data/site_lon_lat_hgt_date1_date2_az_rg.txt @@ -0,0 +1 @@ +SERF 123.45 -12.34 67.89 20180701 99999999 110 87 \ No newline at end of file diff --git a/data/site_lon_lat_hgt_date1_date2_az_rg_initial.txt b/data/site_lon_lat_hgt_date1_date2_az_rg_initial.txt new file mode 100644 index 0000000..5a8c9ec --- /dev/null +++ b/data/site_lon_lat_hgt_date1_date2_az_rg_initial.txt @@ -0,0 +1 @@ +SERF 123.45 -12.34 67.89 20180701 99999999 111 88 \ No newline at end of file diff --git a/data/site_lon_lat_hgt_date1_date2_az_rg_initial_Asc.txt b/data/site_lon_lat_hgt_date1_date2_az_rg_initial_Asc.txt new file mode 100644 index 0000000..5a8c9ec --- /dev/null +++ b/data/site_lon_lat_hgt_date1_date2_az_rg_initial_Asc.txt @@ -0,0 +1 @@ +SERF 123.45 -12.34 67.89 20180701 99999999 111 88 \ No newline at end of file diff --git a/tests/ref_images/clutter_SERF.png b/tests/ref_images/clutter_SERF.png new file mode 100644 index 0000000..92e6a52 Binary files /dev/null and b/tests/ref_images/clutter_SERF.png differ diff --git a/tests/ref_images/clutter_SERF2.png b/tests/ref_images/clutter_SERF2.png new file mode 100644 index 0000000..35217d0 Binary files /dev/null and b/tests/ref_images/clutter_SERF2.png differ diff --git a/tests/ref_images/mean_intensity_SERF.png b/tests/ref_images/mean_intensity_SERF.png new file mode 100644 index 0000000..1eb7aa8 Binary files /dev/null and b/tests/ref_images/mean_intensity_SERF.png differ diff --git a/tests/ref_images/mean_intensity_SERF2.png b/tests/ref_images/mean_intensity_SERF2.png new file mode 100644 index 0000000..dd6a18b Binary files /dev/null and b/tests/ref_images/mean_intensity_SERF2.png differ diff --git a/tests/ref_images/rcs_SERF.png b/tests/ref_images/rcs_SERF.png new file mode 100644 index 0000000..c148306 Binary files /dev/null and b/tests/ref_images/rcs_SERF.png differ diff --git a/tests/ref_images/rcs_SERF2.png b/tests/ref_images/rcs_SERF2.png new file mode 100644 index 0000000..6f16ddb Binary files /dev/null and b/tests/ref_images/rcs_SERF2.png differ diff --git a/tests/ref_images/scr_SERF.png b/tests/ref_images/scr_SERF.png new file mode 100644 index 0000000..bc433ad Binary files /dev/null and b/tests/ref_images/scr_SERF.png differ diff --git a/tests/ref_images/scr_SERF2.png b/tests/ref_images/scr_SERF2.png new file mode 100644 index 0000000..1495e36 Binary files /dev/null and b/tests/ref_images/scr_SERF2.png differ diff --git a/tests/test_coral.py b/tests/test_coral.py index f1a1747..e4a17b6 100644 --- a/tests/test_coral.py +++ b/tests/test_coral.py @@ -5,7 +5,9 @@ import os.path import numpy as np from coral.corner_reflector import * -from coral.dataio import readpar, readmli +from coral.dataio import readpar, readmli, read_radar_coords, write_radar_coords, read_input_files +from coral import config as cf + class TestCoral(unittest.TestCase): @classmethod @@ -123,7 +125,6 @@ def setUpClass(cls): cls.files = [] for file in glob.glob("data/*.mli"): cls.files.append(file) - cls.files.sort() cls.cr = [ 87, 110] cls.sub_im = 51 @@ -132,7 +133,7 @@ def setUpClass(cls): def test_loop(self): '''test the calculation loop function''' - avgI, rcs, scr, Avg_clt, t = loop(self.files, self.sub_im, self.cr, self.targ_win_sz, self.clt_win_sz) + avgI, rcs, scr, Avg_clt, t, cr_new, cr_pos = loop(self.files, self.sub_im, self.cr, self.targ_win_sz, self.clt_win_sz) # test the mean value of the average intensity image self.assertEqual(round(np.mean(avgI), 6), -11.385613) # -11.385613198856245 @@ -153,10 +154,89 @@ def setUpClass(cls): def test_loop(self): '''test the calculation loop function''' - avgI, rcs, scr, Avg_clt, t = loop(self.files, self.sub_im, self.cr, self.targ_win_sz, self.clt_win_sz) + avgI, rcs, scr, Avg_clt, t, cr_new, cr_pos = loop(self.files, self.sub_im, self.cr, self.targ_win_sz, self.clt_win_sz) # test the mean value of the average intensity image self.assertEqual(round(np.mean(avgI), 6), -11.385613) # -11.385613198856245 + + +class TestShift(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.files = [] + for file in glob.glob("data/*.mli"): + cls.files.append(file) + + cls.files.sort() + cls.sub_im = 51 + cls.targ_win_sz = 5 + cls.clt_win_sz = 9 + cls.cr = [ 86, 111] # coordinates changed to check shift + + def test_loop(self): + '''test the coordinate shift''' + avgI, rcs, scr, Avg_clt, t, cr_new, cr_pos = loop(self.files, self.sub_im, self.cr, self.targ_win_sz, self.clt_win_sz) + + # test array containing values of shift to be applied to radar coordinates + self.assertEqual(cr_pos[0], 52) + self.assertEqual(cr_pos[1], 50) + + +class TestConfigFile(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.file_in = "data/coral_serf.conf" + cls.sites = {'SERF': np.array([[-999, -999], [88, 111]])} + + def test_config_file(self): + '''test function to read parameters from the config-file''' + params = cf.get_config_params(self.file_in) + self.assertEqual(params[cf.SUB_IM], 51) + self.assertEqual(params[cf.ASC_LIST], None) + self.assertEqual(params[cf.DESC_LIST], 'data/mli_desc.list') + + def test_read_input_files(self): + '''test function to read input files from config-file''' + params = cf.get_config_params(self.file_in) + files_a, files_d, sites = read_input_files(params) + self.assertEqual(files_a, None) + self.assertEqual(files_d[0], './data/20180726_VV.mli') + self.assertEqual(files_d[8], './data/20181030_VV.mli') + self.assertEqual(self.sites.keys(), sites.keys()) + array1 = sites.get('SERF') + array2 = self.sites.get('SERF') + np.testing.assert_array_equal(array1, array2) + + +class TestCRfiles(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.file_in1 = "data/site_az_rg.txt" + cls.file_in2 = "data/site_lon_lat_hgt_date1_date2_az_rg.txt" + cls.file_out1 = "data/site_az_rg_new.txt" + cls.file_out2 = "data/site_lon_lat_hgt_date1_date2_az_rg_new.txt" + cls.sites = {'SERF' : np.array([[-999, -999],[ 87, 110]])} + cls.geom = "desc" + + def test_cr_file(self): + '''test function to read and write the radar coordinate files''' + # open reduced CR coordinate file + site, az, rg = read_radar_coords(self.file_in1) + self.assertEqual(int(rg[0]), 87) + self.assertEqual(int(az[0]), 110) + # write to new file + write_radar_coords(self.file_in1, self.file_out1, self.sites, self.geom) + assert os.path.exists(self.file_out1) == 1 + os.remove(self.file_out1) + + # open full CR coordinate file + site, az, rg = read_radar_coords(self.file_in2) + self.assertEqual(int(rg[0]), 87) + self.assertEqual(int(az[0]), 110) + # write to new file + write_radar_coords(self.file_in2, self.file_out2, self.sites, self.geom) + assert os.path.exists(self.file_out2) == 1 + os.remove(self.file_out2) if __name__ == '__main__': diff --git a/tests/test_plotting.py b/tests/test_plotting.py new file mode 100644 index 0000000..c80592f --- /dev/null +++ b/tests/test_plotting.py @@ -0,0 +1,155 @@ +""" +This Python module contains unittests for testing CoRAL algorithms and workflow +""" +import unittest, glob +import os.path +import numpy as np +import cv2 # requires opencv, e.g. use: pip install --user opencv-python-headless +from coral.corner_reflector import * +from coral.plot import * +from coral.plot2 import * + + +class TestPlotting(unittest.TestCase): + @classmethod + def setUpClass(cls): + file_in = "data/coral_serf.conf" + cls.name = 'SERF' + cls.geom = 'Descending' + cr = [88, 111] # this is used as initial coordinate + files = [] + for file in glob.glob("data/*.mli"): + files.append(file) + files.sort() + cls.params = cf.get_config_params(file_in) + cls.params[cf.OUT_DIR] = './tests/' + cls.avgI, cls.rcs, cls.scr, cls.clt, cls.t, cls.cr_new, cls.cr_pos = loop(files, cls.params[cf.SUB_IM], cr, \ + cls.params[cf.TARG_WIN_SZ], cls.params[cf.CLT_WIN_SZ]) + start_time = cls.t[0] + end_time = cls.t[-1] + margin = (end_time - start_time) / 50 + cls.start = start_time - margin + cls.end = end_time + margin + + def test_plot_mean_intensity(self): + '''test function to plot mean intensity (single geometry)''' + plot_mean_intensity(self.avgI, self.cr_pos, self.name, self.params) + ref_image = './tests/ref_images/mean_intensity_SERF.png' + act_image = './tests/mean_intensity_SERF.png' + imageA = cv2.imread(ref_image) + imageB = cv2.imread(act_image) + err = np.sum((imageA.astype("float") - imageB.astype("float")) ** 2) # squared sum of differences + self.assertEqual(err, 0.0) + os.remove(act_image) + + def test_plot_clutter(self): + '''test function to plot clutter time series (single geometry)''' + plot_clutter(self.t, self.clt, self.start, self.end, self.name, self.geom, self.params) + ref_image = './tests/ref_images/clutter_SERF.png' + act_image = './tests/clutter_SERF.png' + imageA = cv2.imread(ref_image) + imageB = cv2.imread(act_image) + err = np.sum((imageA.astype("float") - imageB.astype("float")) ** 2) # squared sum of differences + self.assertEqual(err, 0.0) + os.remove(act_image) + + def test_plot_scr(self): + '''test function to plot SCR (single geometry)''' + plot_scr(self.t, self.scr, self.start, self.end, self.name, self.geom, self.params) + ref_image = './tests/ref_images/scr_SERF.png' + act_image = './tests/scr_SERF.png' + imageA = cv2.imread(ref_image) + imageB = cv2.imread(act_image) + err = np.sum((imageA.astype("float") - imageB.astype("float")) ** 2) # squared sum of differences + self.assertEqual(err, 0.0) + os.remove(act_image) + + def test_plot_rcs(self): + '''test function to plot RCS (single geometry)''' + plot_rcs(self.t, self.rcs, self.start, self.end, self.name, self.geom, self.params) + ref_image = './tests/ref_images/rcs_SERF.png' + act_image = './tests/rcs_SERF.png' + imageA = cv2.imread(ref_image) + imageB = cv2.imread(act_image) + err = np.sum((imageA.astype("float") - imageB.astype("float")) ** 2) # squared sum of differences + self.assertEqual(err, 0.0) + os.remove(act_image) + + +class TestPlotting2(unittest.TestCase): + @classmethod + def setUpClass(cls): + file_in = "data/coral_serf2.conf" + cls.name = 'SERF' + cr = [88, 111] # this is used as initial coordinate + files_a = [] + for file in glob.glob("data/*.tif"): + files_a.append(file) + files_a.sort() + files_d = [] + for file in glob.glob("data/*.mli"): + files_d.append(file) + files_d.sort() + cls.params = cf.get_config_params(file_in) + cls.params[cf.OUT_DIR] = './tests/' + cls.avgI_a, cls.rcs_a, cls.scr_a, cls.clt_a, cls.t_a, cls.cr_new_a, cls.cr_pos_a = loop(files_a, \ + cls.params[cf.SUB_IM], cr, \ + cls.params[cf.TARG_WIN_SZ], cls.params[cf.CLT_WIN_SZ]) + cls.avgI_d, cls.rcs_d, cls.scr_d, cls.clt_d, cls.t_d, cls.cr_new_d, cls.cr_pos_d = loop(files_d, \ + cls.params[cf.SUB_IM], cr, \ + cls.params[cf.TARG_WIN_SZ], cls.params[cf.CLT_WIN_SZ]) + start_time = min(cls.t_a[0], cls.t_d[0]) + end_time = max(cls.t_a[-1], cls.t_d[-1]) + margin = (end_time - start_time) / 50 + cls.start = start_time - margin + cls.end = end_time + margin + + def test_plot_mean_intensity2(self): + '''test function to plot mean intensity (both geometries)''' + plot_mean_intensity2(self.avgI_a, self.avgI_d, self.cr_pos_a, self.cr_pos_d, self.name, self.params) + ref_image = './tests/ref_images/mean_intensity_SERF2.png' + act_image = './tests/mean_intensity_SERF.png' + imageA = cv2.imread(ref_image) + imageB = cv2.imread(act_image) + err = np.sum((imageA.astype("float") - imageB.astype("float")) ** 2) # squared sum of differences + self.assertEqual(err, 0.0) + os.remove(act_image) + + def test_plot_clutter2(self): + '''test function to plot clutter time series (single geometry)''' + plot_clutter2(self.t_a, self.t_d, self.clt_a, self.clt_d, self.start, self.end, self.name, self.params) + ref_image = './tests/ref_images/clutter_SERF2.png' + act_image = './tests/clutter_SERF.png' + imageA = cv2.imread(ref_image) + imageB = cv2.imread(act_image) + print(imageA.shape) + print(imageB.shape) + err = np.sum((imageA.astype("float") - imageB.astype("float")) ** 2) # squared sum of differences + self.assertEqual(err, 0.0) + os.remove(act_image) + + def test_plot_scr2(self): + '''test function to plot SCR (single geometry)''' + plot_scr2(self.t_a, self.t_d, self.scr_a, self.scr_d, self.start, self.end, self.name, self.params) + ref_image = './tests/ref_images/scr_SERF2.png' + act_image = './tests/scr_SERF.png' + imageA = cv2.imread(ref_image) + imageB = cv2.imread(act_image) + err = np.sum((imageA.astype("float") - imageB.astype("float")) ** 2) # squared sum of differences + self.assertEqual(err, 0.0) + os.remove(act_image) + + def test_plot_rcs2(self): + '''test function to plot RCS (single geometry)''' + plot_rcs2(self.t_a, self.t_d, self.rcs_a, self.rcs_d, self.start, self.end, self.name, self.params) + ref_image = './tests/ref_images/rcs_SERF2.png' + act_image = './tests/rcs_SERF.png' + imageA = cv2.imread(ref_image) + imageB = cv2.imread(act_image) + err = np.sum((imageA.astype("float") - imageB.astype("float")) ** 2) # squared sum of differences + self.assertEqual(err, 0.0) + os.remove(act_image) + + +if __name__ == '__main__': + unittest.main() diff --git a/utils/get_CR_position_rdc.py b/utils/get_CR_position_rdc.py new file mode 100644 index 0000000..334187f --- /dev/null +++ b/utils/get_CR_position_rdc.py @@ -0,0 +1,215 @@ +""" +Main script +- converts lat/lon coordinates of corner reflectors into azimuth and range + pixel positions +- can be run after SLC images have been coregistered to a DEM and + lat/lon coordinates have been derived (using calc_lat_lon_TxxxX) + +INPUT: +txtfile: CR_site_lon_lat_hgt_date.txt located in ../GAMMA/TxxxX/ + *sar_latlon.txt file containing lon/lat values, located in ../GAMMA/TxxxX/DEM + ./DEM/diff*.par file containing length/width of radar-coded DEM, located in ../GAMMA/TxxxX/DEM + +OUTPUT: +textfile: CR_site_lon_lat_hgt_date_az_rg.txt + same as Input file with two columns added containing + azimuth and range position of central pixel at CR location + +@author: Thomas Fuhrmann @ GA June, 2018 + +usage: +this script is currently located and executed in the TxxxX directory +Load python3 module on NCI: module load python3 +then execute, e.g.: python3 get_CR_position_rdc.py +for effective parallel processing and if memory is exceeded, start an interactive session: +qsub -I -q express -lstorage=gdata/dg9,walltime=1:00:00,mem=32Gb,ncpus=7,wd +""" + + +# import modules +import sys +import os +import os.path +import math +import fnmatch +import numpy as np +# for parallel processing: module joblib has to be installed +from joblib import Parallel, delayed +import multiprocessing + + +####################### +### Path and file names +####################### +# give general processing path here +io_path = "/g/data/dg9/INSAR_ANALYSIS/CAMDEN/S1/GAMMA/T147D/" +# note that the GAMMA "DEM" folder is assumed inside this path +dem_path = io_path + "DEM/" +# give name of CR input file consisting of site name, longitude, latitude (+ optional column fields) +cr_filename = "CR_site_lon_lat_hgt_date1_date2.txt" +# name of CR output file (automatic naming): +cr_filename_out = cr_filename[:-4] + "_az_rg_initial.txt" +######## +# set number of processors here: +nCPU = 7 +######## + + + +# Welome +print() +print("########################################") +print("# Get CR position in radar-coordinates #") +print("########################################") + +print() + + +# Read CR_site_lon_lat_hgt_date1_date2.txt +print("Reading textfile with CR positions (latitude longitude height date)\ +...") +filename_in = io_path + cr_filename +if os.path.isfile(filename_in) and os.access(filename_in, os.R_OK): + print("File", filename_in, "exists and is readable.") + f = open(filename_in) + sites = f.readlines() + sites[:] = [line.split("\t")[0] for line in sites] + f = open(filename_in) + lons = f.readlines() + lons[:] = [float(line.split("\t")[1]) for line in lons] + f = open(filename_in) + lats = f.readlines() + lats[:] = [float(line.split("\t")[2]) for line in lats] +else: + print("ERROR: Can't read file", filename) + sys.exit() +print() + + +# Read diff.par file +for file in os.listdir(dem_path): + if fnmatch.fnmatch(file, "diff*.par"): + filename = dem_path + file + if os.path.isfile(filename) and os.access(filename, os.R_OK): + f = open(filename) + lines = f.readlines() + for line in lines: + if "range_samp_1" in line: + width = int(line.split()[1]) + if "az_samp_1" in line: + length = int(line.split()[1]) + else: + print("ERROR: Can't read file", filename) + sys.exit() +print("Width of radar-coded data (i.e. number of range samples):", width) +print("Length of radar-coded data (i.e. number of azimuth lines):", length) +print() + + +# Read lat/lon coordinate file +for file in os.listdir(dem_path): + if fnmatch.fnmatch(file, "*sar_latlon.txt"): + filename = dem_path + file + if os.path.isfile(filename) and os.access(filename, os.R_OK): + print("File", filename, "exists and is readable.") + else: + print("ERROR: Can't read file", filename) + sys.exit() +print() + + +# Convert lat/lon diff to metric using the CR lat/lon values +# Constants +a = 6378137 +f = 0.0033528107 +e2 = 1/(1-f)**2-1 +c = a*math.sqrt(1+e2) + + +def _inner(cr): + lonCR = cr[0][0] + latCR = cr[0][1] + f = open(filename) + lines = f.readlines() + dist_sq_min = 999 # in deg + az = 0 + rg = 0 + for line in lines: + lat = float(line.split(" ")[0]) + lon = float(line.split(" ")[1]) + dist_sq = (lon-lonCR)**2 + (lat-latCR)**2 + if dist_sq < dist_sq_min: + dist_sq_min = dist_sq + az_CR = az + rg_CR = rg + if rg < width-1: + rg = rg + 1 + else: + rg = 0 + az = az + 1 + # save az and rg coordinates to list + #az_CRs.append(az_CR) + #rg_CRs.append(rg_CR) + # print distance of pixel coordinate to CR coordinate + dist_min = math.sqrt(dist_sq_min) + V = math.sqrt((1+(e2*math.cos(latCR/180*math.pi)**2))) + N = c/V + M = c/V**3 + dist_min_m = dist_min/180*math.pi*math.sqrt(M*N) + print("Pixel with minimum distance to CR location is:") + print("azimuth: %d, range: %d" % (az_CR, rg_CR)) + print("Distance is %8.6f degree or %3.1f m" % (dist_min, dist_min_m)) + return az_CR, rg_CR +print() + + +################################## +# function for parallel loop processing +keys = sites +values = (np.array([[lons[i], lats[i]]], dtype=float) \ + for i in range(0,len(sites))) +sites = dict(zip(keys, values)) +names = sites.keys() +# note that dictionaries are unsorted and the variable name is hence not ordered +def processInput(name): + cr = sites[name] + az, rg = _inner(cr) + return az, rg + +# 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=nCPU)(delayed(processInput)(name) for name in names) +az_CRs = [] +rg_CRs = [] +for i in range(0,len(names)): + az_CRs.append(results[i][0]) + rg_CRs.append(results[i][1]) +# print(az_CRs) +# print(rg_CRs) + +# save the range and azimuth numbers of pixels in a new txt file +print("Reading textfile with CR positions (latitude longitude height date)\ +...") +filename_out = io_path + cr_filename_out +fout = open(filename_out,'w') +if os.path.isfile(filename_in) and os.access(filename_in, os.R_OK): + print("File", filename_in, "exists and is readable.") + f = open(filename_in) + lines = f.readlines() + idx = 0 + for line in lines: + out_line = line.strip("\n") + "\t" + str(az_CRs[idx]) + "\t" + \ + str(rg_CRs[idx]) + "\n" + fout.write(out_line) + idx = idx + 1 +else: + print("ERROR: Can't read file", filename_in) + sys.exit() +fout.close() +print() +print("Azimuth and Range number of CR pixels has been written to " + \ +filename_out + ".") +print() +