Repository: RadioAstronomySoftwareGroup/pyuvdata

Issue:
uv beam frequency interpolation with new object returns tuple
calling _interp_freq with new_object=true gives a tuple with a uv beam object as the zeroth entry and 'None' as the first entry. Is this the desired behavior? 

Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/RadioAstronomySoftwareGroup/pyuvdata/pull/631?src=pr&el=h1) Report
> Merging [#631](https://codecov.io/gh/RadioAstronomySoftwareGroup/pyuvdata/pull/631?src=pr&el=desc) into [master](https://codecov.io/gh/RadioAstronomySoftwareGroup/pyuvdata/commit/67eca3a770ec22814fac626fab59f6834f136eb8?src=pr&el=desc) will **decrease** coverage by `0.09%`.
> The diff coverage is `93.54%`.

[![Impacted file tree graph](https://codecov.io/gh/RadioAstronomySoftwareGroup/pyuvdata/pull/631/graphs/tree.svg?width=650&token=4QykcCeB1H&height=150&src=pr)](https://codecov.io/gh/RadioAstronomySoftwareGroup/pyuvdata/pull/631?src=pr&el=tree)

```diff
@@            Coverage Diff            @@
##           master     #631     +/-   ##
=========================================
- Coverage   98.87%   98.78%   -0.1%     
=========================================
  Files          20       20             
  Lines        8062     8120     +58     
=========================================
+ Hits         7971     8021     +50     
- Misses         91       99      +8
```


| [Impacted Files](https://codecov.io/gh/RadioAstronomySoftwareGroup/pyuvdata/pull/631?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [pyuvdata/uvbeam.py](https://codecov.io/gh/RadioAstronomySoftwareGroup/pyuvdata/pull/631/diff?src=pr&el=tree#diff-cHl1dmRhdGEvdXZiZWFtLnB5) | `98.63% <93.54%> (-0.65%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/RadioAstronomySoftwareGroup/pyuvdata/pull/631?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/RadioAstronomySoftwareGroup/pyuvdata/pull/631?src=pr&el=footer). Last update [67eca3a...36c82ec](https://codecov.io/gh/RadioAstronomySoftwareGroup/pyuvdata/pull/631?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

jpober: Sam is on vacation until next week.  We'll take a look then and see if this is a clean fix for our use case.  (PS: Why can't I tag her here?  She opened an issue, but maybe you can only tag people that are part of the team?)
mkolopanis: It's possible you might be able to tag but it won't autocomplete for you? let's try: @samanthamcgraw
bhazelton: Great. @aelanman can you check to make sure this doesn't break anything in healvis or pyuvsim?
aelanman: @bhazelton It does seem to have broken some pyuvsim tests, but this might just be issues that can be addressed in pyuvsim. One is because `UVBeam.interp` no longer accepts the keyword `freq_interp_kind`. I've fixed that by setting the attribute on each UVBeam object.

There's another that's a little harder to track down. One test is bringing up an error from scipy.interp1d:
`ValueError: The number of derivatives at boundaries does not match: expected 2, got 0+0`

This is coming from the frequency interpolation step.

`healvis` on the other hand needs to be updated to call the new functions, instead of the private `_interp_freq` method.
aelanman: Update -- This last error is caused by using cubic frequency interpolation instead of linear. Curious that this test passed before...

I realized that although the antenna class functions accept `freq_interp_kind` as a keyword, there is no way for the user to set interpolation methods on the UVBeams. I'll make an issue on pyuvsim for that.
bhazelton: @jpober @samanthamcgraw I think I fixed the bug you told me about on this branch. Sorry for the delay! Let me know if you still have trouble.

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index ba71c363..9d146b08 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,11 +4,15 @@ All notable changes to this project will be documented in this file.
 ## [Unreleased]
 
 ### Added
+- option for `UVBeam.interp` to return a new beam object.
 - `UVFlag` information on Read The Docs
 
 ### Changed
 - `UVFlag` inherits from `UVBase` object.
 
+### Fixed
+- A bug that gave the wrong error message when calling `UVData.phase_to_time` without an Astropy Time object.
+
 ## [1.4.1] - 2019-08-2
 
 ### Added
@@ -23,6 +27,9 @@ All notable changes to this project will be documented in this file.
 - Fixed a bug in select that caused bls and antenna_names/numbers to be or'ed rather than and'ed together.
 - Fixed a bug where `baseline_to_antnums` could accept a numpy array as input but not other array_like objects.
 
+### Changed
+- removed `new_object` keyword from `UVBeam._interp_freq` in favor of new functionality in `UVBeam.interp`
+
 ## [1.4.0] - 2019-05-23
 
 ### Added
diff --git a/docs/make_flag_parameters.py b/docs/make_flag_parameters.py
index d80c3d9f..742323a9 100644
--- a/docs/make_flag_parameters.py
+++ b/docs/make_flag_parameters.py
@@ -13,7 +13,7 @@ from astropy.time import Time
 def write_flagparams_rst(write_file=None):
     test_file = os.path.join(DATA_PATH, 'zen.2457698.40355.xx.HH.uvcAA.testuvflag.h5')
     UV = UVFlag(test_file)
-    out = 'UVFLag Parameters\n=================\n'
+    out = 'UVFlag Parameters\n=================\n'
     out += ("These are the standard attributes of UVFlag objects.\n\nUnder the hood "
             "they are actually properties based on UVParameter objects.\n\n")
     out += 'Required\n----------------\n'
diff --git a/pyuvdata/uvbeam.py b/pyuvdata/uvbeam.py
index 36b69263..63ad448d 100644
--- a/pyuvdata/uvbeam.py
+++ b/pyuvdata/uvbeam.py
@@ -737,40 +737,37 @@ class UVBeam(UVBase):
         if not inplace:
             return beam_object
 
-    def _interp_freq(self, freq_array, kind='linear', tol=1.0, new_object=False):
+    def _interp_freq(self, freq_array, kind='linear', tol=1.0):
         """
         Simple interpolation function for frequency axis.
 
-        Args:
-            freq_array: frequency values [Hz] to interpolate to
-            kind: str, interpolation method, see scipy.interpolate.interp1d
-            tol: float, frequency distance tolerance [Hz] of nearest neighbors.
-                If *all* elements in freq_array have nearest neighbor distances within
-                the specified tolerance then return the beam at each nearest neighbor,
-                otherwise interpolate the beam.
-            new_object: bool, if True return a new UVBeam object, else return just the interpolated data
-
-        Returns:
-            interpolated beam values, or UVBeam object if new_object==True
-                shape: (Naxes_vec, Nspws, Nfeeds or Npols, freq_array.size, Npixels or (Naxis2, Naxis1))
-            interpolated bandpass values or None if new_object==True
-                shape: (Nspws, freq_array.size)
+        Parameters
+        ----------
+        freq_array : array_like of floats
+            Frequency values to interpolate to.
+        kind : str
+            Interpolation method to use frequency.
+            See scipy.interpolate.interp1d for details.
+
+        Returns
+        --------
+        interp_data : array_like of float or complex
+            The array of interpolated data values,
+            shape: (Naxes_vec, Nspws, Nfeeds or Npols, freq_array.size, Npixels or (Naxis2, Naxis1))
+        interp_bandpass : array_like of float
+            The interpolated bandpass. shape: (Nspws, freq_array.size)
         """
         assert(isinstance(freq_array, np.ndarray))
         assert(freq_array.ndim == 1)
 
         nfreqs = freq_array.size
 
-        # get frequency distances
-        freq_dists = np.abs(self.freq_array - freq_array.reshape(-1, 1))
-        nearest_dist = np.min(freq_dists, axis=1)
-        nearest_inds = np.argmin(freq_dists, axis=1)
-        interp_bool = np.any(nearest_dist >= tol)
-
-        # use the beam at nearest neighbors if not interp_bool
-        if not interp_bool:
-            interp_arrays = [self.data_array[:, :, :, nearest_inds, :], self.bandpass_array[:, nearest_inds]]
-            kind = 'nearest'
+        # use the beam at nearest neighbors if not kind is 'nearest'
+        if kind == 'nearest':
+            freq_dists = np.abs(self.freq_array - freq_array.reshape(-1, 1))
+            nearest_inds = np.argmin(freq_dists, axis=1)
+            interp_arrays = [self.data_array[:, :, :, nearest_inds, :],
+                             self.bandpass_array[:, nearest_inds]]
 
         # otherwise interpolate the beam
         else:
@@ -801,62 +798,62 @@ class UVBeam(UVBase):
 
                 interp_arrays.append(lut(freq_array))
 
-        # return just the interpolated arrays
-        if not new_object:
-            return tuple(interp_arrays)
-
-        # return a new UVBeam object with interpolated data
-        else:
-            # make a new object
-            new_uvb = self.select(freq_chans=np.arange(np.min([self.Nfreqs, len(freq_array)])), inplace=False)
-            new_uvb.data_array = interp_arrays[0]
-            new_uvb.Nfreqs = new_uvb.data_array.shape[3]
-            new_uvb.freq_array = freq_array.reshape(1, -1)
-            new_uvb.bandpass_array = interp_arrays[1]
-            new_uvb.freq_interp_kind = kind
-            if hasattr(new_uvb, 'saved_interp_functions'):
-                delattr(new_uvb, 'saved_interp_functions')
-
-            new_uvb.check()
-
-            return new_uvb, None
+        return tuple(interp_arrays)
 
-    def _interp_az_za_rect_spline(self, az_array, za_array, freq_array, freq_interp_kind='linear',
-                                  freq_interp_tol=1.0, reuse_spline=False, polarizations=None, **kwargs):
+    def _interp_az_za_rect_spline(self, az_array, za_array, freq_array,
+                                  freq_interp_kind='linear', freq_interp_tol=1.0,
+                                  polarizations=None, reuse_spline=False):
         """
         Simple interpolation function for az_za coordinate system.
 
-        Args:
-            az_array: az values to interpolate to (same length as za_array)
-            za_array: za values to interpolate to (same length as az_array)
-            freq_array: frequency values to interpolate to
-            freq_interp_kind: str, interpolation method across frequency. See scipy.interpolate.interp1d for details.
-            freq_interp_tol: float, frequency distance tolerance [Hz] of nearest neighbors.
-                If *all* elements in freq_array have nearest neighbor distances within
-                the specified tolerance then return the beam at each nearest neighbor,
-                otherwise interpolate the beam.
-            reuse_spline: Save the interpolation functions for reuse.
-            polarizations: list of str, polarizations to interpolate if beam_type is 'power'.
-                Default is all polarizations in self.polarization_array.
+        Parameters
+        ----------
+        az_array : array_like of floats
+            Azimuth values to interpolate to in radians, specifying the azimuth
+            positions for every interpolation point (same length as `za_array`).
+        za_array : array_like of floats
+            Zenith values to interpolate to in radians, specifying the zenith
+            positions for every interpolation point (same length as `az_array`).
+        freq_array : array_like of floats
+            Frequency values to interpolate to.
+        freq_interp_kind : str
+            Interpolation method to use frequency.
+            See scipy.interpolate.interp1d for details.
+        polarizations : list of str
+            polarizations to interpolate if beam_type is 'power'.
+            Default is all polarizations in self.polarization_array.
+        reuse_spline : bool
+            Save the interpolation functions for reuse.
+
+        Returns
+        -------
+        interp_data : array_like of float or complex
+            The array of interpolated data values,
+            shape: (Naxes_vec, Nspws, Nfeeds or Npols, Nfreqs, az_array.size)
+        interp_basis_vector : array_like of float
+            The array of interpolated basis vectors,
+            shape: (Naxes_vec, Ncomponents_vec, az_array.size)
+        interp_bandpass : array_like of float
+            The interpolated bandpass. shape: (Nspws, freq_array.size)
 
-        Returns:
-            an array of interpolated values, shape: (Naxes_vec, Nspws, Nfeeds or Npols, Nfreqs, az_array.size)
-            an array of interpolated basis vectors, shape: (Naxes_vec, Ncomponents_vec, az_array.size)
         """
         if self.pixel_coordinate_system != 'az_za':
             raise ValueError('pixel_coordinate_system must be "az_za"')
 
         if freq_array is not None:
             assert(isinstance(freq_array, np.ndarray))
-            input_data_array, _ = self._interp_freq(freq_array, kind=freq_interp_kind, tol=freq_interp_tol)
+            input_data_array, interp_bandpass = self._interp_freq(freq_array,
+                                                                  kind=freq_interp_kind,
+                                                                  tol=freq_interp_tol)
             input_nfreqs = freq_array.size
         else:
             input_data_array = self.data_array
             input_nfreqs = self.Nfreqs
             freq_array = self.freq_array[0]
+            interp_bandpass = self.bandpass_array[0]
 
         if az_array is None or za_array is None:
-            return input_data_array, self.basis_vector_array
+            return input_data_array, self.basis_vector_array, interp_bandpass
 
         assert(isinstance(az_array, np.ndarray))
         assert(isinstance(za_array, np.ndarray))
@@ -887,7 +884,8 @@ class UVBeam(UVBase):
                                       phi_vals,
                                       phi_vals[:, -1 * extend_length:] + extend_length * axis1_diff),
                                      axis=1)
-            theta_use = np.concatenate((theta_vals[:, :extend_length], theta_vals, theta_vals[:, -1 * extend_length:]),
+            theta_use = np.concatenate((theta_vals[:, :extend_length], theta_vals,
+                                        theta_vals[:, -1 * extend_length:]),
                                        axis=1)
 
             low_slice = input_data_array[:, :, :, :, :, :extend_length]
@@ -935,7 +933,8 @@ class UVBeam(UVBase):
                 pol_inds = []
                 for pol in pols:
                     if pol not in self.polarization_array:
-                        raise ValueError("Requested polarization {} not found in self.polarization_array".format(pol))
+                        raise ValueError('Requested polarization {} not found '
+                                         'in self.polarization_array'.format(pol))
                     pol_inds.append(np.where(self.polarization_array == pol)[0][0])
                 pol_inds = np.asarray(pol_inds)
                 Npol_feeds = len(pol_inds)
@@ -966,17 +965,20 @@ class UVBeam(UVBase):
                         if do_interp:
                             if np.iscomplexobj(data_use):
                                 # interpolate real and imaginary parts separately
-                                real_lut = interpolate.RectBivariateSpline(theta_use[:, 0],
-                                                                           phi_use[0, :],
-                                                                           data_use[index0, index1, index2, index3, :].real)
-                                imag_lut = interpolate.RectBivariateSpline(theta_use[:, 0],
-                                                                           phi_use[0, :],
-                                                                           data_use[index0, index1, index2, index3, :].imag)
+                                real_lut = \
+                                    interpolate.RectBivariateSpline(theta_use[:, 0],
+                                                                    phi_use[0, :],
+                                                                    data_use[index0, index1, index2, index3, :].real)
+                                imag_lut = \
+                                    interpolate.RectBivariateSpline(theta_use[:, 0],
+                                                                    phi_use[0, :],
+                                                                    data_use[index0, index1, index2, index3, :].imag)
                                 lut = get_lambda(real_lut, imag_lut)
                             else:
-                                lut = interpolate.RectBivariateSpline(theta_use[:, 0],
-                                                                      phi_use[0, :],
-                                                                      data_use[index0, index1, index2, index3, :])
+                                lut = \
+                                    interpolate.RectBivariateSpline(theta_use[:, 0],
+                                                                    phi_use[0, :],
+                                                                    data_use[index0, index1, index2, index3, :])
                                 lut = get_lambda(lut)
                             if reuse_spline:
                                 self.saved_interp_functions[key] = lut
@@ -989,28 +991,43 @@ class UVBeam(UVBase):
                                                      'the UVBeam pixel coverage.')
                         interp_data[index0, index1, pol_return_ind, index3, :] = lut(za_array, az_array)
 
-        return interp_data, interp_basis_vector
+        return interp_data, interp_basis_vector, interp_bandpass
 
-    def _interp_healpix_bilinear(self, az_array, za_array, freq_array, freq_interp_kind='linear',
-                                 freq_interp_tol=1.0, polarizations=None, **kwargs):
+    def _interp_healpix_bilinear(self, az_array, za_array, freq_array,
+                                 freq_interp_kind='linear',
+                                 freq_interp_tol=1.0, polarizations=None,
+                                 reuse_spline=False):
         """
         Simple bi-linear interpolation wrapper for healpix.
 
-        Args:
-            az_array: azimuth angles to interpolate to [radians]
-            za_array: zenith angles to interpolate to [radians]
-            freq_array: frequency values to interpolate to [Hz]
-            freq_interp_kind: str, interpolation method across frequency. See scipy.interpolate.interp1d for details.
-            freq_interp_tol: float, frequency distance tolerance [Hz] of nearest neighbors.
-                If *all* elements in freq_array have nearest neighbor distances within
-                the specified tolerance then return the beam at each nearest neighbor,
-                otherwise interpolate the beam.
-            polarizations: list of str, polarizations to interpolate if beam_type is 'power'.
-                Default is all polarizations in self.polarization_array.
+        Parameters
+        ----------
+        az_array : array_like of floats
+            Azimuth values to interpolate to in radians, specifying the azimuth
+            positions for every interpolation point (same length as `za_array`).
+        za_array : array_like of floats
+            Zenith values to interpolate to in radians, specifying the zenith
+            positions for every interpolation point (same length as `az_array`).
+        freq_array : array_like of floats
+            Frequency values to interpolate to.
+        freq_interp_kind : str
+            Interpolation method to use frequency.
+            See scipy.interpolate.interp1d for details.
+        polarizations : list of str
+            polarizations to interpolate if beam_type is 'power'.
+            Default is all polarizations in self.polarization_array.
+
+        Returns
+        -------
+        interp_data : array_like of float or complex
+            The array of interpolated data values,
+            shape: (Naxes_vec, Nspws, Nfeeds or Npols, Nfreqs, az_array.size)
+        interp_basis_vector : array_like of float
+            The array of interpolated basis vectors,
+            shape: (Naxes_vec, Ncomponents_vec, az_array.size)
+        interp_bandpass : array_like of float
+            The interpolated bandpass. shape: (Nspws, freq_array.size)
 
-        Returns:
-            an array of interpolated values, shape: (Naxes_vec, Nspws, Nfeeds or Npols, Nfreqs, az_array.size)
-            an array of interpolated basis vectors, shape: (Naxes_vec, Ncomponents_vec, az_array.size)
         """
         try:
             import healpy as hp
@@ -1023,15 +1040,18 @@ class UVBeam(UVBase):
 
         if freq_array is not None:
             assert(isinstance(freq_array, np.ndarray))
-            input_data_array, _ = self._interp_freq(freq_array, kind=freq_interp_kind, tol=freq_interp_tol)
+            input_data_array, interp_bandpass = self._interp_freq(freq_array,
+                                                                  kind=freq_interp_kind,
+                                                                  tol=freq_interp_tol)
             input_nfreqs = freq_array.size
         else:
             input_data_array = self.data_array
             input_nfreqs = self.Nfreqs
             freq_array = self.freq_array[0]
+            interp_bandpass = self.bandpass_array[0]
 
         if az_array is None or za_array is None:
-            return input_data_array, self.basis_vector_array
+            return input_data_array, self.basis_vector_array, interp_bandpass
 
         assert(isinstance(az_array, np.ndarray))
         assert(isinstance(za_array, np.ndarray))
@@ -1047,11 +1067,13 @@ class UVBeam(UVBase):
                 Npol_feeds = self.Npols
                 pol_inds = np.arange(Npol_feeds)
             else:
-                pols = [uvutils.polstr2num(p, x_orientation=self.x_orientation) for p in polarizations]
+                pols = [uvutils.polstr2num(p, x_orientation=self.x_orientation)
+                        for p in polarizations]
                 pol_inds = []
                 for pol in pols:
                     if pol not in self.polarization_array:
-                        raise ValueError("Requested polarization {} not found in self.polarization_array".format(pol))
+                        raise ValueError('Requested polarization {} not found '
+                                         'in self.polarization_array'.format(pol))
                     pol_inds.append(np.where(self.polarization_array == pol)[0][0])
                 pol_inds = np.asarray(pol_inds)
                 Npol_feeds = len(pol_inds)
@@ -1063,7 +1085,8 @@ class UVBeam(UVBase):
             data_type = np.complex
         else:
             data_type = np.float
-        interp_data = np.zeros((self.Naxes_vec, self.Nspws, Npol_feeds, input_nfreqs, len(az_array)), dtype=data_type)
+        interp_data = np.zeros((self.Naxes_vec, self.Nspws, Npol_feeds,
+                                input_nfreqs, len(az_array)), dtype=data_type)
 
         if self.basis_vector_array is not None:
             if (np.any(self.basis_vector_array[0, 1, :] > 0)
@@ -1092,58 +1115,241 @@ class UVBeam(UVBase):
                     for index2 in range(Npol_feeds):
                         if np.iscomplexobj(input_data_array):
                             # interpolate real and imaginary parts separately
-                            real_hmap = hp.get_interp_val(input_data_array[index0, index1, pol_inds[index2], index3, :].real, za_array, az_array)
-                            imag_hmap = hp.get_interp_val(input_data_array[index0, index1, pol_inds[index2], index3, :].imag, za_array, az_array)
+                            real_hmap = hp.get_interp_val(
+                                input_data_array[index0, index1, pol_inds[index2], index3, :].real,
+                                za_array, az_array)
+                            imag_hmap = hp.get_interp_val(
+                                input_data_array[index0, index1, pol_inds[index2], index3, :].imag,
+                                za_array, az_array)
                             hmap = real_hmap + 1j * imag_hmap
                         else:
                             # interpolate once
-                            hmap = hp.get_interp_val(input_data_array[index0, index1, pol_inds[index2], index3, :], za_array, az_array)
+                            hmap = hp.get_interp_val(
+                                input_data_array[index0, index1, pol_inds[index2], index3, :],
+                                za_array, az_array)
 
                         interp_data[index0, index1, index2, index3, :] = hmap
 
-        return interp_data, interp_basis_vector
+        return interp_data, interp_basis_vector, interp_bandpass
 
-    def interp(self, az_array=None, za_array=None, freq_array=None, freq_interp_tol=1.0,
-               polarizations=None, freq_interp_kind='linear', reuse_spline=False):
+    def interp(self, az_array=None, za_array=None, az_za_grid=False,
+               healpix_nside=None, healpix_inds=None,
+               freq_array=None, freq_interp_tol=1.0, polarizations=None,
+               return_bandpass=False, reuse_spline=False, new_object=False,
+               run_check=True, check_extra=True, run_check_acceptability=True):
         """
-        Interpolate beam to given az, za locations (in radians).
-
-        Args:
-            az_array: az values to interpolate to (same length as za_array)
-            za_array: za values to interpolate to (same length as az_array)
-            freq_array: frequency values to interpolate to
-            freq_interp_tol: float, frequency distance tolerance [Hz] of nearest neighbors.
-                If *all* elements in freq_array have nearest neighbor distances within
-                the specified tolerance then return the beam at each nearest neighbor,
-                otherwise interpolate the beam.
-            polarizations: list of str, polarizations to interpolate if beam_type is 'power'.
-                Default is all polarizations in self.polarization_array.
-            freq_interp_kind: str, interpolation method across frequency. See scipy.interpolate.interp1d for details.
-            reuse_spline: Save the interpolation functions for reuse. Only applies for `az_za_simple` interpolation.
-
-        Returns:
-            an array of interpolated values, shape: (Naxes_vec, Nspws, Nfeeds or Npols,
-                Nfreqs or freq_array.size if freq_array is passed,
-                Npixels/(Naxis1, Naxis2) or az_array.size if az/za_arrays are passed)
+        Interpolate beam to given frequency, az & za locations or Healpix pixel centers.
+
+        Parameters
+        ----------
+        az_array : array_like of floats, optional
+            Azimuth values to interpolate to in radians, either specifying the
+            azimuth positions for every interpolation point or specifying the
+            azimuth vector for a meshgrid if az_za_grid is True.
+        za_array : array_like of floats, optional
+            Zenith values to interpolate to in radians, either specifying the
+            zenith positions for every interpolation point or specifying the
+            zenith vector for a meshgrid if az_za_grid is True.
+        az_za_grid : bool
+            Option to treat the `az_array` and `za_array` as the input vectors
+            for points on a mesh grid.
+        healpix_nside : int, optional
+            HEALPix nside parameter if interpolating to HEALPix pixels.
+        healpix_inds : array_like of int, optional
+            HEALPix indices to interpolate to. Defaults to all indices in the
+            map if `healpix_nside` is set and `az_array` and `za_array` are None.
+        freq_array : array_like of floats, optional
+            Frequency values to interpolate to.
+        freq_interp_tol : float
+            Frequency distance tolerance [Hz] of nearest neighbors.
+            If *all* elements in freq_array have nearest neighbor distances within
+            the specified tolerance then return the beam at each nearest neighbor,
+            otherwise interpolate the beam.
+        polarizations : list of str
+            polarizations to interpolate if beam_type is 'power'.
+            Default is all polarizations in self.polarization_array.
+        new_object : bool
+            Option to return a new UVData object with the interpolated data,
+            if possible. Note that this is only possible for Healpix pixels or
+            if az_za_grid is True and `az_array` and `za_array` are evenly spaced
+            or for frequency only interpolation.
+        reuse_spline : bool
+            Save the interpolation functions for reuse. Only applies for
+            `az_za_simple` interpolation.
+        run_check : bool
+            Only used if new_object is True. Option to check for the existence
+            and proper shapes of required parameters on the new object.
+        check_extra : bool
+            Only used if new_object is True. Option to check optional parameters
+            as well as required ones on the new object.
+        run_check_acceptability : bool
+            Only used if new_object is True. Option to check acceptable range
+            of the values of required parameters on the new object.
+
+        Returns
+        -------
+        array_like of float or complex or a UVBeam object
+            Either an array of interpolated values or a UVBeam object if
+            `new_object` is True. The shape of the interpolated data will be:
+            (Naxes_vec, Nspws, Nfeeds or Npols, Nfreqs or freq_array.size if
+            freq_array is passed, Npixels/(Naxis1, Naxis2) or az_array.size if
+            az/za_arrays are passed)
+        interp_basis_vector : array_like of float, optional
             an array of interpolated basis vectors (or self.basis_vector_array
-                if az/za_arrays are not passed), shape: (Naxes_vec, Ncomponents_vec,
-                Npixels/(Naxis1, Naxis2) or az_array.size if az/za_arrays are passed)
+            if az/za_arrays are not passed), shape: (Naxes_vec, Ncomponents_vec,
+            Npixels/(Naxis1, Naxis2) or az_array.size if az/za_arrays are passed)
+        interp_bandpass : array_like of float, optional
+            The interpolated bandpass, only returned if `return_bandpass` is True.
+            shape: (Nspws, freq_array.size)
+
         """
         if self.interpolation_function is None:
             raise ValueError('interpolation_function must be set on object first')
         if self.freq_interp_kind is None:
             raise ValueError('freq_interp_kind must be set on object first')
 
+        if new_object:
+            if not az_za_grid and az_array is not None:
+                raise ValueError('A new object can only be returned if '
+                                 'az_za_grid is True or for Healpix pixels or '
+                                 'for frequency only interpolation.')
+
+        kind_use = self.freq_interp_kind
+        if freq_array is not None:
+            # get frequency distances
+            freq_dists = np.abs(self.freq_array - freq_array.reshape(-1, 1))
+            nearest_dist = np.min(freq_dists, axis=1)
+            interp_bool = np.any(nearest_dist >= freq_interp_tol)
+
+            # use the beam at nearest neighbors if not interp_bool
+            if not interp_bool:
+                kind_use = 'nearest'
+
+        if az_za_grid:
+            az_array_use, za_array_use = np.meshgrid(az_array, za_array)
+            az_array_use = az_array_use.flatten()
+            za_array_use = za_array_use.flatten()
+        else:
+            az_array_use = copy.copy(az_array)
+            za_array_use = copy.copy(za_array)
+
+        if healpix_nside is not None or healpix_inds is not None:
+            if healpix_nside is None:
+                raise ValueError('healpix_nside must be set if healpix_inds is set.')
+            if az_array is not None or za_array is not None:
+                raise ValueError('healpix_nside and healpix_inds can not be '
+                                 'set if az_array or za_array is set.')
+            try:
+                import healpy as hp
+            except ImportError:  # pragma: no cover
+                uvutils._reraise_context('healpy is not installed but is required for '
+                                         'healpix functionality')
+            if healpix_inds is None:
+                healpix_inds = np.arange(hp.nside2npix(healpix_nside))
+            za_array_use, az_array_use = hp.pix2ang(healpix_nside, healpix_inds)
+
         interp_func = self.interpolation_function_dict[self.interpolation_function]['func']
-        return getattr(self, interp_func)(az_array, za_array, freq_array,
-                                          freq_interp_kind=self.freq_interp_kind,
-                                          freq_interp_tol=freq_interp_tol,
-                                          polarizations=polarizations,
-                                          reuse_spline=reuse_spline)
+
+        extra_keyword_dict = {}
+        if interp_func == '_interp_az_za_rect_spline':
+            extra_keyword_dict['reuse_spline'] = reuse_spline
+
+        interp_data, interp_basis_vector, interp_bandpass = \
+            getattr(self, interp_func)(az_array_use, za_array_use, freq_array,
+                                       freq_interp_kind=kind_use,
+                                       polarizations=polarizations,
+                                       **extra_keyword_dict)
+
+        # return just the interpolated arrays
+        if not new_object:
+            if return_bandpass:
+                return interp_data, interp_basis_vector, interp_bandpass
+            else:
+                return interp_data, interp_basis_vector
+
+        # return a new UVBeam object with interpolated data
+        else:
+            # make a new object
+            new_uvb = copy.deepcopy(self)
+
+            history_update_string = (' Interpolated')
+            if freq_array is not None:
+                history_update_string += (' in frequency')
+                new_uvb.Nfreqs = freq_array.size
+                new_uvb.freq_array = freq_array.reshape(1, -1)
+                new_uvb.bandpass_array = interp_bandpass
+                new_uvb.freq_interp_kind = kind_use
+
+            if az_array is not None:
+                if freq_array is not None:
+                    history_update_string += (' and')
+                if new_uvb.pixel_coordinate_system != 'az_za':
+                    input_desc = self.coordinate_system_dict[new_uvb.pixel_coordinate_system]['description']
+                    output_desc = self.coordinate_system_dict['az_za']['description']
+                    history_update_string += (' from ' + input_desc
+                                              + ' to ' + output_desc)
+                    new_uvb.pixel_coordinate_system = 'az_za'
+                    new_uvb.Npixels = None
+                    new_uvb.pixel_array = None
+                    new_uvb.nside = None
+                    new_uvb.ordering = None
+                else:
+                    history_update_string += (' to a new azimuth/zenith angle grid')
+
+                interp_data = interp_data.reshape(interp_data.shape[:-1]
+                                                  + (za_array.size, az_array.size))
+                if interp_basis_vector is not None:
+                    interp_basis_vector = interp_basis_vector.reshape(
+                        interp_basis_vector.shape[:-1] + (za_array.size, az_array.size))
+
+                new_uvb.axis1_array = az_array
+                new_uvb.axis2_array = za_array
+                new_uvb.Naxes1 = new_uvb.axis1_array.size
+                new_uvb.Naxes2 = new_uvb.axis2_array.size
+
+            elif healpix_nside is not None:
+                if freq_array is not None:
+                    history_update_string += (' and')
+                if new_uvb.pixel_coordinate_system != 'healpix':
+                    input_desc = self.coordinate_system_dict[new_uvb.pixel_coordinate_system]['description']
+                    output_desc = self.coordinate_system_dict['healpix']['description']
+                    history_update_string += (' from ' + input_desc
+                                              + ' to ' + output_desc)
+                    new_uvb.pixel_coordinate_system = 'healpix'
+                    new_uvb.Naxes1 = None
+                    new_uvb.axis1_array = None
+                    new_uvb.Naxes2 = None
+                    new_uvb.axis2_array = None
+                else:
+                    history_update_string += (' to a new healpix grid')
+
+                new_uvb.pixel_array = healpix_inds
+                new_uvb.Npixels = new_uvb.pixel_array.size
+                new_uvb.nside = healpix_nside
+                new_uvb.ordering = 'ring'
+
+            history_update_string += (' using pyuvdata with interpolation_function = '
+                                      + new_uvb.interpolation_function)
+            if freq_array is not None:
+                history_update_string += (' and freq_interp_kind = '
+                                          + new_uvb.freq_interp_kind)
+            history_update_string += '.'
+            new_uvb.history = new_uvb.history + history_update_string
+            new_uvb.data_array = interp_data
+            if new_uvb.basis_vector_array is not None:
+                new_uvb.basis_vector_array = interp_basis_vector
+
+            if hasattr(new_uvb, 'saved_interp_functions'):
+                delattr(new_uvb, 'saved_interp_functions')
+
+            new_uvb.set_cs_params()
+            if run_check:
+                new_uvb.check(check_extra=check_extra,
+                              run_check_acceptability=run_check_acceptability)
+            return new_uvb
 
     def to_healpix(self, nside=None, run_check=True, check_extra=True,
-                   run_check_acceptability=True,
-                   inplace=True):
+                   run_check_acceptability=True, inplace=True):
         """
         Convert beam in to healpix coordinates.
         The interpolation is done using the interpolation method specified in
@@ -1153,18 +1359,24 @@ class UVBeam(UVBase):
         and then converting to power gives a different result than converting
         to power and then interpolating at about a 5% level.
 
-        Args:
-            nside: The nside to use for the Healpix map. If not specified, use
-                the nside that gives the closest resolution that is higher than the
-                input resolution.
-            run_check: Option to check for the existence and proper shapes of
-                required parameters after converting to healpix. Default is True.
-            check_extra: Option to check optional parameters as well as
-                required ones. Default is True.
-            run_check_acceptability: Option to check acceptable range of the values of
-                required parameters after combining objects. Default is True.
-            inplace: Option to perform the select directly on self (True, default) or return
-                a new UVBeam object, which is a subselection of self (False)
+        Parameters
+        ----------
+        nside : int
+            The nside to use for the Healpix map. If not specified, use
+            the nside that gives the closest resolution that is higher than the
+            input resolution.
+        run_check : bool
+            Option to check for the existence and proper shapes of required
+            parameters after converting to healpix.
+        check_extra : bool
+            Option to check optional parameters as well as required ones.
+        run_check_acceptability : bool
+            Option to check acceptable range of the values of required parameters
+            after combining objects
+        inplace : bool
+            Option to perform the select directly on self or return a new UVBeam
+            object.
+
         """
         try:
             import healpy as hp
@@ -1172,13 +1384,9 @@ class UVBeam(UVBase):
             uvutils._reraise_context('healpy is not installed but is required for '
                                      'healpix functionality')
 
-        if inplace:
-            beam_object = self
-        else:
-            beam_object = copy.deepcopy(self)
-
         if nside is None:
-            min_res = np.min(np.array([np.diff(beam_object.axis1_array)[0], np.diff(beam_object.axis2_array)[0]]))
+            min_res = np.min(np.array([np.diff(self.axis1_array)[0],
+                                       np.diff(self.axis2_array)[0]]))
             nside_min_res = np.sqrt(3 / np.pi) * np.radians(60.) / min_res
             nside = int(2**np.ceil(np.log2(nside_min_res)))
             assert(hp.pixelfunc.nside2resol(nside) < min_res)
@@ -1186,7 +1394,7 @@ class UVBeam(UVBase):
         npix = hp.nside2npix(nside)
         hpx_res = hp.pixelfunc.nside2resol(nside)
 
-        if np.iscomplexobj(beam_object.data_array):
+        if np.iscomplexobj(self.data_array):
             data_type = np.complex
         else:
             data_type = np.float
@@ -1207,42 +1415,18 @@ class UVBeam(UVBase):
 
         if inds_to_use.size < npix:
             pixels = pixels[inds_to_use]
-            hpx_theta = hpx_theta[inds_to_use]
-            hpx_phi = hpx_phi[inds_to_use]
-
-        interp_data, interp_basis_vector = \
-            self.interp(az_array=hpx_phi, za_array=hpx_theta)
-
-        beam_object.pixel_coordinate_system = 'healpix'
-        beam_object.nside = nside
-        beam_object.Npixels = npix
-        beam_object.ordering = 'ring'
-        beam_object.set_cs_params()
-
-        if beam_object.basis_vector_array is not None:
-            beam_object.basis_vector_array = interp_basis_vector
-
-        beam_object.pixel_array = pixels
-        beam_object.Npixels = beam_object.pixel_array.size
-        beam_object.data_array = interp_data
 
-        beam_object.Naxes1 = None
-        beam_object.Naxes2 = None
-        beam_object.axis1_array = None
-        beam_object.axis2_array = None
-
-        history_update_string = (' Interpolated from regularly gridded '
-                                 'azimuth/zenith_angle to HEALPix using pyuvdata '
-                                 'with interpolation_function = '
-                                 + self.interpolation_function + '.')
-
-        beam_object.history = beam_object.history + history_update_string
+        beam_object = self.interp(healpix_nside=nside, healpix_inds=pixels,
+                                  new_object=True, run_check=run_check,
+                                  check_extra=check_extra,
+                                  run_check_acceptability=run_check_acceptability)
 
-        if run_check:
-            beam_object.check(check_extra=check_extra,
-                              run_check_acceptability=run_check_acceptability)
         if not inplace:
             return beam_object
+        else:
+            for p in beam_object:
+                param = getattr(beam_object, p)
+                setattr(self, p, param)
 
     def __add__(self, other, run_check=True, check_extra=True,
                 run_check_acceptability=True, inplace=False):
diff --git a/pyuvdata/uvdata.py b/pyuvdata/uvdata.py
index d977b48b..e9fafb46 100644
--- a/pyuvdata/uvdata.py
+++ b/pyuvdata/uvdata.py
@@ -946,7 +946,7 @@ class UVData(UVBase):
                              'reflect the phasing status of the data')
 
         if not isinstance(time, Time):
