Merge b10a22d into 8b28feb

docs/tutorial.rst

@@ -327,7 +327,7 @@ a) Getting antenna positions in topocentric frame in units of meters
   >>> uvd = UVData()
   >>> uvd.read_miriad('pyuvdata/data/new.uvA')
   >>> antpos = uvd.antenna_positions + uvd.telescope_location # get antennas positions in ECEF
-  >>> antpos = uvutils.ENU_from_ECEF(antpos.T, *uvd.telescope_location_lat_lon_alt).T # convert to topo (ENU) coords.
+  >>> antpos = uvutils.ENU_from_ECEF(antpos, *uvd.telescope_location_lat_lon_alt) # convert to topo (ENU) coords.
 
 
 UVData: Selecting data

pyuvdata/tests/test_utils.py

@@ -7,12 +7,12 @@
 
 import os
 import nose.tools as nt
+import numpy as np
 from astropy import units
 from astropy.time import Time
 from astropy.coordinates import SkyCoord, Angle
 from astropy.io import fits
 import pyuvdata
-import numpy as np
 from pyuvdata.data import DATA_PATH
 import pyuvdata.utils as uvutils
 import pyuvdata.tests as uvtest
@@ -47,9 +47,20 @@ def test_LatLonAlt_from_XYZ():
     nt.assert_raises(ValueError, pyuvdata.LatLonAlt_from_XYZ, ref_latlonalt)
 
     # test passing multiple values
-    xyz_mult = np.stack((np.array(ref_xyz), np.array(ref_xyz)), axis=1)
+    xyz_mult = np.stack((np.array(ref_xyz), np.array(ref_xyz)))
     lat_vec, lon_vec, alt_vec = pyuvdata.LatLonAlt_from_XYZ(xyz_mult)
     nt.assert_true(np.allclose(ref_latlonalt, (lat_vec[1], lon_vec[1], alt_vec[1]), rtol=0, atol=1e-3))
+    # check warning if array transposed
+    uvtest.checkWarnings(pyuvdata.LatLonAlt_from_XYZ, [xyz_mult.T],
+                         message='The expected shape of ECEF xyz array',
+                         category=PendingDeprecationWarning)
+    # check warning if  3 x 3 array
+    xyz_3 = np.stack((np.array(ref_xyz), np.array(ref_xyz), np.array(ref_xyz)))
+    uvtest.checkWarnings(pyuvdata.LatLonAlt_from_XYZ, [xyz_3],
+                         message='The xyz array in LatLonAlt_from_XYZ is',
+                         category=PendingDeprecationWarning)
+    # check error if only 2 coordinates
+    nt.assert_raises(ValueError, pyuvdata.LatLonAlt_from_XYZ, xyz_mult[:, 0:2])
 
     # test error checking
     nt.assert_raises(ValueError, pyuvdata.LatLonAlt_from_XYZ, ref_xyz[0:1])
@@ -92,25 +103,51 @@ def test_ENU_tofrom_ECEF():
           -0.02019057, 0.16979185, 0.06945155, -0.64058124]
 
     xyz = pyuvdata.XYZ_from_LatLonAlt(lats, lons, alts)
-    nt.assert_true(np.allclose(np.stack((x, y, z)), xyz, atol=1e-3))
+    nt.assert_true(np.allclose(np.stack((x, y, z), axis=1), xyz, atol=1e-3))
 
     enu = pyuvdata.ENU_from_ECEF(xyz, center_lat, center_lon, center_alt)
