Merge pull request #959 from HERA-Team/usenumpy2

fix: allow to use numpy 2 and pyuvdata 3

.github/workflows/ci.yml

@@ -14,7 +14,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, macos-latest]
-        python-version: [3.9, "3.10", "3.11", "3.12"]
+        python-version: ["3.10", "3.11", "3.12"]
       fail-fast: false
 
     steps:
@@ -53,7 +53,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest]
-        python-version: ["3.11"]
+        python-version: ["3.12"]
       fail-fast: false
 
     steps:

hera_cal/abscal.py

@@ -363,6 +363,8 @@ def TT_phs_logcal(model, data, antpos, wgts=None, refant=None, assume_2D=True,
     # angle after divide
     with warnings.catch_warnings():
         warnings.filterwarnings("ignore", message="invalid value encountered in divide")
+        warnings.filterwarnings("ignore", message="divide by zero encountered in divide")
+
         ydata = {k: np.angle(data[k] / model[k]) for k in keys}
 
     # make unit weights if None
@@ -476,6 +478,7 @@ def amp_logcal(model, data, wgts=None, verbose=True):
     # difference of log-amplitudes is ydata independent variable
     with warnings.catch_warnings():
         warnings.filterwarnings("ignore", message="invalid value encountered in divide")
+        warnings.filterwarnings("ignore", message="divide by zero encountered in divide")
         ydata = odict([(k, np.log(np.abs(data[k] / model[k]))) for k in keys])
 
     # make weights if None
@@ -547,6 +550,7 @@ def phs_logcal(model, data, wgts=None, refant=None, verbose=True):
     # angle of visibility ratio is ydata independent variable
     with warnings.catch_warnings():
         warnings.filterwarnings("ignore", message="invalid value encountered in divide")
+        warnings.filterwarnings("ignore", message="divide by zero encountered in divide")
         ydata = odict([(k, np.angle(data[k] / model[k])) for k in keys])
 
     # make weights if None
@@ -2623,7 +2627,7 @@ def __init__(self, model, data, refant=None, wgts=None, antpos=None, freqs=None,
                  #!/usr/bin/env python
                  uvd = pyuvdata.UVData()
                  uvd.read_miriad(<filename>)
-                 antenna_pos, ants = uvd.get_ENU_antpos()
+                 antenna_pos, ants = uvd.telescope.get_enu_antpos()
                  antpos = dict(zip(ants, antenna_pos))
                  ----
                  This is needed only for Tip Tilt, phase slope, and delay slope calibration.
@@ -3554,8 +3558,11 @@ def get_all_times_and_lsts(hd, solar_horizon=90.0, unwrap=True):
 
     # remove times when sun was too high
     if solar_horizon < 90.0:
-        lat, lon, alt = hd.telescope_location_lat_lon_alt_degrees
-        solar_alts = utils.get_sun_alt(all_times, latitude=lat, longitude=lon)
+        solar_alts = utils.get_sun_alt(
+            all_times,
+            latitude=hd.telescope.location.lat.deg,
+            longitude=hd.telescope.location.lon.deg
+        )
         solar_flagged = solar_alts > solar_horizon
         return all_times[~solar_flagged], all_lsts[~solar_flagged]
     else:  # skip this step for speed
@@ -4051,8 +4058,8 @@ def post_redcal_abscal_run(data_file, redcal_file, model_files, raw_auto_file=No
             warnings.warn(f"Warning: Overwriting redcal gain_scale of {hc.gain_scale} with model gain_scale of {hdm.vis_units}", RuntimeWarning)
         hc.gain_scale = hdm.vis_units  # set vis_units of hera_cal based on model files.
         hd_autos = io.HERAData(raw_auto_file)
-        assert hdm.x_orientation.lower() == hd.x_orientation.lower(), 'Data x_orientation, {}, does not match model x_orientation, {}'.format(hd.x_orientation.lower(), hdm.x_orientation.lower())
-        assert hc.x_orientation.lower() == hd.x_orientation.lower(), 'Data x_orientation, {}, does not match redcal x_orientation, {}'.format(hd.x_orientation.lower(), hc.x_orientation.lower())
+        assert hdm.telescope.x_orientation.lower() == hd.telescope.x_orientation.lower(), 'Data x_orientation, {}, does not match model x_orientation, {}'.format(hd.telescope.x_orientation.lower(), hdm.telescope.x_orientation.lower())
+        assert hc.telescope.x_orientation.lower() == hd.telescope.x_orientation.lower(), 'Data x_orientation, {}, does not match redcal x_orientation, {}'.format(hd.telescope.x_orientation.lower(), hc.telescope.x_orientation.lower())
         pol_load_list = [pol for pol in hd.pols if split_pol(pol)[0] == split_pol(pol)[1]]
 
         # get model bls and antpos to use later in baseline matching
@@ -4110,7 +4117,7 @@ def post_redcal_abscal_run(data_file, redcal_file, model_files, raw_auto_file=No
                                 model_flags.select_or_expand_times(model_times_to_load)
                             model_blvecs = {bl: model.antpos[bl[0]] - model.antpos[bl[1]] for bl in model.keys()}
                             utils.lst_rephase(model, model_blvecs, model.freqs, data.lsts - model.lsts,
-                                              lat=hdm.telescope_location_lat_lon_alt_degrees[0], inplace=True)
+                                              lat=hdm.telescope.location.lat.deg, inplace=True)
 
                             # Flag frequencies and times in the data that are entirely flagged in the model
                             model_flag_waterfall = np.all([f for f in model_flags.values()], axis=0)
@@ -4346,19 +4353,21 @@ def run_model_based_calibration(data_file, model_file, output_filename, auto_fil
         refant_init = str(refant)
     # initialize HERACal to store new cal solutions
     hc = UVCal()
-    hc = hc.initialize_from_uvdata(uvdata=hdd, gain_convention='divide', cal_style='sky',
-                                   ref_antenna_name=refant_init, sky_catalog=f'{model_file}',
-                                   metadata_only=False, cal_type='gain',
-                                   future_array_shapes=True)
+    hc = hc.initialize_from_uvdata(
+        uvdata=hdd, gain_convention='divide', cal_style='sky',
+        ref_antenna_name=refant_init, sky_catalog=f'{model_file}',
+        metadata_only=False, cal_type='gain',
+    )
 
     hc = io.to_HERACal(hc)
     hc.update(flags=data_ant_flags)
     # generate cal object from model to hold model flags.
     hcm = UVCal()
-    hcm = hcm.initialize_from_uvdata(uvdata=hdm, gain_convention='divide', cal_style='sky',
-                                     ref_antenna_name=refant_init, sky_catalog=f'{model_file}',
-                                     metadata_only=False, cal_type='gain',
-                                     future_array_shapes=True)
+    hcm = hcm.initialize_from_uvdata(
+        uvdata=hdm, gain_convention='divide', cal_style='sky',
+        ref_antenna_name=refant_init, sky_catalog=f'{model_file}',
+        metadata_only=False, cal_type='gain',
+    )
     hcm = io.to_HERACal(hcm)
     hcm.update(flags=model_ant_flags)
     # init all gains to unity.

hera_cal/chunker.py

@@ -81,7 +81,7 @@ def chunk_files(filenames, inputfile, outputfile, chunk_size, type="data",  # no
         if spw_range is None:
             spw_range = (0, chunked_files.Nfreqs)
         # convert polarizations to jones integers.
-        jones = [uvutils.polstr2num(pol, x_orientation=chunked_files.x_orientation) for pol in polarizations]
+        jones = [uvutils.polstr2num(pol, x_orientation=chunked_files.telescope.x_orientation) for pol in polarizations]
         chunked_files.select(freq_chans=np.arange(spw_range[0], spw_range[1]).astype(int), jones=jones)
     # throw away fully flagged baselines.
     if throw_away_flagged_ants:

hera_cal/flag_utils.py

@@ -63,7 +63,7 @@ def solar_flag(flags, times=None, flag_alt=0.0, longitude=21.42830, latitude=-30
     elif isinstance(flags, UVData):
         if verbose:
             print("Note: using latitude and longitude in given UVData object")
-        latitude, longitude, altitude = flags.telescope_location_lat_lon_alt_degrees
+        latitude, longitude, altitude = flags.telescope._location.lat_lon_alt_degrees()
         times = np.unique(flags.time_array)
         dtype = 'uvd'
     if dtype in ['ndarr', 'DC']:

hera_cal/frf.py

@@ -285,8 +285,9 @@ def sky_frates(uvd, keys=None, frate_standoff=0.0, frate_width_multiplier=1.0,
     """
     if keys is None:
         keys = uvd.get_antpairpols()
-    antpos, antnums = uvd.get_ENU_antpos()
-    sinlat = np.sin(np.abs(uvd.telescope_location_lat_lon_alt[0]))
+    antpos = uvd.telescope.get_enu_antpos()
+    antnums = uvd.telescope.antenna_numbers
+    sinlat = np.sin(np.abs(uvd.telescope.location.lat))
     frate_centers = {}
     frate_half_widths = {}
 
@@ -412,11 +413,10 @@ def build_fringe_rate_profiles(uvd, uvb, keys=None, normed=True, combine_pols=Tr
     if keys is None:
         keys = uvd.get_antpairpols()
 
-    antpos_trf = uvd.antenna_positions  # earth centered antenna positions
-    antnums = uvd.antenna_numbers  # antenna numbers.
+    antpos_trf = uvd.telescope.antenna_positions  # earth centered antenna positions
+    antnums = uvd.telescope.antenna_numbers  # antenna numbers.
 
-    lat, lon, alt = uvd.telescope_location_lat_lon_alt_degrees
-    location = EarthLocation(lon=lon * units.deg, lat=lat * units.deg, height=alt * units.m)
+    location = uvd.telescope.location
 
     # get topocentricl AzEl Beam coordinates.
     hp = HEALPix(nside=uvb.nside, order=uvb.ordering)
@@ -462,7 +462,7 @@ def build_fringe_rate_profiles(uvd, uvb, keys=None, normed=True, combine_pols=Tr
     # for if we are going to sum over polarizations.
     # get redundancies (will only compute fr-profile once for each red group).
     if reds is None:
-        reds = _get_key_reds(dict(zip(*uvd.get_ENU_antpos()[::-1])), keys)
+        reds = _get_key_reds(utils.get_ENU_antpos(uvd, asdict=True), keys)
 
     for redgrp in reds:
         # only explicitly calculate fr profile for the first vis in each redgroup.
@@ -610,7 +610,7 @@ def get_fringe_rate_limits(uvd, uvb=None, frate_profiles=None, percentile_low=5.
             dfr = 1. / (nfr * np.mean(np.diff(np.unique(uvd.time_array))) * 1e-3 * SDAY_SEC)
         fr_grid = np.arange(-nfr // 2, nfr // 2) * dfr
 
-    reds = _get_key_reds(dict(zip(*uvd.get_ENU_antpos()[::-1])), keys)
+    reds = _get_key_reds(utils.get_ENU_antpos(uvd, asdict=True), keys)
 
     for redgrp in reds:
         bl0 = redgrp[0]
@@ -1500,8 +1500,8 @@ def time_avg_data_and_write(input_data_list, output_data, t_avg, baseline_list=N
                               times=avg_times, lsts=avg_lsts, filetype=filetype, overwrite=clobber)
                 if flag_output is not None:
                     uv_avg = UVData()
-                    uv_avg.read(output_data_name, use_future_array_shapes=True)
-                    uvf = UVFlag(uv_avg, mode='flag', copy_flags=True, use_future_array_shapes=True)
+                    uv_avg.read(output_data_name)
+                    uvf = UVFlag(uv_avg, mode='flag', copy_flags=True)
                     uvf.to_waterfall(keep_pol=False, method='and')
                     uvf.write(os.path.join(os.path.dirname(flag_output),
                                            os.path.basename(flag_output).replace('.h5', f'.interleave_{inum}.h5')), clobber=clobber)
@@ -1512,8 +1512,8 @@ def time_avg_data_and_write(input_data_list, output_data, t_avg, baseline_list=N
                           nsamples=fr.avg_nsamples, times=fr.avg_times, lsts=fr.avg_lsts)
             if flag_output is not None:
                 uv_avg = UVData()
-                uv_avg.read(output_data, use_future_array_shapes=True)
-                uvf = UVFlag(uv_avg, mode='flag', copy_flags=True, use_future_array_shapes=True)
+                uv_avg.read(output_data)
+                uvf = UVFlag(uv_avg, mode='flag', copy_flags=True)
                 uvf.to_waterfall(keep_pol=False, method='and')
                 uvf.write(flag_output, clobber=clobber)
 
@@ -1817,7 +1817,7 @@ def load_tophat_frfilter_and_write(
         # Read in the beam file if provided.
         if beamfitsfile is not None:
             uvb = UVBeam()
-            uvb.read_beamfits(beamfitsfile, use_future_array_shapes=True)
+            uvb.read_beamfits(beamfitsfile)
         else:
             uvb = None