-            raise(TypeError, "time must be an astropy.time.Time object")
+            raise TypeError("time must be an astropy.time.Time object")
 
         # Generate ra/dec of zenith at time in the phase_frame coordinate system
         # to use for phasing
@@ -1963,7 +1963,7 @@ class UVData(UVBase):
         if bls is not None:
             if isinstance(bls, tuple) and (len(bls) == 2 or len(bls) == 3):
                 bls = [bls]
-            if not all(isinstance(item, tuple) for item in bls):
+            if len(bls) == 0 or not all(isinstance(item, tuple) for item in bls):
                 raise ValueError(
                     'bls must be a list of tuples of antenna numbers (optionally with polarization).')
             if not all([isinstance(item[0], six.integer_types + (np.integer,)) for item in bls]
diff --git a/pyuvdata/uvflag.py b/pyuvdata/uvflag.py
index e183352b..06bafec3 100644
--- a/pyuvdata/uvflag.py
+++ b/pyuvdata/uvflag.py
@@ -612,7 +612,7 @@ class UVFlag(UVBase):
         if self.type == "waterfall":
             if antenna_nums is not None:
                 raise ValueError('Cannot select on antenna_nums with waterfall type '
-                                 'UVFLag objects.')
+                                 'UVFlag objects.')
             if bls is not None:
                 raise ValueError('Cannot select on bls with waterfall type '
                                  'UVFlag objects.')