-    nt.assert_true(np.allclose(np.stack((east, north, up)), enu, atol=1e-3))
+    nt.assert_true(np.allclose(np.stack((east, north, up), axis=1), enu, atol=1e-3))
+    # check warning if array transposed
+    uvtest.checkWarnings(pyuvdata.ENU_from_ECEF, [xyz.T, center_lat, center_lon,
+                                                  center_alt],
+                         message='The expected shape of ECEF xyz array',
+                         category=PendingDeprecationWarning)
+    # check warning if  3 x 3 array
+    uvtest.checkWarnings(pyuvdata.ENU_from_ECEF, [xyz[0:3], center_lat, center_lon,
+                                                  center_alt],
+                         message='The xyz array in ENU_from_ECEF is',
+                         category=PendingDeprecationWarning)
+    # check error if only 2 coordinates
+    nt.assert_raises(ValueError, pyuvdata.ENU_from_ECEF, xyz[:, 0:2],
+                     center_lat, center_lon, center_alt)
 
     # check that a round trip gives the original value.
     xyz_from_enu = pyuvdata.ECEF_from_ENU(enu, center_lat, center_lon, center_alt)
     nt.assert_true(np.allclose(xyz, xyz_from_enu, atol=1e-3))
+    # check warning if array transposed
+    uvtest.checkWarnings(pyuvdata.ECEF_from_ENU, [enu.T, center_lat, center_lon,
+                                                  center_alt],
+                         message='The expected shape the ENU array',
+                         category=PendingDeprecationWarning)
+    # check warning if  3 x 3 array
+    uvtest.checkWarnings(pyuvdata.ECEF_from_ENU, [enu[0:3], center_lat, center_lon,
+                                                  center_alt],
+                         message='The enu array in ECEF_from_ENU is',
+                         category=PendingDeprecationWarning)
+    # check error if only 2 coordinates
+    nt.assert_raises(ValueError, pyuvdata.ENU_from_ECEF, enu[:, 0:2], center_lat,
+                     center_lon, center_alt)
 
     # check passing a single value
-    enu_single = pyuvdata.ENU_from_ECEF(xyz[:, 0], center_lat, center_lon, center_alt)
-    nt.assert_true(np.allclose(np.stack((east[0], north[0], up[0])), enu[:, 0], atol=1e-3))
+    enu_single = pyuvdata.ENU_from_ECEF(xyz[0, :], center_lat, center_lon, center_alt)
+    nt.assert_true(np.allclose(np.array((east[0], north[0], up[0])), enu[0, :], atol=1e-3))
 
     xyz_from_enu = pyuvdata.ECEF_from_ENU(enu_single, center_lat, center_lon, center_alt)
-    nt.assert_true(np.allclose(xyz[:, 0], xyz_from_enu, atol=1e-3))
+    nt.assert_true(np.allclose(xyz[0, :], xyz_from_enu, atol=1e-3))
 
     # error checking
-    nt.assert_raises(ValueError, pyuvdata.ENU_from_ECEF, xyz[0:1, :], center_lat, center_lon, center_alt)
-    nt.assert_raises(ValueError, pyuvdata.ECEF_from_ENU, enu[0:1, :], center_lat, center_lon, center_alt)
+    nt.assert_raises(ValueError, pyuvdata.ENU_from_ECEF, xyz[:, 0:1], center_lat, center_lon, center_alt)
+    nt.assert_raises(ValueError, pyuvdata.ECEF_from_ENU, enu[:, 0:1], center_lat, center_lon, center_alt)
     nt.assert_raises(ValueError, pyuvdata.ENU_from_ECEF, xyz / 2., center_lat, center_lon, center_alt)
 
 
@@ -140,7 +177,6 @@ def test_mwa_ecef_conversion():
     # transpose these arrays to get them into the right shape
     txt_topo = txt_topo.T
     uvw_topo = uvw_topo.T
-    enh = enh.T
 
     # ARRAYX, ARRAYY, ARRAYZ in ECEF frame from Cotter file
     arrcent = f['arrcent']
@@ -152,7 +188,7 @@ def test_mwa_ecef_conversion():
     # add in array center to get real ECEF
     ecef_xyz = new_xyz + arrcent
 
-    enu = uvutils.ENU_from_ECEF(ecef_xyz.T, lat, lon, alt)
+    enu = uvutils.ENU_from_ECEF(ecef_xyz, lat, lon, alt)
 
     nt.assert_true(np.allclose(enu, enh))
 
