286 extract holog may fail when pointing offsets contain nans

src/astrohack/core/extract_holog.py

@@ -590,20 +590,34 @@ def _extract_pointing_chunk(map_ant_ids, time_vis, pnt_ant_dict):
     Returns:
         dict:  Dictionary of directional cosine data mapped to nearest MAIN table sample times.
     """
-
     pnt_map_dict = {}
     coords = {"time": time_vis}
     for antenna in map_ant_ids:
         pnt_xds = pnt_ant_dict[antenna]
-        avg_dir, avg_dir_cos, avg_enc, avg_pnt_off, avg_tgt = \
-            _time_avg_pointing_jit(time_vis,
-                                   pnt_xds.time.values,
-                                   pnt_xds['DIRECTION'].values,
-                                   pnt_xds['DIRECTIONAL_COSINES'].values,
-                                   pnt_xds['ENCODER'].values,
-                                   pnt_xds['POINTING_OFFSET'].values,
-                                   pnt_xds['TARGET'].values,
-                                   )
+        pnt_time = pnt_xds.time.values
+        pnt_int = np.average(np.diff(pnt_time))
+        vis_int = time_vis[1] - time_vis[0]
+
+        if pnt_int < vis_int:
+            avg_dir, avg_dir_cos, avg_enc, avg_pnt_off, avg_tgt = \
+                _time_avg_pointing_jit(time_vis,
+                                       pnt_xds.time.values,
+                                       pnt_xds['DIRECTION'].values,
+                                       pnt_xds['DIRECTIONAL_COSINES'].values,
+                                       pnt_xds['ENCODER'].values,
+                                       pnt_xds['POINTING_OFFSET'].values,
+                                       pnt_xds['TARGET'].values,
+                )
+        else:
+            avg_dir, avg_dir_cos, avg_enc, avg_pnt_off, avg_tgt = \
+                _interpolate_pointing(time_vis,
+                                       pnt_xds.time.values,
+                                       pnt_xds['DIRECTION'].values,
+                                       pnt_xds['DIRECTIONAL_COSINES'].values,
+                                       pnt_xds['ENCODER'].values,
+                                       pnt_xds['POINTING_OFFSET'].values,
+                                       pnt_xds['TARGET'].values,
+                )
 
         new_pnt_xds = xr.Dataset()
         new_pnt_xds.assign_coords(coords)
@@ -638,15 +652,16 @@ def _time_avg_pointing_jit(time_vis, pnt_time, dire, dir_cos, enc, pnt_off, tgt)
             continue
         else:
             i_time = _get_time_index(pnt_time[i_row], i_time, time_vis, half_int)
-        if i_time < 0:
-            break
+            if i_time < 0:
+                break
+
         avg_dir[i_time] += dire[i_row]
         avg_dir_cos[i_time] += dir_cos[i_row]
         avg_enc[i_time] += enc[i_row]
         avg_pnt_off[i_time] += pnt_off[i_row]
         avg_tgt[i_time] += tgt[i_row]
-        avg_wgt[i_time] += 1
-
+        avg_wgt[i_time] += 1            
+            
     avg_dir /= avg_wgt
     avg_dir_cos /= avg_wgt
     avg_enc /= avg_wgt
@@ -656,6 +671,34 @@ def _time_avg_pointing_jit(time_vis, pnt_time, dire, dir_cos, enc, pnt_off, tgt)
     return avg_dir, avg_dir_cos, avg_enc, avg_pnt_off, avg_tgt
 
 
+@njit(cache=False, nogil=True)
+def _interpolate_pointing(time_vis, pnt_time, dire, dir_cos, enc, pnt_off, tgt):
+    n_samples = time_vis.shape[0]
+    the_shape = (n_samples, 2)
+    n_row = pnt_time.shape[0]
+    pnt_int = np.average(np.diff(pnt_time))
+    half_int = (time_vis[1] - time_vis[0]) / 2
+
+    avg_dir = np.zeros(the_shape)
+    avg_dir_cos = np.zeros(the_shape)
+    avg_enc = np.zeros(the_shape)
+    avg_pnt_off = np.zeros(the_shape)
+    avg_tgt = np.zeros(the_shape)
+    avg_wgt = np.zeros(the_shape)
+
+    for i_time in range(n_samples):
+        i_row = int(np.floor((time_vis[i_time]-half_int-pnt_time[0])/pnt_int))
+
+        avg_dir[i_time] += dire[i_row]
+        avg_dir_cos[i_time] += dir_cos[i_row]
+        avg_enc[i_time] += enc[i_row]
+        avg_pnt_off[i_time] += pnt_off[i_row]
+        avg_tgt[i_time] += tgt[i_row]
+        avg_wgt[i_time] += 1         
+
+    return avg_dir, avg_dir_cos, avg_enc, avg_pnt_off, avg_tgt
+
+
 def create_holog_meta_data(holog_file, holog_dict, input_params):
     """Save holog file meta information to json file with the transformation
         of the ordering (ddi, holog_map, ant) --> (ant, ddi, holog_map).

src/astrohack/extract_holog.py