Test patch:
diff --git a/pyuvdata/tests/__init__.py b/pyuvdata/tests/__init__.py
index b82d74cb..9ee709f3 100644
--- a/pyuvdata/tests/__init__.py
+++ b/pyuvdata/tests/__init__.py
@@ -86,10 +86,42 @@ def clearWarnings():
             reg.clear()
 
 
-def checkWarnings(func, func_args=[], func_kwargs={},
-                  category=UserWarning,
-                  nwarnings=1, message=None, known_warning=None):
-    """Function to check expected warnings."""
+def checkWarnings(func, func_args=[], func_kwargs={}, nwarnings=1,
+                  category=UserWarning, message=None, known_warning=None):
+    """Function to check expected warnings in tests.
+
+    Useful for checking that appropriate warnings are raised and to capture
+    (and silence) warnings in tests.
+
+    Parameters
+    ----------
+    func : function
+        Function or method to check warnings for.
+    func_args : list, optional
+        List of positional parameters to pass `func`
+    func_kwargs : dict, optional
+        Dict of keyword parameter to pass func. Keys are the parameter names,
+        values are the values to pass to the parameters.
+    nwarnings : int
+        Number of expected warnings.
+    category : warning type or list of warning types
+        Expected warning type(s). If a scalar is passed and `nwarnings` is
+        greater than one, the same category will be expected for all warnings.
+    message : str or list of str
+        Expected warning string(s). If a scalar is passed and `nwarnings` is
+        greater than one, the same warning string will be expected for all warnings.
+    known_warning : {'miriad', 'paper_uvfits', 'fhd'}, optional
+        Shorthand way to specify one of a standard set of warnings.
+
+    Returns
+    -------
+    Value returned by `func`
+
+    Raises
+    ------
+    AssertionError
+        If the warning(s) raised by func do not match the expected values.
+    """
 
     if (not isinstance(category, list) or len(category) == 1) and nwarnings > 1:
         if isinstance(category, list):
