Merge pull request #402 from HERA-Team/enu_from_ecef_api_fix

Fix api call to ENU_from_ECEF for pyuvdata 3.0

.github/workflows/ci.yml

@@ -16,7 +16,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:

examples/UVWindow_tutorial.ipynb

@@ -119,13 +119,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from astropy.io import fits\n",
     "from astropy import units\n",
-    "from scipy.interpolate import interp2d, interp1d\n",
-    "from scipy import integrate\n",
-    "import sys, matplotlib, os, h5py\n",
     "import numpy as np\n",
-    "from itertools import product, combinations\n",
     "from astropy import units\n",
     "from multiprocessing import Pool\n",
     "import time, h5py\n",

hera_pspec/pspecbeam.py

@@ -79,7 +79,7 @@ def _compute_pspec_scalar(cosmo, beam_freqs, omega_ratio, pspec_freqs,
     X2Y = np.array([cosmo.X2Y(z, little_h=little_h) for z in redshifts])
 
     if exact_norm: #Beam and spectral tapering are already taken into account in normalization. We only use averaged X2Y
-        scalar = integrate.trapz(X2Y, x=integration_freqs)/(np.abs(integration_freqs[-1]-integration_freqs[0]))
+        scalar = integrate.trapezoid(X2Y, x=integration_freqs)/(np.abs(integration_freqs[-1]-integration_freqs[0]))
         return scalar
 
     # Use linear interpolation to interpolate the frequency-dependent
@@ -105,7 +105,7 @@ def _compute_pspec_scalar(cosmo, beam_freqs, omega_ratio, pspec_freqs,
 
     # Integrate to get scalar
     d_inv_scalar = dBpp_over_BpSq * dOpp_over_Op2 / X2Y
-    scalar = 1. / integrate.trapz(d_inv_scalar, x=integration_freqs)
+    scalar = 1. / integrate.trapezoid(d_inv_scalar, x=integration_freqs)
     return scalar
 
 

hera_pspec/pspecdata.py

@@ -198,9 +198,6 @@ def add(self, dsets, wgts, labels=None, dsets_std=None, cals=None, cal_flag=True
         for d, w, s in zip(dsets, wgts, dsets_std):
             if not isinstance(d, UVData):
                 raise TypeError("Only UVData objects can be used as datasets.")
-            elif not d.future_array_shapes:
-                warnings.warn('Converting data to future_array_shapes...')
-                d.use_future_array_shapes()
             if not isinstance(w, UVData) and w is not None:
                 raise TypeError("Only UVData objects (or None) can be used as "
                                 "weights.")
@@ -322,7 +319,13 @@ def validate_datasets(self, verbose=True):
 
         # Check phase type
         phase_types = []
-        for d in self.dsets: phase_types.append(d.phase_type)
+        for d in self.dsets: 
+            if len(d.phase_center_catalog) > 1:
+                raise ValueError(
+                    "phase_center_catalog should contain only one entry per dataset"
+                )
+
+            phase_types.append(next(iter(d.phase_center_catalog.values()))['cat_type'])
         if np.unique(phase_types).size > 1:
             raise ValueError("all datasets must have the same phase type "
                              "(i.e. 'drift', 'phased', ...)\ncurrent phase "
@@ -362,7 +365,7 @@ def check_key_in_dset(self, key, dset_ind):
         """
         #FIXME: Fix this to enable label keys
         # get iterable
-        key = uvutils._get_iterable(key)
+        key = uvutils.tools._get_iterable(key)
         if isinstance(key, str):
             key = (key,)
 
@@ -1721,7 +1724,7 @@ def get_unnormed_V(self, key1, key2, model='empirical', exact_norm=False,
 
     def get_analytic_covariance(self, key1, key2, M=None, exact_norm=False,
                                 pol=False, model='empirical', known_cov=None):
-        """
+        r"""
         Calculates the auto-covariance matrix for both the real and imaginary
         parts of bandpowers (i.e., the q vectors and the p vectors).
 
@@ -2242,7 +2245,7 @@ def get_Q_alt(self, mode: int, allow_fft=True, include_extension=False):
         )
 
     def get_Q_alt_tensor(self, allow_fft=True, include_extension=False):
-        """
+        r"""
         Gets the (Ndly, Nfreq, Nfreq) tensor of Q_alt matrices.
 
         Parameters
@@ -2349,7 +2352,7 @@ def get_Q(self, mode):
         return Q_alt
 
     def p_hat(self, M, q):
-        """
+        r"""
         Optimal estimate of bandpower p_alpha, defined as p_hat = M q_hat.
 
         Parameters
@@ -2526,7 +2529,7 @@ def delays(self):
 
     def scalar(self, polpair, little_h=True, num_steps=2000, beam=None,
                taper_override='no_override', exact_norm=False):
-        """
+        r"""
         Computes the scalar function to convert a power spectrum estimate
         in "telescope units" to cosmological units, using self.spw_range to set
         spectral window.
@@ -2610,7 +2613,7 @@ def scalar(self, polpair, little_h=True, num_steps=2000, beam=None,
 
     def scalar_delay_adjustment(self, key1=None, key2=None, sampling=False,
                                 Gv=None, Hv=None):
-        """
+        r"""
         Computes an adjustment factor for the pspec scalar. There are
         two reasons why this might be needed:
 
@@ -3426,15 +3429,17 @@ def pspec(self, bls1, bls2, dsets, pols, n_dlys=None,
                     # insert time and blpair info only once per blpair
                     if i < 1 and j < 1:
                         # insert time info
+                        # inds is a slice if dset1 has rectangular blts.
                         inds1 = dset1.antpair2ind(bl1, ordered=False)
                         inds2 = dset2.antpair2ind(bl2, ordered=False)
+
                         time1.extend(dset1.time_array[inds1])
                         time2.extend(dset2.time_array[inds2])
                         lst1.extend(dset1.lst_array[inds1])
                         lst2.extend(dset2.lst_array[inds2])
 
                         # insert blpair info
-                        blp_arr.extend(np.ones_like(inds1, int) \
+                        blp_arr.extend(np.ones_like(dset1.time_array[inds1], int) \
                                        * uvputils._antnums_to_blpair(blp))
 
                 # insert into data and wgts integrations dictionaries
@@ -3644,16 +3649,19 @@ def rephase_to_dset(self, dset_index=0, inplace=True):
 
         # iterate over dsets
         for i, dset in enumerate(dsets):
+            if len(dset.phase_center_catalog) > 1:
+                raise ValueError("Cannot deal with datasets with more than one phase center")
+
             # don't rephase dataset we are using as our LST anchor
             if i == dset_index:
                 # even though not phasing this dset, must set to match all other
                 # dsets due to phasing-check validation
-                dset.phase_type = 'unknown'
+                dset.phase_center_catalog[0]['cat_type'] = 'unknown'
                 continue
 
             # skip if dataset is not drift phased
-            if dset.phase_type != 'drift':
-                warnings.warn(f"Skipping dataset {i} b/c it isn't drift phased")
+            if dset.phase_center_catalog[0]['cat_type'] != 'unprojected':
+                warnings.warn(f"Skipping dataset {i} b/c it isn't unprojected")
 
             # convert UVData to DataContainers. Note this doesn't make
             # a copy of the data
@@ -3685,7 +3693,7 @@ def rephase_to_dset(self, dset_index=0, inplace=True):
 
             # set phasing in UVData object to unknown b/c there isn't a single
             # consistent phasing for the entire data set.
-            dset.phase_type = 'unknown'
+            dset.phase_center_catalog[0]['cat_type'] = 'unknown'
 
         if inplace is False:
             return dsets
@@ -4579,9 +4587,9 @@ def _load_cals(cnames, logf=None, verbose=True):
         # read data
         uvc = UVCal()
         if isinstance(cfile, str):
-            uvc.read_calfits(glob.glob(cfile), use_future_array_shapes=True)
+            uvc.read(glob.glob(cfile))
         else:
-            uvc.read_calfits(cfile, use_future_array_shapes=True)
+            uvc.read(cfile)
         uvc.extra_keywords['filename'] = json.dumps(cfile)
         cals.append(uvc)
 

hera_pspec/testing.py

@@ -447,8 +447,8 @@ def noise_sim(
     Nextend = int(Nextend)
     if Nextend > 0:
         assert (
-            data.phase_type == "drift"
-        ), "data must be drift phased in order to extend along time axis"
+            data.phase_center_catalog[0]['cat_type'] == "unprojected"
+        ), "data must be unprojected in order to extend along time axis"
         data = copy.deepcopy(data)
         _data = copy.deepcopy(data)
         dt = np.median(np.diff(np.unique(_data.time_array)))

hera_pspec/tests/Scalar_dev2.ipynb

@@ -232,7 +232,7 @@
     }
    ],
    "source": [
-    "scalar = 1 / integrate.trapz(d_inv_scalar, pspec_freqs_Hz)\n",
+    "scalar = 1 / integrate.trapezoid(d_inv_scalar, pspec_freqs_Hz)\n",
     "print scalar#* (pspec_freqs_Hz[-1] - pspec_freqs_Hz[0])"
    ]
   },

hera_pspec/tests/test_pspecdata.py

@@ -3,7 +3,7 @@
 import numpy as np
 import pyuvdata as uv
 import os, copy, sys
-from scipy.integrate import simps, trapz
+from scipy.integrate import simpson, trapezoid
 from .. import pspecdata, pspecbeam, conversions, container, utils, testing, uvwindow
 from hera_pspec.data import DATA_PATH
 from pyuvdata import UVData, UVCal, utils as uvutils
@@ -192,25 +192,25 @@ def test_init(self):
         pytest.raises(ValueError, ds.get_G, key, key)
         pytest.raises(ValueError, ds.get_H, key, key)
 
-        # Test conversion if not future_array_shapes
-        uvd_old = uv.UVData()
-        uvd_old.read_miriad(os.path.join(DATA_PATH,
-                                          "zen.2458042.17772.xx.HH.uvXA"),
-                             use_future_array_shapes=False)
-        ds = pspecdata.PSpecData()
-        pytest.warns(UserWarning, ds.add, [uvd_old, uvd_old], [None, None])
 
     def test_add_data(self):
         """
         Test PSpecData add()
         """
         uv = self.d[0]
+        # proper usage
+        ds1 = copy.deepcopy(self.ds)
+        ds1.add(dsets=[uv], wgts=None, labels=None)
         # test adding non list objects
         pytest.raises(TypeError, self.ds.add, 1, 1)
         # test adding non UVData objects
         pytest.raises(TypeError, self.ds.add, [1], None)
         pytest.raises(TypeError, self.ds.add, [uv], [1])
         pytest.raises(TypeError, self.ds.add, [uv], None, dsets_std=[1])
+        # test adding UVData object with old array shape
+        ds2 = copy.deepcopy(uv)
+        ds2.data_array = np.tile(uv.data_array, [1,1,1,1])
+        pytest.raises(TypeError, self.ds.add, [uv], [1])
         # test adding non UVCal for cals
         pytest.raises(TypeError, self.ds.add, [uv], None, cals=[1])
         # test TypeError if dsets is dict but other inputs are not
@@ -1077,9 +1077,9 @@ def test_parseval(self):
         ps_avg = np.fft.fftshift( np.mean(ps[0], axis=1) )
 
         # Calculate integrals for Parseval's theorem
-        parseval_real = simps(gsq, x)
-        parseval_ft = dx**2. * simps(fsq, k)
-        parseval_phat = simps(ps_avg, tau)
+        parseval_real = simpson(gsq, x)
+        parseval_ft = dx**2. * simpson(fsq, k)
+        parseval_phat = simpson(ps_avg, tau)
 
         # Report on results for different ways of calculating Parseval integrals
         print "Parseval's theorem:"
@@ -1186,6 +1186,14 @@ def test_validate_datasets(self):
         pytest.raises(ValueError, ds.validate_datasets)
         uvd2.phase_to_time(Time(2458042.5, format='jd'))
         ds.validate_datasets()
+        # phase_center_catalog should contain only one entry per dataset
+        uvd3 = copy.deepcopy(self.d[0])
+        uvd3.phase_center_catalog[1] = uvd3.phase_center_catalog[0]
+        ds2 = pspecdata.PSpecData(dsets=[uvd, uvd3], wgts=[None, None])
+        pytest.raises(ValueError, ds2.validate_datasets)
+        # phased data
+        uvd4 = copy.deepcopy(self.d[0])
+
 
         # test polarization
         ds.validate_pol((0,1), ('xx', 'xx'))
@@ -1209,6 +1217,7 @@ def test_rephase_to_dset(self):
         uvp1 = ds.pspec(bls, bls, (0, 1), pols=('xx','xx'), verbose=False)
         # rephase and get pspec
         ds.rephase_to_dset(0)
+        ds2 = ds.rephase_to_dset(0, inplace=False)
         uvp2 = ds.pspec(bls, bls, (0, 1), pols=('xx','xx'), verbose=False)
         blp = (0, ((37,39),(37,39)), ('xx','xx'))
         assert np.isclose(np.abs(uvp2.get_data(blp)/uvp1.get_data(blp)), 1.0).min()
@@ -1375,7 +1384,7 @@ def test_get_analytic_covariance(self):
         for i in range(2):
             new = copy.deepcopy(uvd)
             new.time_array += new.Ntimes * np.diff(np.unique(new.time_array))[0]
-            new.lst_array = uvutils.get_lst_for_time(new.time_array, *new.telescope_location_lat_lon_alt_degrees)
+            new.lst_array = uvutils.get_lst_for_time(new.time_array, telescope_loc=new.telescope.location)
             uvd += new
 
         # get redundant baselines
@@ -1774,7 +1783,7 @@ def test_normalization(self):
 
         # taper
         window = windows.blackmanharris(len(freqs))
-        NEB = Bp / trapz(window**2, x=freqs)
+        NEB = Bp / trapezoid(window**2, x=freqs)
         scalar = cosmo.X2Y(np.mean(cosmo.f2z(freqs))) * np.mean(OmegaP**2/OmegaPP) * Bp * NEB
         data1 = d1.get_data(bls1[0])
         data2 = d2.get_data(bls2[0])
@@ -1812,13 +1821,13 @@ def test_broadcast_dset_flags(self):
 
         # test single integration being flagged within spw
         ds = pspecdata.PSpecData(dsets=[copy.deepcopy(uvd), copy.deepcopy(uvd)], wgts=[None, None])
-        ds.dsets[0].flag_array[ds.dsets[0].antpair2ind(24, 25, ordered=False)[3], 600, 0] = True
+        ds.dsets[0].flag_array[ds.dsets[0].antpair2ind(24, 25, ordered=False), 600, 0][3] = True
         ds.broadcast_dset_flags(spw_ranges=[(400, 800)], time_thresh=0.25, unflag=False)
         assert ds.dsets[0].get_flags(24, 25)[3, 400:800].all()
         assert ds.dsets[0].get_flags(24, 25)[3, :].all() == False
 
         # test pspec run sets flagged integration to have zero weight
-        uvd.flag_array[uvd.antpair2ind(24, 25, ordered=False)[3], 400, :] = True
+        uvd.flag_array[uvd.antpair2ind(24, 25, ordered=False), 400, :][3] = True
         ds = pspecdata.PSpecData(dsets=[copy.deepcopy(uvd), copy.deepcopy(uvd)], wgts=[None, None])
         ds.broadcast_dset_flags(spw_ranges=[(400, 450)], time_thresh=0.25)
         uvp = ds.pspec([(24, 25), (37, 38), (38, 39)], [(24, 25), (37, 38), (38, 39)], (0, 1), ('xx', 'xx'),
@@ -1830,7 +1839,7 @@ def test_broadcast_dset_flags(self):
         # average spectra
         avg_uvp = uvp.average_spectra(blpair_groups=[sorted(np.unique(uvp.blpair_array))], time_avg=True, inplace=False)
         # repeat but change data in flagged portion
-        ds.dsets[0].data_array[uvd.antpair2ind(24, 25, ordered=False)[3], 400:450, :] *= 100
+        ds.dsets[0].data_array[uvd.antpair2ind(24, 25, ordered=False), 400:450, :][3] *= 100
         uvp2 = ds.pspec([(24, 25), (37, 38), (38, 39)], [(24, 25), (37, 38), (38, 39)], (0, 1), ('xx', 'xx'),
                         spw_ranges=[(400, 450)], verbose=False)
         avg_uvp2 = uvp.average_spectra(blpair_groups=[sorted(np.unique(uvp.blpair_array))], time_avg=True, inplace=False)
@@ -2004,8 +2013,8 @@ def test_pspec_run():
     uvd.read_miriad(fnames[0], use_future_array_shapes=True)
     # interleave the data by hand, and add some flags in
     uvd.flag_array[:] = False
-    uvd.flag_array[uvd.antpair2ind(37, 38, ordered=False)[0], 10, 0] = True
-    uvd.flag_array[uvd.antpair2ind(37, 38, ordered=False)[:3], 15, 0] = True
+    uvd.flag_array[uvd.antpair2ind(37, 38, ordered=False), 10, 0][0] = True
+    uvd.flag_array[uvd.antpair2ind(37, 38, ordered=False), 15, 0][:3] = True
     uvd1 = uvd.select(times=np.unique(uvd.time_array)[::2], inplace=False)
     uvd2 = uvd.select(times=np.unique(uvd.time_array)[1::2], inplace=False)
     if os.path.exists("./out2.h5"):