@@ -259,7 +259,7 @@ def extract_holog(
         holog_obs_dict: HologObsDict = None,
         ddi: Union[int, List[int], str] = 'all',
         baseline_average_distance: Union[float, str] = 'all',
-        baseline_average_nearest: Union[float, str] = 'all',
+        baseline_average_nearest: Union[float, str] = 1,
         exclude_antennas: Union[list[str], str] = None,
         data_column: str = "CORRECTED_DATA",
         time_smoothing_interval: float = None,
@@ -275,7 +275,7 @@ def extract_holog(
     :param point_name: Name of *<point_name>.point.zarr* file to use. This is must be provided.
     :type holog_name: str
 
-    :param holog_name: Name of *<holog_name>.holog.zarr* file to create. Defaults to measurement set name with
+    :param holog_name: Name of *<holog_name>.holog.zarr* file to create. Defaults to measurement set name with \
     *holog.zarr* extension.
     :type holog_name: str, optional
 
@@ -303,8 +303,8 @@ def extract_holog(
     baseline_average_nearest can not be used together.
     :type baseline_average_nearest: int, optional
 
-    :param exclude_antennas: If an antenna is given for exclusion it will not be processed as a reference or a
-    mapping antenna. This can be used to exclude antennas that have bad data for whatever reason. Default is None,
+    :param exclude_antennas: If an antenna is given for exclusion it will not be processed as a reference or a \
+    mapping antenna. This can be used to exclude antennas that have bad data for whatever reason. Default is None, \
     meaning no antenna is excluded.
     :type exclude_antennas: str | list, optional
 
@@ -529,28 +529,24 @@ def extract_holog(
     )
 
     telescope_name = obs_ctb.getcol("TELESCOPE_NAME")[0]
-    start_time_unix = obs_ctb.getcol('TIME_RANGE')[0][0] - 3506716800.0
-    time = Time(start_time_unix, format='unix').jyear
+    # start_time_unix = obs_ctb.getcol('TIME_RANGE')[0][0] - 3506716800.0
+    # time = Time(start_time_unix, format='unix').jyear
 
     # If we have an EVLA run from before 2023 the pointing table needs to be fixed.
-    if telescope_name == "EVLA" and time < 2023:
-
-        # Convert from casa epoch to unix time
-        his_ctb = ctables.table(
-            os.path.join(extract_holog_params['ms_name'], "HISTORY"),
-            readonly=True,
-            lockoptions={"option": "usernoread"},
-            ack=False,
-        )
-
-        if "pnt_tbl:fixed" not in his_ctb.getcol("MESSAGE"):
-            logger.error(
-                "Pointing table not corrected, users should apply function astrohack.dio.fix_pointing_table() to "
-                "remedy this.")
-
-            return None
-
-        his_ctb.close()
+    if telescope_name == "EVLA": # and time < 2023:
+        n_mapping = 0
+        for ddi_key, ddi_dict in holog_obs_dict.items():
+            n_map_ddi = 0
+            for map_dict in ddi_dict.values():
+                n_map_ddi += len(map_dict['ant'])
+            if n_map_ddi == 0:
+                logger.warning(f'DDI {ddi_key} has 0 mapping antennas')
+            n_mapping += n_map_ddi
+
+        if n_mapping == 0:
+            msg = 'No mapping antennas to process, maybe you need to fix the pointing table?'
+            logger.error(msg)
+            raise Exception(msg)
 
     count = 0
     delayed_list = []

src/astrohack/utils/algorithms.py

@@ -309,6 +309,11 @@ def _least_squares_fit_block(system, vector):
     return results, variances
 
 
+def _nan_statistics(array):
+    nnans = np.count_nonzero(np.isnan(array))
+    print(f'N elements: {array.size}, N NaNs: {nnans}, NaN %: {100*nnans/array.size:.1f}')
+
+
 def calculate_optimal_grid_parameters(pnt_map_dict, antenna_name, telescope_diameter, chan_freq, ddi):
     reference_frequency = np.median(chan_freq)
     reference_lambda = scipy.constants.speed_of_light / reference_frequency
@@ -322,18 +327,21 @@ def calculate_optimal_grid_parameters(pnt_map_dict, antenna_name, telescope_diam
     data_range = \
         (pnt_map_dict[antenna_name].POINTING_OFFSET.values[:, 1].max()
          - pnt_map_dict[antenna_name].POINTING_OFFSET.values[:, 1].min())
+    pnt_off = pnt_map_dict[antenna_name].POINTING_OFFSET.values
 
     logger.info(f'{create_dataset_label(antenna_name, ddi)}: Cell size {format_angular_distance(cell_size)}, '
                 f'FOV: {format_angular_distance(data_range)}')
     # logger.info(f"cell_size: {cell_size}")
     # logger.info(f"data_range: {data_range}")
 
+    n_pix = 1
     try:
         n_pix = int(np.ceil(data_range / cell_size)) ** 2
 
     except ZeroDivisionError:
-        logger.error(f"Zero division error, there was likely a problem calculating the data range.", verbose=True)
-        raise ZeroDivisionError
+        logger.warning(f"{antenna_name} has a zero sized sky cell", verbose=True)
+    except ValueError:
+        logger.warning(f"{antenna_name} has a NaN valued data range", verbose=True)
 
     return n_pix, cell_size