diff --git a/pyuvdata/tests/test_uvbeam.py b/pyuvdata/tests/test_uvbeam.py
index 60cade96..1d52d14c 100644
--- a/pyuvdata/tests/test_uvbeam.py
+++ b/pyuvdata/tests/test_uvbeam.py
@@ -161,9 +161,9 @@ def test_properties(uvbeam_data):
         this_param = getattr(uvbeam_data.beam_obj, v)
         try:
             assert rand_num == this_param.value
-        except(AssertionError):
+        except AssertionError:
             print('setting {prop_name} to a random number failed'.format(prop_name=k))
-            raise(AssertionError)
+            raise
 
 
 def test_errors():
@@ -391,31 +391,72 @@ def test_freq_interpolation():
                              feed_version='0.1', feed_pol=['x'],
                              model_name='E-field pattern - Rigging height 4.9m',
                              model_version='1.0')
+    power_beam.interpolation_function = 'az_za_simple'
 
     # test frequency interpolation returns data arrays for small and large tolerances
     freq_orig_vals = np.array([123e6, 150e6])
-    interp_data, interp_bandpass = power_beam._interp_freq(freq_orig_vals, tol=0.0, new_object=False)
+    interp_data, interp_basis_vector, interp_bandpass = \
+        power_beam.interp(freq_array=freq_orig_vals, freq_interp_tol=0.0,
+                          return_bandpass=True)
     assert isinstance(interp_data, np.ndarray)
     assert isinstance(interp_bandpass, np.ndarray)
