Merge remote-tracking branch 'origin/astrohack-dev' into astrohack-dev

# Conflicts:
#	pyproject.toml

src/astrohack/_utils/_algorithms.py

@@ -6,6 +6,33 @@
 
 from astrohack._utils._logger._astrohack_logger import _get_astrohack_logger
 
+def _calculate_suggested_grid_paramater(parameter, quantile=0.01):
+    import scipy
+
+    logger = _get_astrohack_logger()
+
+    # Determine skew properties and return median
+    if np.abs(scipy.stats.skew(parameter)) > 0.5:
+        if scipy.stats.skew(parameter) > 0:
+            cutoff = np.quantile(parameter, 1 - quantile)
+            parameter = parameter[parameter <= cutoff]
+
+        else:
+            cutoff = np.quantile(parameter, quantile)
+            parameter = parameter[parameter >= cutoff]
+
+        return np.median(parameter)
+
+    # Process as mean
+    else:
+        upper_cutoff = np.quantile(parameter, 1 - quantile)
+        lower_cutoff = np.quantile(parameter, quantile)
+
+        parameter = parameter[parameter >= lower_cutoff]
+        parameter = parameter[parameter <= upper_cutoff]    
+
+        return np.mean(parameter)
+
 def _apply_mask(data, scaling=0.5):
     """ Applies a cropping mask to the input data according to the scale factor
 
@@ -225,8 +252,8 @@ def _get_grid_params(vis_map_dict, pnt_map_dict, ant_names):
         n_pix_x = np.sum([abs_diff[:,0] > max_dis_x]) + 1
         n_pix_y = np.sum([abs_diff[:,1] > max_dis_y]) + 1
 
-        cell_size_x = np.mean(abs_diff[abs_diff[:,0] > max_dis_x,0])
-        cell_size_y = np.mean(abs_diff[abs_diff[:,1] > max_dis_y,1])
+        cell_size_x = np.mean(abs_diff[abs_diff[:,0] > max_dis_x, 0])
+        cell_size_y = np.mean(abs_diff[abs_diff[:,1] > max_dis_y, 1])
 
         if n_pix_x < n_pix_y:
             n_pix = n_pix_x**2

src/astrohack/_utils/_extract_holog.py

@@ -4,6 +4,8 @@
 import numpy as np
 import xarray as xr
 import astropy
+import astrohack
+
 from numba import njit
 from numba.core import types
 
@@ -386,11 +388,14 @@ def _create_holog_file(
             )
 
 
-def _create_holog_obs_dict(pnt_dict, baseline_average_distance, baseline_average_nearest, ant_names, ant_pos, ant_names_main):
+def _create_holog_obs_dict(pnt_dict, baseline_average_distance, baseline_average_nearest, ant_names, ant_pos, ant_names_main, write_distance_matrix=False):
     '''
     Generate holog_obs_dict.
     '''
 
+    import pandas as pd
+    from scipy.spatial import distance_matrix
+
     logger = _get_astrohack_logger()
     mapping_scans_dict = {}
     holog_obs_dict = {}
@@ -419,17 +424,13 @@ def _create_holog_obs_dict(pnt_dict, baseline_average_distance, baseline_average
                                 holog_obs_dict[ddi][map_key] = {'scans':np.array(scan_list),'ant':{}}
 
                             holog_obs_dict[ddi][map_key]['ant'][ant_name] = []
-
-    # If users specifies a baseline_average_distance we need to create an antenna distance matrix.
-    if (baseline_average_distance != 'all') or (baseline_average_nearest != 'all'):
-
-        import pandas as pd
-        from scipy.spatial import distance_matrix
-
-        df = pd.DataFrame(ant_pos, columns=['x', 'y', 'z'], index=ant_names)
-        df_mat = pd.DataFrame(distance_matrix(df.values, df.values), index=df.index, columns=df.index)
-        logger.debug('Antenna distance matrix in meters: \n' + str(df_mat))
 
+    df = pd.DataFrame(ant_pos, columns=['x', 'y', 'z'], index=ant_names)
+    df_mat = pd.DataFrame(distance_matrix(df.values, df.values), index=df.index, columns=df.index)
+    if write_distance_matrix:
+        df_mat.to_csv(path_or_buf="{base}/.baseline_distance_matrix.csv".format(base=os.getcwd()), sep="\t")
+        logger.info("Writing distance matrix to {base}/.baseline_distance_matrix.csv ...".format(base=os.getcwd()))
+
 
     if (baseline_average_distance != 'all') and (baseline_average_nearest != 'all'):
         logger.error('[{function_name}]: baseline_average_distance and baseline_average_nearest can not both be specified.'.format(function_name=function_name))
@@ -556,14 +557,22 @@ def _create_holog_meta_data(holog_file, holog_dict, input_params):
         'n_pix': n_pixs[0],
         'telescope_name': telescope_names[0]
     }
-
+    
     if not (len(set(cell_sizes_sigfigs)) == 1):
-        logger.error('Cell size not consistent: ' + str(cell_sizes))
-        meta_data['cell_size'] = None
+        logger.warning('Cell size not consistent: ' + str(cell_sizes))
+        logger.warning('Calculating suggested cell size ...')
+
+        meta_data["cell_size"] = astrohack._utils._algorithms._calculate_suggested_grid_paramater(parameter=np.array(cell_sizes))
+
+        logger.info("The suggested cell size is calculated to be: {cell_size}".format(cell_size=meta_data["cell_size"]))
 
     if not (len(set(n_pixs)) == 1):
-        logger.error('Number of pixels not consistent: ' + str(n_pixs))
-        meta_data['n_pix'] = None
+        logger.warning('Number of pixels not consistent: ' + str(n_pixs))
+        logger.warning('Calculating suggested number of pixels ...')
+
+        meta_data['n_pix'] = int(astrohack._utils._algorithms._calculate_suggested_grid_paramater(parameter=np.array(n_pixs)))
+
+        logger.info("The suggested number of pixels is calculated to be: {n_pix}".format(n_pix=meta_data["n_pix"]))
 
     if not (len(set(telescope_names)) == 1):
         logger.error('Telescope name not consistent: ' + str(telescope_names))

src/astrohack/_utils/_extract_locit.py

@@ -1,6 +1,4 @@
 import numpy as np
-from matplotlib import pyplot as plt
-from matplotlib.patches import Rectangle
 
 from casacore import tables as ctables
 from astropy.coordinates import SkyCoord, CIRS
@@ -12,9 +10,11 @@
 from astrohack._utils._logger._astrohack_logger import _get_astrohack_logger
 from astrohack._utils._tools import _casa_time_to_mjd, _rad_to_deg_str
 from astrohack._utils._conversion import _convert_unit
-from astrohack._utils._constants import figsize, twopi, fontsize, notavail
+from astrohack._utils._constants import figsize, twopi, notavail
 from astrohack._utils._dio import _write_meta_data
-from astrohack._utils._locit import _open_telescope
+from astrohack._utils._locit_commons import _open_telescope, _compute_antenna_relative_off, _get_telescope_lat_lon_rad
+from astrohack._utils._locit_commons import _create_figure_and_axes, _plot_boxes_limits_and_labels
+from astrohack._utils._locit_commons import _plot_antenna_position, _close_figure, _scatter_plot
 
 
 def _extract_antenna_data(fname, extract_locit_params):
@@ -287,7 +287,17 @@ def _extract_antenna_phase_gains(fname, extract_locit_params):
 
 def _plot_source_table(filename, src_dict, label=True, precessed=False, obs_midpoint=None, display=True,
                        figure_size=figsize, dpi=300):
-    """ Backend function for plotting the source table"""
+    """ Backend function for plotting the source table
+    Args:
+        filename: Name for the png plot file
+        src_dict: The dictionary containing the observed sources
+        label: Add source labels
+        precessed: Plot sources with precessed coordinates
+        obs_midpoint: Time to which precesses the coordiantes
+        display: Display plots in matplotlib
+        figure_size: plot dimensions in inches
+        dpi: Dots per inch (plot resolution)
+    """
     logger = _get_astrohack_logger()
     n_src = len(src_dict)
     radec = np.ndarray((n_src, 2))
@@ -308,35 +318,32 @@ def _plot_source_table(filename, src_dict, label=True, precessed=False, obs_midp
         radec[int(i_src)] = src[coorkey]
         name.append(src['name'])
 
-    if figure_size is None or figure_size == 'None':
-        fig, ax = plt.subplots(1, 1, figsize=figsize)
-    else:
-        fig, ax = plt.subplots(1, 1, figsize=figure_size)
+    fig, ax  = _create_figure_and_axes(figure_size, [1, 1])
     radec[:, 0] *= _convert_unit('rad', 'hour', 'trigonometric')
     radec[:, 1] *= _convert_unit('rad', 'deg', 'trigonometric')
 
+    xlabel = 'Right Ascension [h]'
+    ylabel = 'Declination [\u00b0]'
     if label:
-        for i_src in range(n_src):
-            ax.plot(radec[i_src, 0], radec[i_src, 1], marker='+', ls='', color='red')
-            ax.text(radec[i_src, 0]+0.05, radec[i_src, 1], name[i_src], fontsize=.8*fontsize, ha='left', va='center',
-                    rotation=20)
+        labels = name
     else:
-        ax.plot(radec[:, 0], radec[:, 1], marker='+', ls='', color='red')
-    ax.set_xlim([-0.5, 24.5])
-    ax.set_ylim([-95, 95])
-    ax.set_xlabel('Right Ascension [h]')
-    ax.set_ylabel('Declination [\u00b0]')
-
-    fig.suptitle(title)
-    fig.tight_layout()
-    plt.savefig(filename, dpi=dpi)
-    if not display:
-        plt.close()
+        labels = None
+
+    _scatter_plot(ax, radec[:, 0], xlabel, radec[:, 1], ylabel, title=None, labels=labels, xlim=[-0.5, 24.5],
+                  ylim=[-95, 95])
+
+    _close_figure(fig, title, filename, dpi, display)
     return
 
 
-def _plot_antenna_table(ant_dict, telescope_name, param_dict):
-    """Plot antenna positions"""
+def _plot_array_configuration(ant_dict, telescope_name, parm_dict):
+    """ backend for plotting array configuration
+
+    Args:
+        ant_dict: Dictionary containing antenna information
+        telescope_name: Name of the telescope used in observations
+        parm_dict: Parameter dictionary crafted by the calling function
+    """
 
     telescope = _open_telescope(telescope_name)
     stations = param_dict['stations']
@@ -348,79 +355,43 @@ def _plot_antenna_table(ant_dict, telescope_name, param_dict):
     box_size = param_dict['box_size']  # In user input unit
     plot_zoff = param_dict['zoff']
 
-    if figure_size is None or figure_size == 'None':
-        fig, axes = plt.subplots(1, 2, figsize=[10, 5])
-    else:
-        fig, axes = plt.subplots(1, 2, figsize=figure_size)
+    fig, axes = _create_figure_and_axes(figure_size, [1, 2], default_figsize=[10, 5])
 
     len_fac = _convert_unit('m', length_unit, 'length')
-    half_box = box_size / 2
 
-    ax_box = axes[1]
-    ax_all = axes[0]
+    inner_ax = axes[1]
+    outer_ax = axes[0]
 
-    tel_lon = telescope.array_center['m0']['value']
-    tel_lat = telescope.array_center['m1']['value']
-    tel_rad = telescope.array_center['m2']['value']
+    tel_lon, tel_lat, tel_rad = _get_telescope_lat_lon_rad(telescope)
 
     for antenna in ant_dict.values():
-        ew_off, ns_off, el_off, _ = _compute_antenna_relative_off(antenna, tel_lon, tel_lat, tel_rad)
-        ew_off *= len_fac
-        ns_off *= len_fac
-        el_off *= len_fac
+        ew_off, ns_off, el_off, _ = _compute_antenna_relative_off(antenna, tel_lon, tel_lat, tel_rad, len_fac)
         text = f'  {antenna["name"]}'
         if stations:
             text += f'@{antenna["station"]}'
         if plot_zoff:
             text += f' {el_off:.1f} {length_unit}'
-        if abs(ew_off) > half_box or abs(ns_off) > half_box:
-            ax_all.plot(ew_off, ns_off, marker='+', color='black')
-            ax_all.text(ew_off, ns_off, text, fontsize=fontsize, ha='left', va='center')
-        else:
-            ax_box.plot(ew_off, ns_off, marker='+', color='black')
-            ax_box.text(ew_off, ns_off, text, fontsize=fontsize, ha='left', va='center')
+        _plot_antenna_position(outer_ax, inner_ax, ew_off, ns_off, text, box_size)
 
     # axes labels
     xlabel = f'East [{length_unit}]'
     ylabel = f'North [{length_unit}]'
 
-    # Larger box limits and labels
-    x_lim, y_lim = ax_all.get_xlim(), ax_all.get_ylim()
-    x_half, x_mid = (x_lim[1] - x_lim[0])/2, (x_lim[1] + x_lim[0]) / 2
-    y_half, y_mid = (y_lim[1] - y_lim[0])/2, (y_lim[1] + y_lim[0]) / 2
-    if x_half > y_half:
-        y_lim = [y_mid-x_half, y_mid+x_half]
-    else:
-        x_lim = [x_mid-y_half, x_mid+y_half]
-    ax_all.set_xlim(x_lim)
-    ax_all.set_ylim(y_lim)
-    ax_all.set_xlabel(xlabel)
-    ax_all.set_ylabel(ylabel)
-    ax_all.plot(0, 0, marker='x', color='blue')
-    box = Rectangle([-half_box, -half_box], 2*half_box, 2*half_box, linewidth=0.5, edgecolor='red', facecolor='none')
-    ax_all.add_patch(box)
-    ax_all.set_title('Whole array')
-    ax_all.set_aspect(1)
-
-    # Smaller box limits and labels
-    ax_box.set_xlim([-half_box, half_box])
-    ax_box.set_ylim([-half_box, half_box])
-    ax_box.set_xlabel(xlabel)
-    ax_box.set_ylabel(ylabel)
-    ax_box.plot(0, 0, marker='x', color='blue')
-    ax_box.set_title('Inner array')
-    ax_box.set_aspect(1)
+    _plot_boxes_limits_and_labels(outer_ax, inner_ax, xlabel, ylabel, box_size, 'Outer array', 'Inner array')
 
     title = 'Antenna positions during observation'
-    fig.suptitle(title)
-    fig.tight_layout()
-    plt.savefig(filename, dpi=dpi)
-    if not display:
-        plt.close()
+    _close_figure(fig, title, filename, dpi, display)
     return
 
 
-def _print_antenna_table(params, ant_dict, telescope_name):
+def _print_array_configuration(params, ant_dict, telescope_name):
+    """ Backend for printing the array configuration onto a table
+
+    Args:
+        params: Parameter dictionary crafted by the calling function
+        ant_dict: Parameter dictionary crafted by the calling function
+        telescope_name: Name of the telescope used in observations
+    """
     telescope = _open_telescope(telescope_name)
     relative = params['relative']
 
@@ -430,9 +401,7 @@ def _print_antenna_table(params, ant_dict, telescope_name):
     if relative:
         nfields = 5
         table.field_names = ['Name', 'Station', 'East [m]', 'North [m]', 'Elevation [m]', 'Distance [m]']
-        tel_lon = telescope.array_center['m0']['value']
-        tel_lat = telescope.array_center['m1']['value']
-        tel_rad = telescope.array_center['m2']['value']
+        tel_lon, tel_lat, tel_rad = _get_telescope_lat_lon_rad(telescope)
     else:
         nfields = 4
         table.field_names = ['Name', 'Station', 'Longitude', 'Latitude', 'Radius [m]']
@@ -461,9 +430,3 @@ def _print_antenna_table(params, ant_dict, telescope_name):
     return
 
 
-def _compute_antenna_relative_off(antenna, tel_lon, tel_lat, tel_rad):
-    antenna_off_east = tel_rad * (antenna['longitude'] - tel_lon) * np.cos(tel_lat)
-    antenna_off_north = tel_rad * (antenna['latitude'] - tel_lat)
-    antenna_off_ele = antenna['radius'] - tel_rad
-    antenna_dist = np.sqrt(antenna_off_east ** 2 + antenna_off_north ** 2 + antenna_off_ele ** 2)
-    return antenna_off_east, antenna_off_north, antenna_off_ele, antenna_dist