@@ -179,7 +215,7 @@ def test_phasing_funcs():
     # in east/north/up frame (relative to array center) in meters: (Nants, 3)
     ants_enu = np.array([-101.94, 0156.41, 0001.24])
 
-    ant_xyz_abs = uvutils.ECEF_from_ENU(ants_enu.T, lat_lon_alt[0], lat_lon_alt[1], lat_lon_alt[2]).T
+    ant_xyz_abs = uvutils.ECEF_from_ENU(ants_enu, lat_lon_alt[0], lat_lon_alt[1], lat_lon_alt[2])
     ant_xyz_rel_itrs = ant_xyz_abs - array_center_xyz
     ant_xyz_rel_rot = uvutils.rotECEF_from_ECEF(ant_xyz_rel_itrs, lat_lon_alt[1])
 

pyuvdata/tests/test_uvdata.py

@@ -362,8 +362,8 @@ def test_phase_unphaseHERA():
 
     # check that they match if you phase & unphase using antenna locations
     # first replace the uvws with the right values
-    antenna_enu = uvutils.ENU_from_ECEF((UV_raw.antenna_positions + UV_raw.telescope_location).T,
-                                        *UV_raw.telescope_location_lat_lon_alt).T
+    antenna_enu = uvutils.ENU_from_ECEF((UV_raw.antenna_positions + UV_raw.telescope_location),
+                                        *UV_raw.telescope_location_lat_lon_alt)
     uvw_calc = np.zeros_like(UV_raw.uvw_array)
     unique_times, unique_inds = np.unique(UV_raw.time_array, return_index=True)
     for ind, jd in enumerate(unique_times):

pyuvdata/utils.py

@@ -37,35 +37,50 @@ def LatLonAlt_from_XYZ(xyz):
     Calculate lat/lon/alt from ECEF x,y,z.
 
     Args:
-        xyz: numpy array, shape (3, Npts), with ECEF x,y,z coordinates
+        xyz: numpy array, shape (Npts, 3), with ECEF x,y,z coordinates
 
     Returns:
-        tuple of latitude, longitude, altitude numpy arrays (if Npts > 1) or values (if Npts = 1) in radians & meters
+        tuple of latitude, longitude, altitude numpy arrays (if Npts > 1) or
+            values (if Npts = 1) in radians & meters
     """
     # convert to a numpy array
     xyz = np.array(xyz)
-    if xyz.shape[0] != 3:
-        raise ValueError(
-            'The first dimension of the ECEF xyz array must be length 3')
+    if xyz.ndim > 1 and xyz.shape[1] != 3:
+        if xyz.shape[0] == 3:
+            warnings.warn('The expected shape of ECEF xyz array is (Npts, 3). '
+                          'Support for arrays shaped (3, Npts) will go away in a '
+                          'future version.', PendingDeprecationWarning)
+            xyz_use = xyz.T
+        else:
+            raise ValueError('The expected shape of ECEF xyz array is (Npts, 3).')
+
+    else:
+        xyz_use = xyz
+
+    if xyz.shape == (3, 3):
+        warnings.warn('The xyz array in LatLonAlt_from_XYZ is being '
+                      'interpreted as (Npts, 3). Historically this function '
+                      'has supported (3, Npts) arrays, please verify that '
+                      'array ordering is as expected.', PendingDeprecationWarning)
 
-    xyz_use = xyz
     if xyz_use.ndim == 1:
-        xyz_use = xyz_use[:, np.newaxis]
+        xyz_use = xyz_use[np.newaxis, :]
 
     # checking for acceptable values
-    if np.any(np.linalg.norm(xyz, axis=0) < 6.35e6) or np.any(np.linalg.norm(xyz, axis=0) > 6.39e6):
+    if (np.any(np.linalg.norm(xyz_use, axis=1) < 6.35e6)
+            or np.any(np.linalg.norm(xyz_use, axis=1) > 6.39e6)):
         raise ValueError(
             'xyz values should be ECEF x, y, z coordinates in meters')
 
     # see wikipedia geodetic_datum and Datum transformations of
     # GPS positions PDF in docs/references folder
-    gps_p = np.sqrt(xyz_use[0, :]**2 + xyz_use[1, :]**2)
-    gps_theta = np.arctan2(xyz_use[2, :] * gps_a, gps_p * gps_b)
-    latitude = np.arctan2(xyz_use[2, :] + e_prime_squared * gps_b
+    gps_p = np.sqrt(xyz_use[:, 0]**2 + xyz_use[:, 1]**2)
+    gps_theta = np.arctan2(xyz_use[:, 2] * gps_a, gps_p * gps_b)
+    latitude = np.arctan2(xyz_use[:, 2] + e_prime_squared * gps_b
                           * np.sin(gps_theta)**3, gps_p - e_squared * gps_a
                           * np.cos(gps_theta)**3)
 
-    longitude = np.arctan2(xyz_use[1, :], xyz_use[0, :])
+    longitude = np.arctan2(xyz_use[:, 1], xyz_use[:, 0])
     gps_N = gps_a / np.sqrt(1 - e_squared * np.sin(latitude)**2)
     altitude = ((gps_p / np.cos(latitude)) - gps_N)
 
@@ -86,7 +101,7 @@ def XYZ_from_LatLonAlt(latitude, longitude, altitude):
         altitude: altitude in meters, can be a single value or a vector of length Npts
 
     Returns:
-        numpy array, shape (3, Npts) (if Npts > 1) or (3,) (if Npts = 1), with ECEF x,y,z coordinates
+        numpy array, shape (Npts, 3) (if Npts > 1) or (3,) (if Npts = 1), with ECEF x,y,z coordinates
     """
     latitude = np.array(latitude)
     longitude = np.array(longitude)
@@ -102,10 +117,10 @@ def XYZ_from_LatLonAlt(latitude, longitude, altitude):
     # see wikipedia geodetic_datum and Datum transformations of
     # GPS positions PDF in docs/references folder
     gps_N = gps_a / np.sqrt(1 - e_squared * np.sin(latitude)**2)
-    xyz = np.zeros((3, Npts))
-    xyz[0, :] = ((gps_N + altitude) * np.cos(latitude) * np.cos(longitude))
-    xyz[1, :] = ((gps_N + altitude) * np.cos(latitude) * np.sin(longitude))
-    xyz[2, :] = ((gps_b**2 / gps_a**2 * gps_N + altitude) * np.sin(latitude))
+    xyz = np.zeros((Npts, 3))
+    xyz[:, 0] = ((gps_N + altitude) * np.cos(latitude) * np.cos(longitude))
+    xyz[:, 1] = ((gps_N + altitude) * np.cos(latitude) * np.sin(longitude))
+    xyz[:, 2] = ((gps_b**2 / gps_a**2 * gps_N + altitude) * np.sin(latitude))
 
     xyz = np.squeeze(xyz)
     return xyz
@@ -155,47 +170,66 @@ def ENU_from_ECEF(xyz, latitude, longitude, altitude):
     Calculate local ENU (east, north, up) coordinates from ECEF coordinates.
 
     Args:
-        xyz: numpy array, shape (3, Npts), with ECEF x,y,z coordinates
+        xyz: numpy array, shape (Npts, 3), with ECEF x,y,z coordinates
         latitude: latitude of center of ENU coordinates in radians
         longitude: longitude of center of ENU coordinates in radians
         altitude: altitude of center of ENU coordinates in radians
 
     Returns:
-        numpy array, shape (3, Npts), with local ENU coordinates
+        numpy array, shape (Npts, 3), with local ENU coordinates
     """
-    if xyz.shape[0] != 3:
-        raise ValueError(
-            'The first dimension of the ECEF xyz array must be length 3')
+    xyz = np.array(xyz)
+    if xyz.ndim > 1 and xyz.shape[1] != 3:
+        if xyz.shape[0] == 3:
+            warnings.warn('The expected shape of ECEF xyz array is (Npts, 3). '
+                          'Support for arrays shaped (3, Npts) will go away in a '
+                          'future version.', PendingDeprecationWarning)
+            xyz_in = xyz.T
+            transpose = True
+        else:
+            raise ValueError('The expected shape of ECEF xyz array is (Npts, 3).')
+    else:
+        xyz_in = xyz
+        transpose = False
+
+    if xyz.shape == (3, 3):
+        warnings.warn('The xyz array in ENU_from_ECEF is being '
+                      'interpreted as (Npts, 3). Historically this function '
+                      'has supported (3, Npts) arrays, please verify that '
+                      'array ordering is as expected.', PendingDeprecationWarning)
+
+    if xyz_in.ndim == 1:
+        xyz_in = xyz_in[np.newaxis, :]
 
     # check that these are sensible ECEF values -- their magnitudes need to be
     # on the order of Earth's radius
-    ecef_magnitudes = np.linalg.norm(xyz, axis=0)
+    ecef_magnitudes = np.linalg.norm(xyz_in, axis=1)
     sensible_radius_range = (6.35e6, 6.39e6)
-    if np.any(ecef_magnitudes <= sensible_radius_range[0]) or np.any(ecef_magnitudes >= sensible_radius_range[1]):
+    if (np.any(ecef_magnitudes <= sensible_radius_range[0])
+            or np.any(ecef_magnitudes >= sensible_radius_range[1])):
         raise ValueError(
             'ECEF vector magnitudes must be on the order of the radius of the earth')
 
     xyz_center = XYZ_from_LatLonAlt(latitude, longitude, altitude)
 
-    xyz_in = xyz
-    if xyz_in.ndim == 1:
-        xyz_in = xyz_in[:, np.newaxis]
     xyz_use = np.zeros_like(xyz_in)
-    xyz_use[0, :] = xyz_in[0, :] - xyz_center[0]
-    xyz_use[1, :] = xyz_in[1, :] - xyz_center[1]
-    xyz_use[2, :] = xyz_in[2, :] - xyz_center[2]
+    xyz_use[:, 0] = xyz_in[:, 0] - xyz_center[0]
+    xyz_use[:, 1] = xyz_in[:, 1] - xyz_center[1]
+    xyz_use[:, 2] = xyz_in[:, 2] - xyz_center[2]
 
     enu = np.zeros_like(xyz_use)
-    enu[0, :] = (-np.sin(longitude) * xyz_use[0, :]
-                 + np.cos(longitude) * xyz_use[1, :])
-    enu[1, :] = (-np.sin(latitude) * np.cos(longitude) * xyz_use[0, :]
-                 - np.sin(latitude) * np.sin(longitude) * xyz_use[1, :]
-                 + np.cos(latitude) * xyz_use[2, :])
-    enu[2, :] = (np.cos(latitude) * np.cos(longitude) * xyz_use[0, :]
-                 + np.cos(latitude) * np.sin(longitude) * xyz_use[1, :]
-                 + np.sin(latitude) * xyz_use[2, :])
+    enu[:, 0] = (-np.sin(longitude) * xyz_use[:, 0]
+                 + np.cos(longitude) * xyz_use[:, 1])
+    enu[:, 1] = (-np.sin(latitude) * np.cos(longitude) * xyz_use[:, 0]
+                 - np.sin(latitude) * np.sin(longitude) * xyz_use[:, 1]
+                 + np.cos(latitude) * xyz_use[:, 2])
+    enu[:, 2] = (np.cos(latitude) * np.cos(longitude) * xyz_use[:, 0]
+                 + np.cos(latitude) * np.sin(longitude) * xyz_use[:, 1]
+                 + np.sin(latitude) * xyz_use[:, 2])
     if len(xyz.shape) == 1:
         enu = np.squeeze(enu)
+    elif transpose:
+        return enu.T
 
     return enu
 
@@ -205,36 +239,54 @@ def ECEF_from_ENU(enu, latitude, longitude, altitude):
     Calculate ECEF coordinates from local ENU (east, north, up) coordinates.
 
     Args:
-        enu: numpy array, shape (3, Npts), with local ENU coordinates
+        enu: numpy array, shape (Npts, 3), with local ENU coordinates
         latitude: latitude of center of ENU coordinates in radians
         longitude: longitude of center of ENU coordinates in radians
 
     Returns:
-        numpy array, shape (3, Npts), with ECEF x,y,z coordinates
-    """
-    if enu.shape[0] != 3:
-        raise ValueError(
-            'The first dimension of the local ENU array must be length 3')
+        numpy array, shape (Npts, 3), with ECEF x,y,z coordinates
+    """
+    enu = np.array(enu)
+    if enu.ndim > 1 and enu.shape[1] != 3:
+        if enu.shape[0] == 3:
+            warnings.warn('The expected shape of the ENU array is (Npts, 3). '
+                          'Support for arrays shaped (3, Npts) will go away in a '
+                          'future version.', PendingDeprecationWarning)
+            enu_use = enu.T
+            transpose = True
+        else:
+            raise ValueError('The expected shape of the ENU array array is (Npts, 3).')
+    else:
+        enu_use = enu
+        transpose = False
+
+    if enu.shape == (3, 3):
+        warnings.warn('The enu array in ECEF_from_ENU is being '
+                      'interpreted as (Npts, 3). Historically this function '
+                      'has supported (3, Npts) arrays, please verify that '
+                      'array ordering is as expected.', PendingDeprecationWarning)
 
-    enu_use = enu
     if enu_use.ndim == 1:
-        enu_use = enu_use[:, np.newaxis]
+        enu_use = enu_use[np.newaxis, :]
+
     xyz = np.zeros_like(enu_use)
-    xyz[0, :] = (-np.sin(latitude) * np.cos(longitude) * enu_use[1, :]
-                 - np.sin(longitude) * enu_use[0, :]
-                 + np.cos(latitude) * np.cos(longitude) * enu_use[2, :])
-    xyz[1, :] = (-np.sin(latitude) * np.sin(longitude) * enu_use[1, :]
-                 + np.cos(longitude) * enu_use[0, :]
-                 + np.cos(latitude) * np.sin(longitude) * enu_use[2, :])
-    xyz[2, :] = (np.cos(latitude) * enu_use[1, :]
-                 + np.sin(latitude) * enu_use[2, :])
+    xyz[:, 0] = (-np.sin(latitude) * np.cos(longitude) * enu_use[:, 1]
+                 - np.sin(longitude) * enu_use[:, 0]
+                 + np.cos(latitude) * np.cos(longitude) * enu_use[:, 2])
+    xyz[:, 1] = (-np.sin(latitude) * np.sin(longitude) * enu_use[:, 1]
+                 + np.cos(longitude) * enu_use[:, 0]
+                 + np.cos(latitude) * np.sin(longitude) * enu_use[:, 2])
+    xyz[:, 2] = (np.cos(latitude) * enu_use[:, 1]
+                 + np.sin(latitude) * enu_use[:, 2])
 
     xyz_center = XYZ_from_LatLonAlt(latitude, longitude, altitude)
-    xyz[0, :] = xyz[0, :] + xyz_center[0]
-    xyz[1, :] = xyz[1, :] + xyz_center[1]
-    xyz[2, :] = xyz[2, :] + xyz_center[2]
+    xyz[:, 0] = xyz[:, 0] + xyz_center[0]
+    xyz[:, 1] = xyz[:, 1] + xyz_center[1]
+    xyz[:, 2] = xyz[:, 2] + xyz_center[2]
     if len(enu.shape) == 1:
         xyz = np.squeeze(xyz)
+    elif transpose:
+        return xyz.T
 
     return xyz