+    np.testing.assert_array_almost_equal(power_beam.bandpass_array, interp_bandpass)
+    np.testing.assert_array_almost_equal(power_beam.data_array, interp_data)
+    assert interp_basis_vector is None
 
-    interp_data, interp_bandpass = power_beam._interp_freq(freq_orig_vals, tol=1.0, new_object=False)
+    interp_data, interp_basis_vector, interp_bandpass = \
+        power_beam.interp(freq_array=freq_orig_vals, freq_interp_tol=1.0,
+                          return_bandpass=True)
     assert isinstance(interp_data, np.ndarray)
     assert isinstance(interp_bandpass, np.ndarray)
+    np.testing.assert_array_almost_equal(power_beam.bandpass_array, interp_bandpass)
+    np.testing.assert_array_almost_equal(power_beam.data_array, interp_data)
+    assert interp_basis_vector is None
 
     # test frequency interpolation returns new UVBeam for small and large tolerances
     power_beam.saved_interp_functions = {}
-    interp_data, interp_bandpass = power_beam._interp_freq(freq_orig_vals, tol=0.0, new_object=True)
-    assert isinstance(interp_data, UVBeam)
-    assert interp_bandpass is None
-    np.testing.assert_array_almost_equal(interp_data.freq_array[0], freq_orig_vals)
-    assert interp_data.freq_interp_kind == 'linear'
-    assert not hasattr(interp_data, 'saved_interp_functions')  # test that saved functions are erased in new obj
-
-    interp_data, interp_bandpass = power_beam._interp_freq(freq_orig_vals, tol=1.0, new_object=True)
-    assert isinstance(interp_data, UVBeam)
-    assert interp_bandpass is None
-    np.testing.assert_array_almost_equal(interp_data.freq_array[0], freq_orig_vals)
-    assert interp_data.freq_interp_kind == 'nearest'  # assert interp kind is 'nearest' when within tol
+    new_beam_obj = power_beam.interp(freq_array=freq_orig_vals, freq_interp_tol=0.0,
+                                     new_object=True)
+    assert isinstance(new_beam_obj, UVBeam)
+    np.testing.assert_array_almost_equal(new_beam_obj.freq_array[0], freq_orig_vals)
+    assert new_beam_obj.freq_interp_kind == 'linear'
+    # test that saved functions are erased in new obj
+    assert not hasattr(new_beam_obj, 'saved_interp_functions')
+    assert power_beam.history != new_beam_obj.history
+    new_beam_obj.history = power_beam.history
+    assert power_beam == new_beam_obj
+
+    new_beam_obj = power_beam.interp(freq_array=freq_orig_vals, freq_interp_tol=1.0,
+                                     new_object=True)
+    assert isinstance(new_beam_obj, UVBeam)
+    np.testing.assert_array_almost_equal(new_beam_obj.freq_array[0], freq_orig_vals)
+    # assert interp kind is 'nearest' when within tol
+    assert new_beam_obj.freq_interp_kind == 'nearest'
+    new_beam_obj.freq_interp_kind = 'linear'
+    assert power_beam.history != new_beam_obj.history
+    new_beam_obj.history = power_beam.history
+    assert power_beam == new_beam_obj
+
+    # test frequency interpolation returns valid new UVBeam for different
+    # number of freqs from input
+    power_beam.saved_interp_functions = {}
+    new_beam_obj = power_beam.interp(freq_array=np.linspace(123e6, 150e6, num=5),
+                                     freq_interp_tol=0.0, new_object=True)
+
+    assert isinstance(new_beam_obj, UVBeam)
+    np.testing.assert_array_almost_equal(new_beam_obj.freq_array[0],
+                                         np.linspace(123e6, 150e6, num=5))
+    assert new_beam_obj.freq_interp_kind == 'linear'
+    # test that saved functions are erased in new obj
+    assert not hasattr(new_beam_obj, 'saved_interp_functions')
+    assert power_beam.history != new_beam_obj.history
+    new_beam_obj.history = power_beam.history
+
+    # down select to orig freqs and test equality
+    new_beam_obj.select(frequencies=freq_orig_vals)
+    assert power_beam.history != new_beam_obj.history
+    new_beam_obj.history = power_beam.history
+    assert power_beam == new_beam_obj
 
     # using only one freq chan should trigger a ValueError if interp_bool is True
     # unless requesting the original frequency channel such that interp_bool is False.
@@ -424,13 +465,13 @@ def test_freq_interpolation():
     # Other ways of testing this (e.g. interp_data_array.flags['OWNDATA']) does not work
     _pb = power_beam.select(frequencies=power_beam.freq_array[0, :1], inplace=False)
     try:
-        interp_data_array, interp_bandp = _pb._interp_freq(_pb.freq_array[0], kind='linear')
+        interp_data, interp_basis_vector = _pb.interp(freq_array=_pb.freq_array[0])
     except ValueError:
-        raise AssertionError("UVBeam._interp_freq didn't return an array slice as expected")
+        raise AssertionError("UVBeam.interp didn't return an array slice as expected")
 
     # test errors if one frequency
     power_beam_singlef = power_beam.select(freq_chans=[0], inplace=False)
-    pytest.raises(ValueError, power_beam_singlef._interp_freq, np.array([150e6]))
+    pytest.raises(ValueError, power_beam_singlef.interp, freq_array=np.array([150e6]))
 
     # assert freq_interp_kind ValueError
     power_beam.interpolation_function = 'az_za_simple'
@@ -469,6 +510,42 @@ def test_power_spatial_interpolation():
 
     assert np.allclose(data_array_compare, interp_data_array)
 
+    # test that new object from interpolation is identical
+    new_power_beam = power_beam.interp(az_array=power_beam.axis1_array,
+                                       za_array=power_beam.axis2_array,
+                                       az_za_grid=True,
+                                       freq_array=freq_orig_vals,
+                                       new_object=True)
+    assert new_power_beam.freq_interp_kind == 'nearest'
+    assert new_power_beam.history == (power_beam.history + ' Interpolated in '
+                                      'frequency and to a new azimuth/zenith '
+                                      'angle grid using pyuvdata with '
+                                      'interpolation_function = az_za_simple '
+                                      'and freq_interp_kind = nearest.')
+    # make histories & freq_interp_kind equal
+    new_power_beam.history = power_beam.history
+    new_power_beam.freq_interp_kind = 'linear'
+    assert new_power_beam == power_beam
+
+    # test that interp to every other point returns an object that matches a select
+    axis1_inds = np.arange(0, power_beam.Naxes1, 2)
+    axis2_inds = np.arange(0, power_beam.Naxes2, 2)
+
+    select_beam = power_beam.select(axis1_inds=axis1_inds, axis2_inds=axis2_inds,
+                                    inplace=False)
+    interp_beam = power_beam.interp(az_array=power_beam.axis1_array[axis1_inds],
+                                    za_array=power_beam.axis2_array[axis2_inds],
+                                    az_za_grid=True, new_object=True)
+    assert select_beam.history != interp_beam.history
+    interp_beam.history = select_beam.history
+    assert select_beam == interp_beam
+
+    # test error if new_object set without az_za_grid
+    with pytest.raises(ValueError) as cm:
+        power_beam.interp(az_array=az_orig_vals, za_array=za_orig_vals,
+                          freq_array=freq_orig_vals, new_object=True)
+    assert str(cm.value).startswith('A new object can only be returned')
+
     # test only a single polarization
     interp_data_array, interp_basis_vector = power_beam.interp(az_array=az_orig_vals,
                                                                za_array=za_orig_vals,
@@ -489,11 +566,13 @@ def test_power_spatial_interpolation():
     # Test requesting separate polarizations on different calls while reusing splines.
     interp_data_array, interp_basis_vector = power_beam.interp(az_array=az_interp_vals[:2],
                                                                za_array=za_interp_vals[:2],
-                                                               freq_array=freq_interp_vals, polarizations=['xx'], reuse_spline=True)
+                                                               freq_array=freq_interp_vals,
+                                                               polarizations=['xx'], reuse_spline=True)
 
     interp_data_array, interp_basis_vector = power_beam.interp(az_array=az_interp_vals[:2],
                                                                za_array=za_interp_vals[:2],
-                                                               freq_array=freq_interp_vals, polarizations=['yy'], reuse_spline=True)
+                                                               freq_array=freq_interp_vals,
+                                                               polarizations=['yy'], reuse_spline=True)
 
     # test reusing the spline fit.
     orig_data_array, interp_basis_vector = power_beam.interp(az_array=az_interp_vals,
@@ -549,6 +628,35 @@ def test_efield_spatial_interpolation():
     assert np.allclose(efield_beam.data_array, interp_data_array)
     assert np.allclose(efield_beam.basis_vector_array, interp_basis_vector)
 
+    # test that new object from interpolation is identical
+    new_efield_beam = efield_beam.interp(az_array=efield_beam.axis1_array,
+                                         za_array=efield_beam.axis2_array,
+                                         az_za_grid=True,
+                                         freq_array=freq_orig_vals,
+                                         new_object=True)
+    assert new_efield_beam.freq_interp_kind == 'nearest'
+    assert new_efield_beam.history == (efield_beam.history + ' Interpolated in '
+                                       'frequency and to a new azimuth/zenith '
+                                       'angle grid using pyuvdata with '
+                                       'interpolation_function = az_za_simple '
+                                       'and freq_interp_kind = nearest.')
+    # make histories & freq_interp_kind equal
+    new_efield_beam.history = efield_beam.history
+    new_efield_beam.freq_interp_kind = 'linear'
+    assert new_efield_beam == efield_beam
+
+    # test that interp to every other point returns an object that matches a select
+    axis1_inds = np.arange(0, efield_beam.Naxes1, 2)
+    axis2_inds = np.arange(0, efield_beam.Naxes2, 2)
+    select_beam = efield_beam.select(axis1_inds=axis1_inds, axis2_inds=axis2_inds,
+                                     inplace=False)
+    interp_beam = efield_beam.interp(az_array=efield_beam.axis1_array[axis1_inds],
+                                     za_array=efield_beam.axis2_array[axis2_inds],
+                                     az_za_grid=True, new_object=True)
+    assert select_beam.history != interp_beam.history
+    interp_beam.history = select_beam.history
+    assert select_beam == interp_beam
+
     # test no errors using different points
     az_interp_vals = np.array(np.arange(0, 2 * np.pi, np.pi / 9.0).tolist()
                               + np.arange(0, 2 * np.pi, np.pi / 9.0).tolist())
@@ -577,7 +685,8 @@ def test_efield_spatial_interpolation():
 
     select_data_array_reused, interp_basis_vector = efield_beam.interp(az_array=az_interp_vals[0:1],
                                                                        za_array=za_interp_vals[0:1],
-                                                                       freq_array=np.array([127e6]), reuse_spline=True)
+                                                                       freq_array=np.array([127e6]),
+                                                                       reuse_spline=True)
     assert np.allclose(select_data_array_orig, select_data_array_reused)
     del efield_beam.saved_interp_functions
 
@@ -639,35 +748,103 @@ def test_interp_longitude_branch_cut():
 
 @uvtest.skipIf_no_healpy
 def test_healpix_interpolation():
-    power_beam = UVBeam()
-    power_beam.read_cst_beam(cst_files, beam_type='efield', frequency=[150e6, 123e6],
-                             telescope_name='TEST', feed_name='bob',
-                             feed_version='0.1', feed_pol=['x'],
-                             model_name='E-field pattern - Rigging height 4.9m',
-                             model_version='1.0')
-    power_beam.interpolation_function = 'az_za_simple'
-    power_beam.to_healpix()
+    efield_beam = UVBeam()
+    efield_beam.read_cst_beam(cst_files, beam_type='efield', frequency=[150e6, 123e6],
+                              telescope_name='TEST', feed_name='bob',
+                              feed_version='0.1', feed_pol=['x'],
+                              model_name='E-field pattern - Rigging height 4.9m',
+                              model_version='1.0')
+    efield_beam.interpolation_function = 'az_za_simple'
+    orig_efield_beam = copy.deepcopy(efield_beam)
+
+    # test calling interp with healpix parameters directly gives same result
+    min_res = np.min(np.array([np.diff(efield_beam.axis1_array)[0],
+                               np.diff(efield_beam.axis2_array)[0]]))
+    nside_min_res = np.sqrt(3 / np.pi) * np.radians(60.) / min_res
+    nside = int(2**np.ceil(np.log2(nside_min_res)))
+
+    new_efield_beam = efield_beam.interp(healpix_nside=nside,
+                                         new_object=True)
+    efield_beam.to_healpix()
+    assert new_efield_beam == efield_beam
 
     # check that interpolating to existing points gives the same answer
-    power_beam.interpolation_function = 'healpix_simple'
-    za_orig_vals, az_orig_vals = hp.pix2ang(power_beam.nside, np.arange(hp.nside2npix(power_beam.nside)))
+    efield_beam.interpolation_function = 'healpix_simple'
+    za_orig_vals, az_orig_vals = hp.pix2ang(efield_beam.nside, np.arange(hp.nside2npix(efield_beam.nside)))
     az_orig_vals = az_orig_vals.ravel(order='C')
     za_orig_vals = za_orig_vals.ravel(order='C')
     freq_orig_vals = np.array([123e6, 150e6])
 
-    interp_data_array, interp_basis_vector = power_beam.interp(az_array=az_orig_vals,
-                                                               za_array=za_orig_vals,
-                                                               freq_array=freq_orig_vals)
-    data_array_compare = power_beam.data_array
+    interp_data_array, interp_basis_vector = efield_beam.interp(az_array=az_orig_vals,
+                                                                za_array=za_orig_vals,
+                                                                freq_array=freq_orig_vals)
+    data_array_compare = efield_beam.data_array
     interp_data_array = interp_data_array.reshape(data_array_compare.shape, order='F')
     assert np.allclose(data_array_compare, interp_data_array)
 
+    # test calling interp with healpix parameters directly gives same result
+    new_efield_beam = efield_beam.interp(healpix_nside=efield_beam.nside,
+                                         freq_array=freq_orig_vals,
+                                         new_object=True)
+    assert new_efield_beam.freq_interp_kind == 'nearest'
+    assert new_efield_beam.history == (efield_beam.history + ' Interpolated in '
+                                       'frequency and to a new healpix grid '
+                                       'using pyuvdata with '
+                                       'interpolation_function = healpix_simple '
+                                       'and freq_interp_kind = nearest.')
+    # make histories & freq_interp_kind equal
+    new_efield_beam.history = efield_beam.history
+    new_efield_beam.freq_interp_kind = 'linear'
+    assert new_efield_beam == efield_beam
+    del(new_efield_beam)
+
+    # test that interp to every other point returns an object that matches a select
+    pixel_inds = np.arange(0, efield_beam.Npixels, 2)
+    select_beam = efield_beam.select(pixels=pixel_inds, inplace=False)
+    interp_beam = efield_beam.interp(healpix_inds=efield_beam.pixel_array[pixel_inds],
+                                     healpix_nside=efield_beam.nside, new_object=True)
+    assert select_beam.history != interp_beam.history
+    interp_beam.history = select_beam.history
+    assert select_beam == interp_beam
+
+    # test interp from healpix to regular az/za grid
+    new_reg_beam = efield_beam.interp(az_array=orig_efield_beam.axis1_array,
+                                      za_array=orig_efield_beam.axis2_array,
+                                      az_za_grid=True, new_object=True)
+
+    # this diff is pretty large. 2 rounds of interpolation is not a good thing.
+    # but we can check that the rest of the object makes sense
+    diff = new_reg_beam.data_array - orig_efield_beam.data_array
+    diff_ratio = diff / orig_efield_beam.data_array
+    assert np.all(np.abs(diff_ratio) < 3)
+    # set data_array tolerances higher to test the rest of the object
+    # tols are (relative, absolute)
+    tols = [3, 0]
+    new_reg_beam._data_array.tols = tols
+    assert new_reg_beam.history != orig_efield_beam.history
+    new_reg_beam.history = orig_efield_beam.history
+    new_reg_beam.interpolation_function = 'az_za_simple'
+    assert new_reg_beam == orig_efield_beam
+
+    # test errors with specifying healpix_inds without healpix_nside
+    with pytest.raises(ValueError) as cm:
+        efield_beam.interp(healpix_inds=np.arange(hp.nside2npix(efield_beam.nside)),
+                           freq_array=freq_orig_vals)
+    assert str(cm.value).startswith('healpix_nside must be set if healpix_inds is set')
+
+    # test error setting both healpix_nside and az_array
+    with pytest.raises(ValueError) as cm:
+        efield_beam.interp(healpix_nside=efield_beam.nside, az_array=az_orig_vals,
+                           za_array=za_orig_vals, freq_array=freq_orig_vals)
+    assert str(cm.value).startswith('healpix_nside and healpix_inds can not be')
+
     # basis_vector exception
-    power_beam.basis_vector_array[0, 1, :] = 10.0
-    pytest.raises(NotImplementedError, power_beam.interp, az_array=az_orig_vals, za_array=za_orig_vals)
+    efield_beam.basis_vector_array[0, 1, :] = 10.0
+    pytest.raises(NotImplementedError, efield_beam.interp, az_array=az_orig_vals, za_array=za_orig_vals)
 
     # now convert to power beam
-    power_beam.efield_to_power()
+    power_beam = efield_beam.efield_to_power(inplace=False)
+    del(efield_beam)
     interp_data_array, interp_basis_vector = power_beam.interp(az_array=az_orig_vals,
                                                                za_array=za_orig_vals,
                                                                freq_array=freq_orig_vals)
@@ -675,6 +852,15 @@ def test_healpix_interpolation():
     interp_data_array = interp_data_array.reshape(data_array_compare.shape, order='F')
     assert np.allclose(data_array_compare, interp_data_array)
 
+    # test that interp to every other point returns an object that matches a select
+    pixel_inds = np.arange(0, power_beam.Npixels, 2)
+    select_beam = power_beam.select(pixels=pixel_inds, inplace=False)
+    interp_beam = power_beam.interp(healpix_inds=power_beam.pixel_array[pixel_inds],
+                                    healpix_nside=power_beam.nside, new_object=True)
+    assert select_beam.history != interp_beam.history
+    interp_beam.history = select_beam.history
+    assert select_beam == interp_beam
+
     # assert not feeding frequencies gives same answer
     interp_data_array2, interp_basis_vector2 = power_beam.interp(az_array=az_orig_vals, za_array=za_orig_vals)
     assert np.allclose(interp_data_array, interp_data_array2)
@@ -689,6 +875,7 @@ def test_healpix_interpolation():
     assert np.allclose(interp_data_array[:, :, 1:2], interp_data_array2[:, :, :1])
 
     # change complex data_array to real data_array and test again
+    assert power_beam.data_array.dtype == np.complex
     power_beam.data_array = np.abs(power_beam.data_array)
     interp_data_array, interp_basis_vector = power_beam.interp(az_array=az_orig_vals,
                                                                za_array=za_orig_vals,
@@ -726,9 +913,12 @@ def test_to_healpix():
 
     # check that history is updated appropriately
     assert power_beam_healpix.history == (power_beam.history
-                                          + ' Interpolated from regularly gridded '
-                                          'azimuth/zenith_angle to HEALPix using pyuvdata '
-                                          'with interpolation_function = az_za_simple.')
+                                          + ' Interpolated from '
+                                          + power_beam.coordinate_system_dict['az_za']['description']
+                                          + ' to '
+                                          + power_beam.coordinate_system_dict['healpix']['description']
+                                          + ' using pyuvdata with '
+                                          'interpolation_function = az_za_simple.')
 
     npix = hp.nside2npix(power_beam_healpix.nside)
     assert power_beam_healpix.Npixels <= npix * 0.55
@@ -767,9 +957,11 @@ def test_to_healpix():
     sq_then_interp._data_array.tols = tols
 
     # check history changes
-    interp_history_add = (' Interpolated from regularly gridded '
-                          'azimuth/zenith_angle to HEALPix using pyuvdata '
-                          'with interpolation_function = az_za_simple.')
+    interp_history_add = (' Interpolated from '
+                          + power_beam.coordinate_system_dict['az_za']['description']
+                          + ' to '
+                          + power_beam.coordinate_system_dict['healpix']['description']
+                          + ' using pyuvdata with interpolation_function = az_za_simple.')
     sq_history_add = ' Converted from efield to power using pyuvdata.'
     assert sq_then_interp.history == efield_beam.history + sq_history_add + interp_history_add
     assert interp_then_sq.history == efield_beam.history + interp_history_add + sq_history_add
diff --git a/pyuvdata/tests/test_uvcal.py b/pyuvdata/tests/test_uvcal.py
index 5f8b14bc..c4b12ee3 100644
--- a/pyuvdata/tests/test_uvcal.py
+++ b/pyuvdata/tests/test_uvcal.py
@@ -128,9 +128,9 @@ def test_properties(uvcal_data):
         this_param = getattr(uv_cal_object, v)
         try:
             assert rand_num == this_param.value
-        except(AssertionError):
+        except AssertionError:
             print('setting {prop_name} to a random number failed'.format(prop_name=k))
-            raise(AssertionError)
+            raise
 
 
 @pytest.fixture(scope='function')
diff --git a/pyuvdata/tests/test_uvdata.py b/pyuvdata/tests/test_uvdata.py
index 4ee34a7c..65e1dbbb 100644
--- a/pyuvdata/tests/test_uvdata.py
+++ b/pyuvdata/tests/test_uvdata.py
@@ -164,9 +164,9 @@ def test_properties(uvdata_props):
         this_param = getattr(uvdata_props.uv_object, v)
         try:
             assert rand_num == this_param.value
-        except(AssertionError):
+        except AssertionError:
             print('setting {prop_name} to a random number failed'.format(prop_name=k))
-            raise(AssertionError)
+            raise
 
 
 @pytest.fixture(scope='function')
@@ -494,7 +494,9 @@ def test_phase_unphaseHERA():
     assert np.allclose(UV_phase_simple_small.uvw_array, UV_raw_small.uvw_array, atol=1e-2)
 
     # check error if not passing a Time object to phase_to_time
-    pytest.raises(TypeError, UV_raw.phase_to_time, UV_raw.time_array[0])
+    with pytest.raises(TypeError) as cm:
+        UV_raw.phase_to_time(UV_raw.time_array[0])
+    assert str(cm.value).startswith('time must be an astropy.time.Time object')
 
     # check errors when trying to unphase drift or unknown data
     pytest.raises(ValueError, UV_raw.unphase_to_drift)
@@ -893,15 +895,36 @@ def test_select_bls():
     assert list(set(sorted_pairs_object2)) == [(0, 6)]
 
     # check for errors associated with antenna pairs not included in data and bad inputs
-    pytest.raises(ValueError, uv_object.select,
-                  bls=list(zip(first_ants, second_ants)) + [0, 6])
-    pytest.raises(ValueError, uv_object.select,
-                  bls=[(uv_object.antenna_names[0], uv_object.antenna_names[1])])
-    pytest.raises(ValueError, uv_object.select, bls=(5, 1))
-    pytest.raises(ValueError, uv_object.select, bls=(0, 5))
-    pytest.raises(ValueError, uv_object.select, bls=(27, 27))
-    pytest.raises(ValueError, uv_object.select, bls=(6, 0, 'RR'), polarizations='RR')
-    pytest.raises(ValueError, uv_object.select, bls=(6, 0, 8))
+    with pytest.raises(ValueError) as cm:
+        uv_object.select(bls=list(zip(first_ants, second_ants)) + [0, 6])
+    assert str(cm.value).startswith('bls must be a list of tuples of antenna numbers')
+
+    with pytest.raises(ValueError) as cm:
+        uv_object.select(bls=[(uv_object.antenna_names[0], uv_object.antenna_names[1])])
+    assert str(cm.value).startswith('bls must be a list of tuples of antenna numbers')
+
+    with pytest.raises(ValueError) as cm:
+        uv_object.select(bls=(5, 1))
+    assert str(cm.value).startswith('Antenna number 5 is not present in the '
+                                    'ant_1_array or ant_2_array')
+    with pytest.raises(ValueError) as cm:
+        uv_object.select(bls=(0, 5))
+    assert str(cm.value).startswith('Antenna number 5 is not present in the '
+                                    'ant_1_array or ant_2_array')
+    with pytest.raises(ValueError) as cm:
+        uv_object.select(bls=(27, 27))
+    assert str(cm.value).startswith('Antenna pair (27, 27) does not have any data')
+    with pytest.raises(ValueError) as cm:
+        uv_object.select(bls=(6, 0, 'RR'), polarizations='RR')
+    assert str(cm.value).startswith('Cannot provide length-3 tuples and also '
+                                    'specify polarizations.')
+    with pytest.raises(ValueError) as cm:
+        uv_object.select(bls=(6, 0, 8))
+    assert str(cm.value).startswith('The third element in each bl must be a '
+                                    'polarization string')
+    with pytest.raises(ValueError) as cm:
+        uv_object.select(bls=[])
+    assert str(cm.value).startswith('bls must be a list of tuples of antenna numbers')
 
 
 def test_select_times():
diff --git a/pyuvdata/tests/test_uvflag.py b/pyuvdata/tests/test_uvflag.py
index 2003515d..16f7730a 100644
--- a/pyuvdata/tests/test_uvflag.py
+++ b/pyuvdata/tests/test_uvflag.py
@@ -66,7 +66,7 @@ def uvf_from_waterfall():
 
 
 # Try to import `pytest-cases` and define decorators used to
-# iterate over the three main types of UVFLag objects
+# iterate over the three main types of UVFlag objects
 # otherwise make the decorators skip the tests that use these iterators
 try:
     import pytest_cases