https://github.com/RadioAstronomySoftwareGroup/pyuvdata
Raw File
Tip revision: 983ae27b57a28550b6ee1a58e13cc9161bb89327 authored by Bryna Hazelton on 15 January 2020, 19:15 UTC
update release date
Tip revision: 983ae27
Line_Profiler_full.out
This file contains the line by line timing results as run on the Oscar cluster at Brown. The task involved reading a 3.1G MIRIAD file, writing it out to miriad, phasing, writing to uvfits, then reading in the uvfits file and unphasing it.


Timer unit: 1e-06 s

Total time: 286.273 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/miriad.py
Function: read_miriad at line 28

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    28                                               @profile
    29                                               def read_miriad(self, filepath, correct_lat_lon=True, run_check=True, run_check_acceptability=True):
    30                                                   """
    31                                                   Read in data from a miriad file.
    32                                           
    33                                                   Args:
    34                                                       filepath: The miriad file directory to read from.
    35                                                       correct_lat_lon: flag -- that only matters if altitude is missing --
    36                                                           to update the latitude and longitude from the known_telescopes list
    37                                                       run_check: Option to check for the existence and proper shapes of
    38                                                           required parameters after reading in the file. Default is True.
    39                                                       run_check_acceptability: Option to check acceptable range of the values of
    40                                                           required parameters after reading in the file. Default is True.
    41                                                   """
    42         1        11715  11715.0      0.0          if not os.path.exists(filepath):
    43                                                       raise(IOError, filepath + ' not found')
    44         1        60578  60578.0      0.0          uv = a.miriad.UV(filepath)
    45                                           
    46         1            5      5.0      0.0          miriad_header_data = {'Nfreqs': 'nchan',
    47         1            4      4.0      0.0                                'Npols': 'npol',
    48         1            4      4.0      0.0                                'integration_time': 'inttime',
    49         1            3      3.0      0.0                                'channel_width': 'sdf',  # in Ghz!
    50         1            4      4.0      0.0                                'object_name': 'source',
    51                                                                         # NB: telescope_name and instrument are treated
    52                                                                         # as the same
    53         1            4      4.0      0.0                                'telescope_name': 'telescop',
    54         1            4      4.0      0.0                                'instrument': 'telescop'
    55                                                                         }
    56         8           32      4.0      0.0          for item in miriad_header_data:
    57         7           48      6.9      0.0              if isinstance(uv[miriad_header_data[item]], str):
    58         3           22      7.3      0.0                  header_value = uv[miriad_header_data[item]].replace('\x00', '')
    59                                                       else:
    60         4           24      6.0      0.0                  header_value = uv[miriad_header_data[item]]
    61         7           45      6.4      0.0              setattr(self, item, header_value)
    62                                           
    63         1            6      6.0      0.0          latitude = uv['latitud']  # in units of radians
    64         1            6      6.0      0.0          longitude = uv['longitu']
    65         1            4      4.0      0.0          try:
    66         1            5      5.0      0.0              altitude = uv['altitude']
    67         1           58     58.0      0.0              self.telescope_location_lat_lon_alt = (latitude, longitude, altitude)
    68                                                   except(KeyError):
    69                                                       # get info from known telescopes. Check to make sure the lat/lon values match reasonably well
    70                                                       telescope_obj = uvtel.get_telescope(self.telescope_name)
    71                                                       if telescope_obj is not False:
    72                                           
    73                                                           tol = 2 * np.pi * 1e-3 / (60.0 * 60.0 * 24.0)  # 1mas in radians
    74                                                           lat_close = np.isclose(telescope_obj.telescope_location_lat_lon_alt[0],
    75                                                                                  latitude, rtol=0, atol=tol)
    76                                                           lon_close = np.isclose(telescope_obj.telescope_location_lat_lon_alt[1],
    77                                                                                  longitude, rtol=0, atol=tol)
    78                                                           if correct_lat_lon:
    79                                                               self.telescope_location_lat_lon_alt = telescope_obj.telescope_location_lat_lon_alt
    80                                                           else:
    81                                                               self.telescope_location_lat_lon_alt = (latitude, longitude, telescope_obj.telescope_location_lat_lon_alt[2])
    82                                                           if lat_close and lon_close:
    83                                                               if correct_lat_lon:
    84                                                                   warnings.warn('Altitude is not present in Miriad file, '
    85                                                                                 'using known location values for '
    86                                                                                 '{telescope_name}.'.format(telescope_name=telescope_obj.telescope_name))
    87                                                               else:
    88                                                                   warnings.warn('Altitude is not present in Miriad file, '
    89                                                                                 'using known location altitude value '
    90                                                                                 'for {telescope_name} and lat/lon from '
    91                                                                                 'file.'.format(telescope_name=telescope_obj.telescope_name))
    92                                                           else:
    93                                                               warn_string = ('Altitude is not present in file ')
    94                                                               if not lat_close and not lon_close:
    95                                                                   warn_string = warn_string + 'and latitude and longitude values do not match values '
    96                                                               else:
    97                                                                   if not lat_close:
    98                                                                       warn_string = warn_string + 'and latitude value does not match value '
    99                                                                   else:
   100                                                                       warn_string = warn_string + 'and longitude value does not match value '
   101                                                               if correct_lat_lon:
   102                                                                   warn_string = (warn_string + 'for {telescope_name} in known '
   103                                                                                  'telescopes. Using values from known '
   104                                                                                  'telescopes.'.format(telescope_name=telescope_obj.telescope_name))
   105                                                                   warnings.warn(warn_string)
   106                                                               else:
   107                                                                   warn_string = (warn_string + 'for {telescope_name} in known '
   108                                                                                  'telescopes. Using altitude value from known '
   109                                                                                  'telescopes and lat/lon from '
   110                                                                                  'file.'.format(telescope_name=telescope_obj.telescope_name))
   111                                                                   warnings.warn(warn_string)
   112                                                       else:
   113                                                           warnings.warn('Altitude is not present in Miriad file, and ' +
   114                                                                         'telescope {telescope_name} is not in ' +
   115                                                                         'known_telescopes. Telescope location not '
   116                                                                         'set.'.format(telescope_name=self.telescope_name))
   117                                           
   118         1         1061   1061.0      0.0          self.history = uv['history']
   119         1           10     10.0      0.0          if self.pyuvdata_version_str not in self.history.replace('\n', ''):
   120         1            8      8.0      0.0              self.history += self.pyuvdata_version_str
   121         1            7      7.0      0.0          self.channel_width *= 1e9  # change from GHz to Hz
   122                                           
   123                                                   # read through the file and get the data
   124         1            8      8.0      0.0          _source = uv['source']  # check source of initial visibility
   125         1            5      5.0      0.0          data_accumulator = {}
   126         1            4      4.0      0.0          pol_list = []
   127    924673     18606339     20.1      6.5          for (uvw, t, (i, j)), d, f in uv.all(raw=True):
   128                                                       # control for the case of only a single spw not showing up in
   129                                                       # the dimension
   130                                                       # Note that the (i, j) tuple is calculated from a baseline number in
   131                                                       # aipy (see miriad_wrap.h). The i, j values are also adjusted by aipy
   132                                                       # to start at 0 rather than 1.
   133    924672      4337237      4.7      1.5              if len(d.shape) == 1:
   134    924672      5044929      5.5      1.8                  d.shape = (1,) + d.shape
   135    924672      6079467      6.6      2.1                  self.Nspws = d.shape[0]
   136    924672      8114680      8.8      2.8                  self.spw_array = np.arange(self.Nspws)
   137                                                       else:
   138                                                           raise(ValueError, """Sorry.  Files with more than one spectral
   139                                                                 window (spw) are not yet supported. A great
   140                                                                 project for the interested student!""")
   141    924672      3594801      3.9      1.3              try:
   142    924672      8281805      9.0      2.9                  cnt = uv['cnt']
   143                                                       except(KeyError):
   144                                                           cnt = np.ones(d.shape, dtype=np.float)
   145    924672      5535689      6.0      1.9              ra = uv['ra']
   146    924672      5230173      5.7      1.8              dec = uv['dec']
   147    924672      5241237      5.7      1.8              lst = uv['lst']
   148    924672      5383193      5.8      1.9              source = uv['source']
   149    924672      3729196      4.0      1.3              if source != _source:
   150                                                           raise(ValueError, 'This appears to be a multi source file, which is not supported.')
   151                                                       else:
   152    924672      3573790      3.9      1.2                  _source = source
   153                                           
   154    924672      3526011      3.8      1.2              try:
   155    924672      5529045      6.0      1.9                  data_accumulator[uv['pol']].append([uvw, t, i, j, d, f, cnt,
   156    924670      5690368      6.2      2.0                                                      ra, dec])
   157         2           10      5.0      0.0              except(KeyError):
   158         2            8      4.0      0.0                  data_accumulator[uv['pol']] = [[uvw, t, i, j, d, f, cnt,
   159         2           25     12.5      0.0                                                  ra, dec]]
   160         2           12      6.0      0.0                  pol_list.append(uv['pol'])
   161                                                           # NB: flag types in miriad are usually ints
   162                                           
   163         3       116059  38686.3      0.0          for pol,data in data_accumulator.iteritems():
   164         2       960005 480002.5      0.3              data_accumulator[pol] = np.array(data)
   165                                           
   166         1           44     44.0      0.0          self.polarization_array = np.array(pol_list)
   167         1           11     11.0      0.0          if len(self.polarization_array) != self.Npols:
   168                                                       warnings.warn('npols={npols} but found {l} pols in data file'.format(
   169                                                           npols=self.Npols, l=len(self.polarization_array)))
   170                                           
   171                                                   # makes a data_array (and flag_array) of zeroes to be filled in by
   172                                                   #   data values
   173                                                   # any missing data will have zeros
   174                                           
   175                                                   # use set to get the unique list of all times ever listed in the file
   176                                                   # iterate over polarizations and all spectra (bls and times) in two
   177                                                   # nested loops, then flatten into a single vector, then set
   178                                                   # then list again.
   179                                           
   180         1            4      4.0      0.0          times = list(set(
   181    924675      3855604      4.2      1.3              np.concatenate([[k[1] for k in d] for d in data_accumulator.values()])))
   182         1         4357   4357.0      0.0          times = np.sort(times)
   183                                           
   184         1            4      4.0      0.0          ant_i_unique = list(set(
   185    924675      3710282      4.0      1.3              np.concatenate([[k[2] for k in d] for d in data_accumulator.values()])))
   186         1            8      8.0      0.0          ant_j_unique = list(set(
   187    924675      3692915      4.0      1.3              np.concatenate([[k[3] for k in d] for d in data_accumulator.values()])))
   188                                           
   189         1           40     40.0      0.0          sorted_unique_ants = sorted(list(set(ant_i_unique + ant_j_unique)))
   190         1           29     29.0      0.0          ant_i_unique = np.array(ant_i_unique)
   191         1           16     16.0      0.0          ant_j_unique = np.array(ant_j_unique)
   192                                           
   193                                                   # Determine maximum digits needed to distinguish different values
   194         1           86     86.0      0.0          ndig_ant = np.ceil(np.log10(sorted_unique_ants[-1])).astype(int) + 1
   195                                                   # Be excessive in precision because we use the floating point values as dictionary keys later
   196         1           39     39.0      0.0          prec_t = - 2 * np.floor(np.log10(self._time_array.tols[-1])).astype(int)
   197         1           15     15.0      0.0          ndig_t = (np.ceil(np.log10(times[-1])).astype(int) + prec_t + 2)
   198         1            4      4.0      0.0          blts = []
   199         3           16      5.3      0.0          for d in data_accumulator.values():
   200    924674      3863501      4.2      1.3              for k in d:
   201    924672      6986079      7.6      2.4                  blt = ["{1:.{0}f}".format(prec_t, k[1]).zfill(ndig_t),
   202    924672      4985296      5.4      1.7                         str(k[2]).zfill(ndig_ant), str(k[3]).zfill(ndig_ant)]
   203    924672      3945795      4.3      1.4                  blt = "_".join(blt)
   204    924672      3696188      4.0      1.3                  blts.append(blt)
   205         1       346633 346633.0      0.1          reverse_inds= np.argsort(np.unique(np.array(blts),return_index=True)[1])
   206         1      1577533 1577533.0      0.6          unique_blts = np.sort(np.unique(np.array(blts)))
   207                                           
   208         1      1632683 1632683.0      0.6          reverse_inds = dict(zip(np.unique(np.array(blts)),reverse_inds))
   209         1           29     29.0      0.0          self.Nants_data = len(sorted_unique_ants)
   210                                           
   211                                                   # Miriad has no way to keep track of antenna numbers, so the antenna
   212                                                   # numbers are simply the index for each antenna in any array that
   213                                                   # describes antenna attributes (e.g. antpos for the antenna_postions).
   214                                                   # Therefore on write, nants (which gives the size of the antpos array)
   215                                                   # needs to be increased to be the max value of antenna_numbers+1 and the
   216                                                   # antpos array needs to be inflated with zeros at locations where we
   217                                                   # don't have antenna information. These inflations need to be undone at
   218                                                   # read. If the file was written by pyuvdata, then the variable antnums
   219                                                   # will be present and we can use it, otherwise we need to test for zeros
   220                                                   # in the antpos array and/or antennas with no visibilities.
   221         1            5      5.0      0.0          try:
   222                                                       # The antnums variable will only exist if the file was written with pyuvdata.
   223                                                       # For some reason Miriad doesn't handle an array of integers properly,
   224                                                       # so we convert to floats on write and back here
   225         1           43     43.0      0.0              self.antenna_numbers = uv['antnums'].astype(int)
   226         1           10     10.0      0.0              self.Nants_telescope = len(self.antenna_numbers)
   227                                                   except(KeyError):
   228                                                       self.antenna_numbers = None
   229                                                       self.Nants_telescope = None
   230                                           
   231         1            7      7.0      0.0          nants = uv['nants']
   232         1            4      4.0      0.0          try:
   233         1           24     24.0      0.0              antpos = uv['antpos'].reshape(3, nants).T
   234                                                       if self.Nants_telescope is not None:
   235                                                           # in this case there is an antnums variable
   236                                                           # (meaning that the file was written with pyuvdata), so we'll use it
   237                                                           if nants == self.Nants_telescope:
   238                                                               # no inflation, so just use the positions
   239                                                               self.antenna_positions = antpos
   240                                                           else:
   241                                                               # there is some inflation, just use the ones that appear in antnums
   242                                                               self.antenna_positions = np.zeros((self.Nants_telescope, 3), dtype=antpos.dtype)
   243                                                               for ai, num in enumerate(self.antenna_numbers):
   244                                                                   self.antenna_positions[ai, :] = antpos[num, :]
   245                                                       else:
   246                                                           # there is no antnums variable (meaning that this file was not
   247                                                           # written by pyuvdata), so we test for antennas with non-zero
   248                                                           # positions and/or that appear in the visibility data
   249                                                           # (meaning that they have entries in ant_1_array or ant_2_array)
   250                                                           antpos_length = np.sqrt(np.sum(np.abs(antpos)**2, axis=1))
   251                                                           good_antpos = np.where(antpos_length > 0)[0]
   252                                                           # take the union of the antennas with good positions (good_antpos)
   253                                                           # and the antennas that have visisbilities (sorted_unique_ants)
   254                                                           # if there are antennas with visibilities but zeroed positions we issue a warning below
   255                                                           ants_use = set(good_antpos).union(sorted_unique_ants)
   256                                                           # ants_use are the antennas we'll keep track of in the UVData
   257                                                           # object, so they dictate Nants_telescope
   258                                                           self.Nants_telescope = len(ants_use)
   259                                                           self.antenna_numbers = np.array(list(ants_use))
   260                                                           self.antenna_positions = np.zeros((self.Nants_telescope, 3), dtype=antpos.dtype)
   261                                                           for ai, num in enumerate(self.antenna_numbers):
   262                                                               self.antenna_positions[ai, :] = antpos[num, :]
   263                                                               if antpos_length[num] == 0:
   264                                                                   warnings.warn('antenna number {n} has visibilities '
   265                                                                                 'associated with it, but it has a position'
   266                                                                                 ' of (0,0,0)'.format(n=num))
   267         1            4      4.0      0.0          except(KeyError):
   268                                                       # there is no antpos variable
   269         1            7      7.0      0.0              self.antenna_positions = None
   270                                           
   271         1            5      5.0      0.0          if self.antenna_numbers is None:
   272                                                       # there are no antenna_numbers or antenna_positions, so just use
   273                                                       # the antennas present in the visibilities
   274                                                       # (Nants_data will therefore match Nants_telescope)
   275                                                       self.antenna_numbers = np.array(sorted_unique_ants)
   276                                                       self.Nants_telescope = len(self.antenna_numbers)
   277                                           
   278                                                   # antenna names is a foreign concept in miriad but required in other formats.
   279         1            4      4.0      0.0          try:
   280                                                       # Here we deal with the way pyuvdata tacks it on to keep the name information if we have it:
   281                                                       # This is a horrible hack to save & recover antenna_names array. Miriad can't handle arrays
   282                                                       # of strings or a long enough single string to put them all into one string
   283                                                       # so we convert them into hex values and then into floats on write and convert
   284                                                       # back to strings here
   285         1            9      9.0      0.0              ant_name_flt = uv['antnames']
   286       129         1309     10.1      0.0              ant_name_list = [('%x' % elem.astype(np.int64)).decode('hex') for elem in ant_name_flt]
   287         1            8      8.0      0.0              self.antenna_names = ant_name_list
   288                                                   except(KeyError):
   289                                                       self.antenna_names = self.antenna_numbers.astype(str).tolist()
   290                                           
   291                                                   # form up a grid which indexes time and baselines along the 'long'
   292                                                   # axis of the visdata array
   293                                           
   294         1      2008998 2008998.0      0.7          tij_grid = np.array(map(lambda x: map(float, x.split("_")), unique_blts))
   295         1           39     39.0      0.0          t_grid, ant_i_grid, ant_j_grid = tij_grid.T
   296                                                   # set the data sizes
   297         1            4      4.0      0.0          try:
   298         1           37     37.0      0.0              self.Nblts = uv['nblts']
   299         1           10     10.0      0.0              if self.Nblts != len(t_grid):
   300                                                           warnings.warn('Nblts does not match the number of unique blts in the data')
   301                                                   except(KeyError):
   302                                                       self.Nblts = len(t_grid)
   303         1            4      4.0      0.0          try:
   304         1           12     12.0      0.0              self.Ntimes = uv['ntimes']
   305         1            6      6.0      0.0              if self.Ntimes != len(times):
   306                                                           warnings.warn('Ntimes does not match the number of unique times in the data')
   307                                                   except(KeyError):
   308                                                       self.Ntimes = len(times)
   309                                           
   310         1            6      6.0      0.0          self.time_array = t_grid
   311         1         1172   1172.0      0.0          self.ant_1_array = ant_i_grid.astype(int)
   312         1         2441   2441.0      0.0          self.ant_2_array = ant_j_grid.astype(int)
   313                                           
   314         1         3470   3470.0      0.0          self.baseline_array = self.antnums_to_baseline(ant_i_grid.astype(int),
   315         1        16148  16148.0      0.0                                                         ant_j_grid.astype(int))
   316         1            5      5.0      0.0          try:
   317         1           18     18.0      0.0              self.Nbls = uv['nbls']
   318         1        21750  21750.0      0.0              if self.Nbls != len(np.unique(self.baseline_array)):
   319                                                           warnings.warn('Nbls does not match the number of unique baselines in the data')
   320                                                   except(KeyError):
   321                                                       self.Nbls = len(np.unique(self.baseline_array))
   322                                           
   323                                                   # slot the data into a grid
   324         1           12     12.0      0.0          self.data_array = np.zeros((self.Nblts, self.Nspws, self.Nfreqs,
   325         1           31     31.0      0.0                                      self.Npols), dtype=np.complex64)
   326         1       107702 107702.0      0.0          self.flag_array = np.ones(self.data_array.shape, dtype=np.bool)
   327         1         1679   1679.0      0.0          self.uvw_array = np.zeros((self.Nblts, 3))
   328                                                   # NOTE: Using our lst calculator, which uses astropy,
   329                                                   # instead of aipy values which come from pyephem.
   330                                                   # The differences are of order 5 seconds.
   331         1      3212863 3212863.0      1.1          self.set_lsts_from_time_array()
   332         1       660518 660518.0      0.2          self.nsample_array = np.ones(self.data_array.shape, dtype=np.float)
   333         1           61     61.0      0.0          self.freq_array = (np.arange(self.Nfreqs) * self.channel_width +
   334         1           31     31.0      0.0                             uv['sfreq'] * 1e9)
   335                                                   # Tile freq_array to shape (Nspws, Nfreqs).
   336                                                   # Currently does not actually support Nspws>1!
   337         1           60     60.0      0.0          self.freq_array = np.tile(self.freq_array, (self.Nspws, 1))
   338                                           
   339                                                   # Temporary arrays to hold polarization axis, which will be collapsed
   340         1         1289   1289.0      0.0          ra_pol_list = np.zeros((self.Nblts, self.Npols))
   341         1          800    800.0      0.0          dec_pol_list = np.zeros((self.Nblts, self.Npols))
   342         1          453    453.0      0.0          uvw_pol_list = np.zeros((self.Nblts, 3, self.Npols))
   343         1         1562   1562.0      0.0          c_ns = const.c.to('m/ns').value
   344         3       119371  39790.3      0.0          for pol, data in data_accumulator.iteritems():
   345         2          217    108.5      0.0              pol_ind = self._pol_to_ind(pol)
   346    924674      4384762      4.7      1.5              for ind, d in enumerate(data):
   347    924672      9734240     10.5      3.4                  blt = ["{1:.{0}f}".format(prec_t, d[1]).zfill(ndig_t),
   348    924672      5654097      6.1      2.0                         str(d[2]).zfill(ndig_ant), str(d[3]).zfill(ndig_ant)]
   349    924672      4244202      4.6      1.5                  blt = "_".join(blt)
   350    924672      4044832      4.4      1.4                  blt_index = reverse_inds[blt]
   351                                           
   352    924672     12337336     13.3      4.3                  self.data_array[blt_index, :, :, pol_ind] = d[4]
   353    924672      9781916     10.6      3.4                  self.flag_array[blt_index, :, :, pol_ind] = d[5]
   354    924672      9661016     10.4      3.4                  self.nsample_array[blt_index, :, :, pol_ind] = d[6]
   355                                           
   356                                                           # because there are uvws/ra/dec for each pol, and one pol may not
   357                                                           # have that visibility, we collapse along the polarization
   358                                                           # axis but avoid any missing visbilities
   359    924672      8283971      9.0      2.9                  uvw = d[0] * c_ns
   360    924672      5304313      5.7      1.9                  uvw.shape = (1, 3)
   361    924672      7964836      8.6      2.8                  uvw_pol_list[blt_index, :, pol_ind] = uvw
   362    924672      6968192      7.5      2.4                  ra_pol_list[blt_index, pol_ind] = d[7]
   363    924672      6485602      7.0      2.3                  dec_pol_list[blt_index, pol_ind] = d[8]
   364                                           
   365                                                   # Collapse pol axis for ra_list, dec_list, and uvw_list
   366         1          655    655.0      0.0          ra_list = np.zeros(self.Nblts)
   367         1          648    648.0      0.0          dec_list = np.zeros(self.Nblts)
   368    462337      1846826      4.0      0.6          for blt_index in xrange(self.Nblts):
   369    462336     11539119     25.0      4.0              test = ~np.all(self.flag_array[blt_index, :, :, :], axis=(0, 1))
   370    462336      2873731      6.2      1.0              good_pol = np.where(test)[0]
   371    462336      2003156      4.3      0.7              if len(good_pol) == 1:
   372                                                           # Only one good pol, use it
   373                                                           self.uvw_array[blt_index, :] = uvw_pol_list[blt_index, :, good_pol]
   374                                                           ra_list[blt_index] = ra_pol_list[blt_index, good_pol]
   375                                                           dec_list[blt_index] = dec_pol_list[blt_index, good_pol]
   376    462336      1887218      4.1      0.7              elif len(good_pol) > 1:
   377                                                           # Multiple good pols, check for consistency. pyuvdata does not
   378                                                           # support pol-dependent uvw, ra, or dec.
   379    216991      6986427     32.2      2.4                  if np.any(np.diff(uvw_pol_list[blt_index, :, good_pol], axis=0)):
   380                                                               raise ValueError('uvw values are different by polarization.')
   381                                                           else:
   382    216991      2026127      9.3      0.7                      self.uvw_array[blt_index, :] = uvw_pol_list[blt_index, :, good_pol[0]]
   383    216991      6339086     29.2      2.2                  if np.any(np.diff(ra_pol_list[blt_index, good_pol])):
   384                                                               raise ValueError('ra values are different by polarization.')
   385                                                           else:
   386    216991      1191884      5.5      0.4                      ra_list[blt_index] = ra_pol_list[blt_index, good_pol[0]]
   387    216991      6089454     28.1      2.1                  if np.any(np.diff(dec_pol_list[blt_index, good_pol])):
   388                                                               raise ValueError('dec values are different by polarization.')
   389                                                           else:
   390    216991      1162411      5.4      0.4                      dec_list[blt_index] = dec_pol_list[blt_index, good_pol[0]]
   391                                                       else:
   392                                                           # No good pols for this blt. Fill with first one.
   393    245345      1995663      8.1      0.7                  self.uvw_array[blt_index, :] = uvw_pol_list[blt_index, :, 0]
   394    245345      1160717      4.7      0.4                  ra_list[blt_index] = ra_pol_list[blt_index, 0]
   395    245345      1047888      4.3      0.4                  dec_list[blt_index] = dec_pol_list[blt_index, 0]
   396                                           
   397                                                   # check if ra is constant throughout file; if it is,
   398                                                   # file is tracking if not, file is drift scanning
   399         1        34598  34598.0      0.0          blt_good = np.where(~np.all(self.flag_array, axis=(1, 2, 3)))
   400         1         2293   2293.0      0.0          if np.isclose(np.mean(np.diff(ra_list[blt_good])), 0.):
   401                                                       self.set_phased()
   402                                                       self.phase_center_ra = ra_list[0]
   403                                                       self.phase_center_dec = dec_list[0]
   404                                                       self.phase_center_epoch = uv['epoch']
   405                                                   else:
   406         1           21     21.0      0.0              self.set_drift()
   407         1            7      7.0      0.0              self.zenith_ra = ra_list
   408         1            6      6.0      0.0              self.zenith_dec = dec_list
   409         1            6      6.0      0.0          self.vis_units = 'UNCALIB'  # assume no calibration
   410                                           
   411         1            4      4.0      0.0          try:
   412         1         1022   1022.0      0.0              self.set_telescope_params()
   413                                                   except ValueError, ve:
   414                                                       warnings.warn(str(ve))
   415                                           
   416                                                   # check if object has all required uv_properties set
   417         1            4      4.0      0.0          if run_check:
   418         1       127303 127303.0      0.0              self.check(run_check_acceptability=run_check_acceptability)

Total time: 76.8959 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/miriad.py
Function: write_miriad at line 419

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   419                                               @profile
   420                                               def write_miriad(self, filepath, run_check=True, run_check_acceptability=True,
   421                                                                clobber=False, no_antnums=False):
   422                                                   """
   423                                                   Write the data to a miriad file.
   424                                           
   425                                                   Args:
   426                                                       filename: The miriad file directory to write to.
   427                                                       run_check: Option to check for the existence and proper shapes of
   428                                                           required parameters before writing the file. Default is True.
   429                                                       run_check_acceptability: Option to check acceptable range of the values of
   430                                                           required parameters before writing the file. Default is True.
   431                                                       clobber: Option to overwrite the filename if the file already exists.
   432                                                           Default is False.
   433                                                       no_antnums: Option to not write the antnums variable to the file.
   434                                                           Should only be used for testing purposes.
   435                                                   """
   436                                                   # check for multiple spws
   437         1            8      8.0      0.0          if self.data_array.shape[1] > 1:
   438                                                       raise ValueError('write_miriad currently only handles single spw files.')
   439                                           
   440         1           85     85.0      0.0          if os.path.exists(filepath):
   441                                                       if clobber:
   442                                                           print 'File exists: clobbering'
   443                                                           shutil.rmtree(filepath)
   444                                                       else:
   445                                                           raise ValueError('File exists: skipping')
   446                                           
   447         1           16     16.0      0.0          freq_spacing = self.freq_array[0, 1:] - self.freq_array[0, :-1]
   448         1           45     45.0      0.0          if not np.isclose(np.min(freq_spacing), np.max(freq_spacing),
   449         1          114    114.0      0.0                            rtol=self._freq_array.tols[0], atol=self._freq_array.tols[1]):
   450                                                       raise ValueError('The frequencies are not evenly spaced (probably '
   451                                                                        'because of a select operation). The miriad format '
   452                                                                        'does not support unevenly spaced frequencies.')
   453         1           16     16.0      0.0          if not np.isclose(np.max(freq_spacing), self.channel_width,
   454         1           81     81.0      0.0                            rtol=self._freq_array.tols[0], atol=self._freq_array.tols[1]):
   455                                                       raise ValueError('The frequencies are separated by more than their '
   456                                                                        'channel width (probably because of a select operation). '
   457                                                                        'The miriad format does not support frequencies '
   458                                                                        'that are spaced by more than their channel width.')
   459                                           
   460         1        20189  20189.0      0.0          uv = a.miriad.UV(filepath, status='new')
   461                                           
   462                                                   # initialize header variables
   463         1           31     31.0      0.0          uv._wrhd('obstype', 'mixed-auto-cross')
   464         1         1498   1498.0      0.0          uv._wrhd('history', self.history + '\n')
   465                                           
   466                                                   # recognized miriad variables
   467         1            8      8.0      0.0          uv.add_var('nchan', 'i')
   468         1           13     13.0      0.0          uv['nchan'] = self.Nfreqs
   469         1            3      3.0      0.0          uv.add_var('npol', 'i')
   470         1            6      6.0      0.0          uv['npol'] = self.Npols
   471         1            2      2.0      0.0          uv.add_var('nspect', 'i')
   472         1            6      6.0      0.0          uv['nspect'] = self.Nspws
   473         1            3      3.0      0.0          uv.add_var('inttime', 'd')
   474         1         1219   1219.0      0.0          uv['inttime'] = self.integration_time
   475         1           17     17.0      0.0          uv.add_var('sdf', 'd')
   476         1           11     11.0      0.0          uv['sdf'] = self.channel_width / 1e9  # in GHz
   477         1            3      3.0      0.0          uv.add_var('source', 'a')
   478         1            6      6.0      0.0          uv['source'] = self.object_name
   479         1            3      3.0      0.0          uv.add_var('telescop', 'a')
   480         1            6      6.0      0.0          uv['telescop'] = self.telescope_name
   481         1            3      3.0      0.0          uv.add_var('latitud', 'd')
   482         1           89     89.0      0.0          uv['latitud'] = self.telescope_location_lat_lon_alt[0]
   483         1            3      3.0      0.0          uv.add_var('longitu', 'd')
   484         1           71     71.0      0.0          uv['longitu'] = self.telescope_location_lat_lon_alt[1]
   485         1            3      3.0      0.0          uv.add_var('nants', 'i')
   486                                           
   487                                                   # Miriad has no way to keep track of antenna numbers, so the antenna
   488                                                   # numbers are simply the index for each antenna in any array that
   489                                                   # describes antenna attributes (e.g. antpos for the antenna_postions).
   490                                                   # Therefore on write, nants (which gives the size of the antpos array)
   491                                                   # needs to be increased to be the max value of antenna_numbers+1 and the
   492                                                   # antpos array needs to be inflated with zeros at locations where we
   493                                                   # don't have antenna information. These inflations need to be undone at
   494                                                   # read. If the file was written by pyuvdata, then the variable antnums
   495                                                   # will be present and we can use it, otherwise we need to test for zeros
   496                                                   # in the antpos array and/or antennas with no visibilities.
   497         1           22     22.0      0.0          nants = np.max(self.antenna_numbers) + 1
   498         1            6      6.0      0.0          uv['nants'] = nants
   499         1            4      4.0      0.0          if self.antenna_positions is not None:
   500                                                       antpos = np.zeros((nants, 3), dtype=self.antenna_positions.dtype)
   501                                                       for ai, num in enumerate(self.antenna_numbers):
   502                                                           antpos[num, :] = self.antenna_positions[ai, :]
   503                                                       uv.add_var('antpos', 'd')
   504                                                       uv['antpos'] = antpos.T.flatten()
   505                                           
   506         1            3      3.0      0.0          uv.add_var('sfreq', 'd')
   507         1            7      7.0      0.0          uv['sfreq'] = self.freq_array[0, 0] / 1e9  # first spw; in GHz
   508         1            3      3.0      0.0          if self.phase_type == 'phased':
   509                                                       uv.add_var('epoch', 'r')
   510                                                       uv['epoch'] = self.phase_center_epoch
   511                                           
   512                                                   # required pyuvdata variables that are not recognized miriad variables
   513         1            3      3.0      0.0          uv.add_var('ntimes', 'i')
   514         1            5      5.0      0.0          uv['ntimes'] = self.Ntimes
   515         1            3      3.0      0.0          uv.add_var('nbls', 'i')
   516         1            5      5.0      0.0          uv['nbls'] = self.Nbls
   517         1            3      3.0      0.0          uv.add_var('nblts', 'i')
   518         1            6      6.0      0.0          uv['nblts'] = self.Nblts
   519         1            3      3.0      0.0          uv.add_var('visunits', 'a')
   520         1            6      6.0      0.0          uv['visunits'] = self.vis_units
   521         1            3      3.0      0.0          uv.add_var('instrume', 'a')
   522         1            5      5.0      0.0          uv['instrume'] = self.instrument
   523         1            3      3.0      0.0          uv.add_var('altitude', 'd')
   524         1           70     70.0      0.0          uv['altitude'] = self.telescope_location_lat_lon_alt[2]
   525                                           
   526         1            2      2.0      0.0          if not no_antnums:
   527                                                       # Add in the antenna_numbers so we have them if we read this file back in
   528                                                       # for some reason Miriad doesn't handle an array of integers properly,
   529                                                       # so convert to floats here and integers on read
   530         1            3      3.0      0.0              uv.add_var('antnums', 'd')
   531         1           12     12.0      0.0              uv['antnums'] = self.antenna_numbers.astype(np.float64)
   532                                           
   533                                                   # antenna names is a foreign concept in miriad but required in other formats.
   534                                                   # This is a horrible hack to save antenna_names array. Miriad can't handle arrays
   535                                                   # of strings or a long enough single string to put them all into one string
   536                                                   # so we convert them into hex values and then into floats and convert
   537                                                   # back to strings on read
   538       129          622      4.8      0.0          ant_name_flt = np.array([int(elem.encode("hex"), 16) for elem in self.antenna_names]).astype(np.float64)
   539         1            3      3.0      0.0          uv.add_var('antnames', 'd')
   540         1            8      8.0      0.0          uv['antnames'] = ant_name_flt
   541                                           
   542                                                   # variables that can get updated with every visibility
   543         1            3      3.0      0.0          uv.add_var('pol', 'i')
   544         1            4      4.0      0.0          uv.add_var('lst', 'd')
   545         1            3      3.0      0.0          uv.add_var('cnt', 'd')
   546         1            2      2.0      0.0          uv.add_var('ra', 'd')
   547         1            3      3.0      0.0          uv.add_var('dec', 'd')
   548                                           
   549                                                   # write data
   550         1          858    858.0      0.0          c_ns = const.c.to('m/ns').value
   551    462337      1317473      2.8      1.7          for viscnt, blt in enumerate(self.data_array):
   552    462336      5310432     11.5      6.9              uvw = (self.uvw_array[viscnt] / c_ns).astype(np.double)  # NOTE issue 50 on conjugation
   553    462336      1944312      4.2      2.5              t = self.time_array[viscnt]
   554    462336      1568615      3.4      2.0              i = self.ant_1_array[viscnt]
   555    462336      1490438      3.2      1.9              j = self.ant_2_array[viscnt]
   556                                           
   557    462336      2884012      6.2      3.8              uv['lst'] = self.lst_array[viscnt]
   558    462336      1600097      3.5      2.1              if self.phase_type == 'phased':
   559                                                           uv['ra'] = self.phase_center_ra
   560                                                           uv['dec'] = self.phase_center_dec
   561    462336      1443938      3.1      1.9              elif self.phase_type == 'drift':
   562    462336      2515366      5.4      3.3                  uv['ra'] = self.zenith_ra[viscnt]
   563    462336      2590451      5.6      3.4                  uv['dec'] = self.zenith_dec[viscnt]
   564                                                       else:
   565                                                           raise ValueError('The phasing type of the data is unknown. '
   566                                                                            'Set the phase_type to "drift" or "phased" to '
   567                                                                            'reflect the phasing status of the data')
   568                                           
   569                                                       # NOTE only writing spw 0, not supporting multiple spws for write
   570   1387008      6176944      4.5      8.0              for polcnt, pol in enumerate(self.polarization_array):
   571    924672      7376416      8.0      9.6                  uv['pol'] = pol.astype(np.int)
   572    924672     11173991     12.1     14.5                  uv['cnt'] = self.nsample_array[viscnt, 0, :, polcnt].astype(np.double)
   573                                           
   574    924672      4099620      4.4      5.3                  data = self.data_array[viscnt, 0, :, polcnt]
   575    924672      3612908      3.9      4.7                  flags = self.flag_array[viscnt, 0, :, polcnt]
   576    924672      2280073      2.5      3.0                  if i > j:
   577                                                               i, j, data = j, i, np.conjugate(data)
   578    924672      2187587      2.4      2.8                  preamble = (uvw, t, (i, j))
   579                                           
   580    924672     17172875     18.6     22.3                  uv.write(preamble, data, flags)
   581         1            2      2.0      0.0          if run_check:
   582                                                       """Check for acceptable units."""
   583         1       125145 125145.0      0.2              self.check(run_check_acceptability=run_check_acceptability)

Total time: 3.82707 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/uvdata.py
Function: set_lsts_from_time_array at line 457

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   457                                               @profile
   458                                               def set_lsts_from_time_array(self):
   459                                                   """Set the lst_array based from the time_array."""
   460         2            3      1.5      0.0          lsts = []
   461         2            9      4.5      0.0          curtime = self.time_array[0]
   462         2         1093    546.5      0.0          self.lst_array = np.zeros(self.Nblts)
   463         2          213    106.5      0.0          latitude, longitude, altitude = self.telescope_location_lat_lon_alt_degrees
   464       114        18347    160.9      0.5          for ind, jd in enumerate(np.unique(self.time_array)):
   465       112       151277   1350.7      4.0              t = Time(jd, format='jd', location=(longitude, latitude))
   466       112      3656128  32644.0     95.5              self.lst_array[np.where(np.isclose(jd, self.time_array, atol=1e-6, rtol=1e-12))] = t.sidereal_time('apparent').radian

Total time: 0.00046 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/uvdata.py
Function: juldate2ephem at line 468

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   468                                               @profile
   469                                               def juldate2ephem(self, num):
   470                                                   """
   471                                                   Convert Julian date to ephem date, measured from noon, Dec. 31, 1899.
   472                                           
   473                                                   Args:
   474                                                       num: Julian date
   475                                           
   476                                                   Returns:
   477                                                       ephem date, measured from noon, Dec. 31, 1899.
   478                                                   """
   479       227          460      2.0    100.0          return ephem.date(num - 2415020.)

Total time: 16.0526 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/uvdata.py
Function: unphase_to_drift at line 481

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   481                                               @profile
   482                                               def unphase_to_drift(self):
   483                                                   """Convert from a phased dataset to a drift dataset."""
   484         1            4      4.0      0.0          if self.phase_type == 'phased':
   485         1            2      2.0      0.0              pass
   486                                                   elif self.phase_type == 'drift':
   487                                                       raise ValueError('The data is already drift scanning; can only ' +
   488                                                                        'unphase phased data.')
   489                                                   else:
   490                                                       raise ValueError('The phasing type of the data is unknown. '
   491                                                                        'Set the phase_type to drift or phased to '
   492                                                                        'reflect the phasing status of the data')
   493                                           
   494         1           69     69.0      0.0          latitude, longitude, altitude = self.telescope_location_lat_lon_alt
   495                                           
   496         1            5      5.0      0.0          obs = ephem.Observer()
   497                                                   # obs inits with default values for parameters -- be sure to replace them
   498         1            4      4.0      0.0          obs.lat = latitude
   499         1            1      1.0      0.0          obs.lon = longitude
   500                                           
   501         1            3      3.0      0.0          phase_center = ephem.FixedBody()
   502         1            5      5.0      0.0          epoch = (self.phase_center_epoch - 2000.) * 365.2422 + ephem.J2000  # convert years to ephemtime
   503         1            2      2.0      0.0          phase_center._epoch = epoch
   504         1            3      3.0      0.0          phase_center._ra = self.phase_center_ra
   505         1            3      3.0      0.0          phase_center._dec = self.phase_center_dec
   506                                           
   507         1          190    190.0      0.0          self.zenith_ra = np.zeros_like(self.time_array)
   508         1          332    332.0      0.0          self.zenith_dec = np.zeros_like(self.time_array)
   509                                           
   510                                                   # apply -w phasor
   511                                                   w_lambda = (self.uvw_array[:, 2].reshape(self.Nblts, 1) /
   512         1       570278 570278.0      3.6                      const.c.to('m/s').value * self.freq_array.reshape(1, self.Nfreqs))
   513         1     11725581 11725581.0     73.0          phs = np.exp(-1j * 2 * np.pi * (-1) * w_lambda[:, None, :, None])
   514         1      3695210 3695210.0     23.0          self.data_array *= phs
   515                                           
   516         1         8690   8690.0      0.1          unique_times = np.unique(self.time_array)
   517        57          144      2.5      0.0          for jd in unique_times:
   518        56        17743    316.8      0.1              inds = np.where(self.time_array == jd)[0]
   519        56          829     14.8      0.0              obs.date, obs.epoch = self.juldate2ephem(jd), self.juldate2ephem(jd)
   520        56          120      2.1      0.0              phase_center.compute(obs)
   521        56         2698     48.2      0.0              phase_center_ra, phase_center_dec = phase_center.a_ra, phase_center.a_dec
   522        56          106      1.9      0.0              zenith_ra = obs.sidereal_time()
   523        56           79      1.4      0.0              zenith_dec = latitude
   524        56         1388     24.8      0.0              self.zenith_ra[inds] = zenith_ra
   525        56         1278     22.8      0.0              self.zenith_dec[inds] = zenith_dec
   526                                           
   527                                                       # generate rotation matrices
   528        56         1780     31.8      0.0              m0 = uvutils.top2eq_m(0., phase_center_dec)
   529        56         1357     24.2      0.0              m1 = uvutils.eq2top_m(phase_center_ra - zenith_ra, zenith_dec)
   530                                           
   531                                                       # rotate and write uvws
   532        56         8072    144.1      0.1              uvw = self.uvw_array[inds, :]
   533        56         4558     81.4      0.0              uvw = np.dot(m0, uvw.T).T
   534        56         3378     60.3      0.0              uvw = np.dot(m1, uvw.T).T
   535        56         8619    153.9      0.1              self.uvw_array[inds, :] = uvw
   536                                           
   537                                                   # remove phase center
   538         1            6      6.0      0.0          self.phase_center_ra = None
   539         1            4      4.0      0.0          self.phase_center_dec = None
   540         1            3      3.0      0.0          self.phase_center_epoch = None
   541         1            9      9.0      0.0          self.set_drift()

Total time: 15.9646 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/uvdata.py
Function: phase_to_time at line 543

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   543                                               @profile
   544                                               def phase_to_time(self, time):
   545                                                   """
   546                                                   Phase a drift scan dataset to the ra/dec of zenith at a particular time.
   547                                           
   548                                                   Args:
   549                                                       time: The time to phase to.
   550                                                   """
   551         1            2      2.0      0.0          if self.phase_type == 'drift':
   552         1            1      1.0      0.0              pass
   553                                                   elif self.phase_type == 'phased':
   554                                                       raise ValueError('The data is already phased; can only phase ' +
   555                                                                        'drift scanning data.')
   556                                                   else:
   557                                                       raise ValueError('The phasing type of the data is unknown. '
   558                                                                        'Set the phase_type to drift or phased to '
   559                                                                        'reflect the phasing status of the data')
   560                                           
   561         1           60     60.0      0.0          obs = ephem.Observer()
   562                                                   # obs inits with default values for parameters -- be sure to replace them
   563         1          119    119.0      0.0          latitude, longitude, altitude = self.telescope_location_lat_lon_alt
   564         1            4      4.0      0.0          obs.lat = latitude
   565         1            1      1.0      0.0          obs.lon = longitude
   566                                           
   567         1           35     35.0      0.0          obs.date, obs.epoch = self.juldate2ephem(time), self.juldate2ephem(time)
   568                                           
   569         1           85     85.0      0.0          ra = obs.sidereal_time()
   570         1            1      1.0      0.0          dec = latitude
   571         1            7      7.0      0.0          epoch = self.juldate2ephem(time)
   572         1     15964312 15964312.0    100.0          self.phase(ra, dec, epoch)

Total time: 15.9553 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/uvdata.py
Function: phase at line 574

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   574                                               @profile
   575                                               def phase(self, ra, dec, epoch):
   576                                                   """
   577                                                   Phase a drift scan dataset to a single ra/dec at a particular epoch.
   578                                           
   579                                                   Will not phase already phased data.
   580                                           
   581                                                   Args:
   582                                                       ra: The ra to phase to in radians.
   583                                                       dec: The dec to phase to in radians.
   584                                                       epoch: The epoch to use for phasing. Should be an ephem date,
   585                                                           measured from noon Dec. 31, 1899.
   586                                                   """
   587         1            3      3.0      0.0          if self.phase_type == 'drift':
   588         1            2      2.0      0.0              pass
   589                                                   elif self.phase_type == 'phased':
   590                                                       raise ValueError('The data is already phased; can only phase '
   591                                                                        'drift scan data. Use unphase_to_drift to '
   592                                                                        'convert to a drift scan.')
   593                                                   else:
   594                                                       raise ValueError('The phasing type of the data is unknown. '
   595                                                                        'Set the phase_type to "drift" or "phased" to '
   596                                                                        'reflect the phasing status of the data')
   597                                           
   598         1            2      2.0      0.0          obs = ephem.Observer()
   599                                                   # obs inits with default values for parameters -- be sure to replace them
   600         1           68     68.0      0.0          latitude, longitude, altitude = self.telescope_location_lat_lon_alt
   601         1            2      2.0      0.0          obs.lat = latitude
   602         1            1      1.0      0.0          obs.lon = longitude
   603                                           
   604                                                   # create a pyephem object for the phasing position
   605         1            2      2.0      0.0          precess_pos = ephem.FixedBody()
   606         1            2      2.0      0.0          precess_pos._ra = ra
   607         1            1      1.0      0.0          precess_pos._dec = dec
   608         1            1      1.0      0.0          precess_pos._epoch = epoch
   609                                           
   610                                                   # calculate RA/DEC in J2000 and write to object
   611         1            2      2.0      0.0          obs.date, obs.epoch = ephem.J2000, ephem.J2000
   612         1            3      3.0      0.0          precess_pos.compute(obs)
   613                                           
   614         1          142    142.0      0.0          self.phase_center_ra = precess_pos.a_ra + 0.0  # force to be a float not ephem.Angle
   615         1            5      5.0      0.0          self.phase_center_dec = precess_pos.a_dec + 0.0  # force to be a float not ephem.Angle
   616                                                   # explicitly set epoch to J2000
   617         1            4      4.0      0.0          self.phase_center_epoch = 2000.0
   618                                           
   619         1        11379  11379.0      0.1          unique_times, unique_inds = np.unique(self.time_array, return_index=True)
   620        57          157      2.8      0.0          for ind, jd in enumerate(unique_times):
   621        56        45881    819.3      0.3              inds = np.where(self.time_array == jd)[0]
   622        56          308      5.5      0.0              lst = self.lst_array[unique_inds[ind]]
   623                                                       # calculate ra/dec of phase center in current epoch
   624        56          784     14.0      0.0              obs.date, obs.epoch = self.juldate2ephem(jd), self.juldate2ephem(jd)
   625        56          117      2.1      0.0              precess_pos.compute(obs)
   626        56         2746     49.0      0.0              ra, dec = precess_pos.a_ra, precess_pos.a_dec
   627                                           
   628                                                       # generate rotation matrices
   629        56         1842     32.9      0.0              m0 = uvutils.top2eq_m(lst - obs.sidereal_time(), latitude)
   630        56         1439     25.7      0.0              m1 = uvutils.eq2top_m(lst - ra, dec)
   631                                           
   632                                                       # rotate and write uvws
   633        56        10616    189.6      0.1              uvw = self.uvw_array[inds, :]
   634        56         3756     67.1      0.0              uvw = np.dot(m0, uvw.T).T
   635        56         3291     58.8      0.0              uvw = np.dot(m1, uvw.T).T
   636        56         8710    155.5      0.1              self.uvw_array[inds, :] = uvw
   637                                           
   638                                                   # calculate data and apply phasor
   639                                                   w_lambda = (self.uvw_array[:, 2].reshape(self.Nblts, 1) /
   640         1       491238 491238.0      3.1                      const.c.to('m/s').value * self.freq_array.reshape(1, self.Nfreqs))
   641         1     11694982 11694982.0     73.3          phs = np.exp(-1j * 2 * np.pi * w_lambda[:, None, :, None])
   642         1      3677817 3677817.0     23.1          self.data_array *= phs
   643                                           
   644         1            4      4.0      0.0          del(obs)
   645         1           14     14.0      0.0          self.set_phased()

Total time: 14.4262 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/uvfits.py
Function: read_uvfits at line 45

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    45                                               @profile
    46                                               def read_uvfits(self, filename, run_check=True, run_check_acceptability=True):
    47                                                   """
    48                                                   Read in data from a uvfits file.
    49                                           
    50                                                   Args:
    51                                                       filename: The uvfits file to read from.
    52                                                       run_check: Option to check for the existence and proper shapes of
    53                                                           required parameters after reading in the file. Default is True.
    54                                                       run_check_acceptability: Option to check acceptable range of the values of
    55                                                           required parameters after reading in the file. Default is True.
    56                                                   """
    57         1        12587  12587.0      0.1          F = fits.open(filename)
    58         1           10     10.0      0.0          D = F[0]  # assumes the visibilities are in the primary hdu
    59         1         2860   2860.0      0.0          hdr = D.header.copy()
    60         1          191    191.0      0.0          hdunames = self._indexhdus(F)  # find the rest of the tables
    61                                           
    62                                                   # astropy.io fits reader scales date according to relevant PZER0 (?)
    63         1      3456621 3456621.0     24.0          time0_array = D.data['DATE']
    64         1            4      4.0      0.0          try:
    65                                                       # uvfits standard is to have 2 DATE parameters, both floats:
    66                                                       # DATE (full day) and _DATE (fractional day)
    67         1          165    165.0      0.0              time1_array = D.data['_DATE']
    68                                                       self.time_array = (time0_array.astype(np.double) +
    69                                                                          time1_array.astype(np.double))
    70         1            4      4.0      0.0          except(KeyError):
    71                                                       # cotter uvfits files have one DATE that is a double
    72         1           11     11.0      0.0              self.time_array = time0_array
    73         1           29     29.0      0.0              if np.finfo(time0_array[0]).precision < 5:
    74                                                           raise ValueError('JDs in this file are not precise to '
    75                                                                            'better than a second.')
    76         1            7      7.0      0.0              if (np.finfo(time0_array[0]).precision > 5 and
    77         1            7      7.0      0.0                      np.finfo(time0_array[0]).precision < 8):
    78                                                           warnings.warn('The JDs in this file have sub-second '
    79                                                                         'precision, but not sub-millisecond. '
    80                                                                         'Use with caution.')
    81                                           
    82         1         9420   9420.0      0.1          self.Ntimes = len(np.unique(self.time_array))
    83                                           
    84                                                   # if antenna arrays are present, use them. otherwise use baseline array
    85         1            3      3.0      0.0          try:
    86                                                       # Note: uvfits antennas are 1 indexed,
    87                                                       # need to subtract one to get to 0-indexed
    88         1       127985 127985.0      0.9              self.ant_1_array = np.int32(D.data.field('ANTENNA1')) - 1
    89         1        27542  27542.0      0.2              self.ant_2_array = np.int32(D.data.field('ANTENNA2')) - 1
    90         1        26001  26001.0      0.2              subarray = np.int32(D.data.field('SUBARRAY')) - 1
    91                                                       # error on files with multiple subarrays
    92         1        40825  40825.0      0.3              if len(set(subarray)) > 1:
    93                                                           raise ValueError('This file appears to have multiple subarray '
    94                                                                            'values; only files with one subarray are '
    95                                                                            'supported.')
    96                                                   except(KeyError):
    97                                                       # cannot set this to be the baseline array because it uses the
    98                                                       # 256 convention, not our 2048 convention
    99                                                       bl_input_array = np.int64(D.data.field('BASELINE'))
   100                                           
   101                                                       # get antenna arrays based on uvfits baseline array
   102                                                       self.ant_1_array, self.ant_2_array = \
   103                                                           self.baseline_to_antnums(bl_input_array)
   104                                                   # check for multi source files
   105         1            3      3.0      0.0          try:
   106         1          170    170.0      0.0              source = D.data.field('SOURCE')
   107                                                       if len(set(source)) > 1:
   108                                                           raise ValueError('This file has multiple sources. Only single '
   109                                                                            'source observations are supported.')
   110         1            3      3.0      0.0          except:
   111         1            2      2.0      0.0              pass
   112                                           
   113                                                   # get self.baseline_array using our convention
   114                                                   self.baseline_array = \
   115         1            9      9.0      0.0              self.antnums_to_baseline(self.ant_1_array,
   116         1         3254   3254.0      0.0                                       self.ant_2_array)
   117         1        23001  23001.0      0.2          self.Nbls = len(np.unique(self.baseline_array))
   118                                           
   119                                                   # initialize internal variables based on the antenna table
   120         1        85620  85620.0      0.6          self.Nants_data = int(len(np.unique(self.ant_1_array.tolist() + self.ant_2_array.tolist())))
   121                                           
   122         1           16     16.0      0.0          self.set_phased()
   123                                                   # check if we have an spw dimension
   124         1          340    340.0      0.0          if hdr.pop('NAXIS') == 7:
   125         1           64     64.0      0.0              if hdr['NAXIS5'] > 1:
   126                                                           raise ValueError('Sorry.  Files with more than one spectral' +
   127                                                                            'window (spw) are not yet supported. A ' +
   128                                                                            'great project for the interested student!')
   129         1          338    338.0      0.0              self.data_array = (D.data.field('DATA')[:, 0, 0, :, :, :, 0] +
   130         1      8472373 8472373.0     58.7                                 1j * D.data.field('DATA')[:, 0, 0, :, :, :, 1])
   131         1       546361 546361.0      3.8              self.flag_array = (D.data.field('DATA')[:, 0, 0, :, :, :, 2] <= 0)
   132         1       781385 781385.0      5.4              self.nsample_array = np.abs(D.data.field('DATA')[:, 0, 0, :, :, :, 2])
   133         1          345    345.0      0.0              self.Nspws = hdr.pop('NAXIS5')
   134         1            9      9.0      0.0              assert(self.Nspws == self.data_array.shape[1])
   135                                           
   136                                                       # the axis number for phase center depends on if the spw exists
   137                                                       # subtract 1 to be zero-indexed
   138         1          614    614.0      0.0              self.spw_array = np.int32(self._gethduaxis(D, 5)) - 1
   139                                           
   140         1          294    294.0      0.0              self.phase_center_ra_degrees = np.float(hdr.pop('CRVAL6'))
   141         1          282    282.0      0.0              self.phase_center_dec_degrees = np.float(hdr.pop('CRVAL7'))
   142                                                   else:
   143                                                       # in many uvfits files the spw axis is left out,
   144                                                       # here we put it back in so the dimensionality stays the same
   145                                                       self.data_array = (D.data.field('DATA')[:, 0, 0, :, :, 0] +
   146                                                                          1j * D.data.field('DATA')[:, 0, 0, :, :, 1])
   147                                                       self.data_array = self.data_array[:, np.newaxis, :, :]
   148                                                       self.flag_array = (D.data.field('DATA')[:, 0, 0, :, :, 2] <= 0)
   149                                                       self.flag_array = self.flag_array[:, np.newaxis, :, :]
   150                                                       self.nsample_array = np.abs(D.data.field('DATA')[:, 0, 0, :, :, 2])
   151                                                       self.nsample_array = (self.nsample_array[:, np.newaxis, :, :])
   152                                           
   153                                                       # the axis number for phase center depends on if the spw exists
   154                                                       self.Nspws = 1
   155                                                       self.spw_array = np.array([0])
   156                                           
   157                                                       self.phase_center_ra_degrees = np.float(hdr.pop('CRVAL5'))
   158                                                       self.phase_center_dec_degrees = np.float(hdr.pop('CRVAL6'))
   159                                           
   160                                                   # get shapes
   161         1          191    191.0      0.0          self.Nfreqs = hdr.pop('NAXIS4')
   162         1            6      6.0      0.0          if self.data_array.shape[2] != self.Nfreqs:
   163                                                       warnings.warn('Nfreqs does not match the number of frequencies in the data')
   164         1          188    188.0      0.0          self.Npols = hdr.pop('NAXIS3')
   165         1            6      6.0      0.0          if self.data_array.shape[3] != self.Npols:
   166                                                       warnings.warn('npols={npols} but found {l} pols in data file'.format(
   167                                                           npols=self.Npols, l=len(self.polarization_array)))
   168         1          184    184.0      0.0          self.Nblts = hdr.pop('GCOUNT')
   169         1            5      5.0      0.0          if self.data_array.shape[0] != self.Nblts:
   170                                                       warnings.warn('Nblts does not match the number of unique blts in the data')
   171                                           
   172                                                   # read baseline vectors in units of seconds, return in meters
   173         1          304    304.0      0.0          self.uvw_array = (np.array(np.stack((D.data.field('UU'),
   174         1          235    235.0      0.0                                               D.data.field('VV'),
   175         1        34285  34285.0      0.2                                               D.data.field('WW')))) *
   176         1        10224  10224.0      0.1                            const.c.to('m/s').value).T
   177                                           
   178         1          618    618.0      0.0          self.freq_array = self._gethduaxis(D, 4)
   179         1          292    292.0      0.0          self.channel_width = hdr.pop('CDELT4')
   180                                           
   181         1            3      3.0      0.0          try:
   182         1          308    308.0      0.0              self.integration_time = float(D.data.field('INTTIM')[0])
   183                                                   except:
   184                                                       if self.Ntimes > 1:
   185                                                           self.integration_time = \
   186                                                               float(np.diff(np.sort(list(set(self.time_array))))[0]) * 86400
   187                                                       else:
   188                                                           raise ValueError('integration time not specified and only '
   189                                                                            'one time present')
   190                                           
   191         1           11     11.0      0.0          self.freq_array.shape = (self.Nspws,) + self.freq_array.shape
   192                                           
   193         1          498    498.0      0.0          self.polarization_array = np.int32(self._gethduaxis(D, 3))
   194                                           
   195                                                   # other info -- not required but frequently used
   196         1          270    270.0      0.0          self.object_name = hdr.pop('OBJECT', None)
   197         1          256    256.0      0.0          self.telescope_name = hdr.pop('TELESCOP', None)
   198         1          245    245.0      0.0          self.instrument = hdr.pop('INSTRUME', None)
   199         1          261    261.0      0.0          latitude_degrees = hdr.pop('LAT', None)
   200         1          258    258.0      0.0          longitude_degrees = hdr.pop('LON', None)
   201         1          252    252.0      0.0          altitude = hdr.pop('ALT', None)
   202         1          319    319.0      0.0          self.history = str(hdr.get('HISTORY', ''))
   203         1            8      8.0      0.0          if self.pyuvdata_version_str not in self.history.replace('\n', ''):
   204                                                       self.history += self.pyuvdata_version_str
   205         6          525     87.5      0.0          while 'HISTORY' in hdr.keys():
   206         5          473     94.6      0.0              hdr.remove('HISTORY')
   207         1          245    245.0      0.0          self.vis_units = hdr.pop('BUNIT', 'UNCALIB')
   208         1          254    254.0      0.0          self.phase_center_epoch = hdr.pop('EPOCH', None)
   209                                           
   210                                                   # remove standard FITS header items that are still around
   211         1            3      3.0      0.0          std_fits_substrings = ['SIMPLE', 'BITPIX', 'EXTEND', 'BLOCKED',
   212         1            3      3.0      0.0                                 'GROUPS', 'PCOUNT', 'BSCALE', 'BZERO', 'NAXIS',
   213         1            3      3.0      0.0                                 'PTYPE', 'PSCAL', 'PZERO', 'CTYPE', 'CRVAL',
   214         1            3      3.0      0.0                                 'CRPIX', 'CDELT', 'CROTA']
   215        57          230      4.0      0.0          for key in hdr.keys():
   216      1008         2689      2.7      0.0              for sub in std_fits_substrings:
   217       952         2834      3.0      0.0                  if key.find(sub) > -1:
   218        55         3645     66.3      0.0                      hdr.remove(key)
   219                                           
   220                                                   # find all the remaining header items and keep them as extra_keywords
   221         2           10      5.0      0.0          for key in hdr:
   222         1            3      3.0      0.0              if key == '':
   223                                                           continue
   224         1            3      3.0      0.0              if key == 'COMMENT':
   225                                                           self.extra_keywords[key] = str(hdr.get(key))
   226                                                       else:
   227         1          147    147.0      0.0                  self.extra_keywords[key] = hdr.get(key)
   228                                           
   229                                                   # READ the antenna table
   230         1           12     12.0      0.0          ant_hdu = F[hdunames['AIPS AN']]
   231                                           
   232                                                   # stuff in columns
   233         1        15194  15194.0      0.1          ant_names = ant_hdu.data.field('ANNAME').tolist()
   234         1            7      7.0      0.0          self.antenna_names = []
   235       129          350      2.7      0.0          for name in ant_names:
   236       128          509      4.0      0.0              self.antenna_names.append(name.replace('\x00!', ''))
   237                                           
   238                                                   # subtract one to get to 0-indexed values rather than 1-indexed values
   239         1          278    278.0      0.0          self.antenna_numbers = ant_hdu.data.field('NOSTA') - 1
   240                                           
   241         1            7      7.0      0.0          self.Nants_telescope = len(self.antenna_numbers)
   242                                           
   243                                                   # stuff in the header
   244         1            4      4.0      0.0          if self.telescope_name is None:
   245                                                       self.telescope_name = ant_hdu.header['ARRNAM']
   246                                           
   247         1            3      3.0      0.0          try:
   248         1           69     69.0      0.0              xyz_telescope_frame = ant_hdu.header['FRAME']
   249                                                   except:
   250                                                       warnings.warn('Required Antenna frame keyword not set, '
   251                                                                     'setting to ????')
   252                                                       xyz_telescope_frame = '????'
   253                                           
   254                                                   # VLA incorrectly sets ARRAYX/ARRAYY/ARRAYZ to 0, and puts array center
   255                                                   # in the antenna positions themselves
   256                                           
   257         1          168    168.0      0.0          if (np.isclose(ant_hdu.header['ARRAYX'], 0) and
   258                                                           np.isclose(ant_hdu.header['ARRAYY'], 0) and
   259                                                           np.isclose(ant_hdu.header['ARRAYZ'], 0)):
   260                                                       x_telescope = np.mean(ant_hdu.data['STABXYZ'][:, 0])
   261                                                       y_telescope = np.mean(ant_hdu.data['STABXYZ'][:, 1])
   262                                                       z_telescope = np.mean(ant_hdu.data['STABXYZ'][:, 2])
   263                                                       self.antenna_positions = (ant_hdu.data.field('STABXYZ') -
   264                                                                                 np.array([x_telescope,
   265                                                                                           y_telescope,
   266                                                                                           z_telescope]))
   267                                           
   268                                                   else:
   269         1           61     61.0      0.0              x_telescope = ant_hdu.header['ARRAYX']
   270         1           55     55.0      0.0              y_telescope = ant_hdu.header['ARRAYY']
   271         1           53     53.0      0.0              z_telescope = ant_hdu.header['ARRAYZ']
   272         1          248    248.0      0.0              self.antenna_positions = ant_hdu.data.field('STABXYZ')
   273                                           
   274         1            3      3.0      0.0          if xyz_telescope_frame == 'ITRF':
   275         1           11     11.0      0.0              self.telescope_location = np.array([x_telescope, y_telescope, z_telescope])
   276                                                   else:
   277                                                       if latitude_degrees is not None and longitude_degrees is not None and altitude is not None:
   278                                                           self.telescope_location_lat_lon_alt_degrees = (latitude_degrees, longitude_degrees, altitude)
   279                                           
   280         1           63     63.0      0.0          self.gst0 = ant_hdu.header['GSTIA0']
   281         1           58     58.0      0.0          self.rdate = ant_hdu.header['RDATE']
   282         1           56     56.0      0.0          self.earth_omega = ant_hdu.header['DEGPDY']
   283         1           56     56.0      0.0          self.dut1 = ant_hdu.header['UT1UTC']
   284         1            3      3.0      0.0          try:
   285         1           24     24.0      0.0              self.timesys = ant_hdu.header['TIMESYS']
   286         1            3      3.0      0.0          except(KeyError):
   287                                                       # CASA misspells this one
   288         1           58     58.0      0.0              self.timesys = ant_hdu.header['TIMSYS']
   289                                           
   290         1            3      3.0      0.0          del(D)
   291                                           
   292         1            3      3.0      0.0          try:
   293         1          942    942.0      0.0              self.set_telescope_params()
   294                                                   except ValueError, ve:
   295                                                       warnings.warn(str(ve))
   296                                           
   297         1       614854 614854.0      4.3          self.set_lsts_from_time_array()
   298                                           
   299                                                   # check if object has all required UVParameters set
   300         1            3      3.0      0.0          if run_check:
   301         1       113947 113947.0      0.8              self.check(run_check_acceptability=run_check_acceptability)

Total time: 32.6351 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/uvfits.py
Function: write_uvfits at line 303

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   303                                               @profile
   304                                               def write_uvfits(self, filename, spoof_nonessential=False,
   305                                                                force_phase=False, run_check=True, run_check_acceptability=True):
   306                                                   """
   307                                                   Write the data to a uvfits file.
   308                                           
   309                                                   Args:
   310                                                       filename: The uvfits file to write to.
   311                                                       spoof_nonessential: Option to spoof the values of optional
   312                                                           UVParameters that are not set but are required for uvfits files.
   313                                                           Default is False.
   314                                                       force_phase: Option to automatically phase drift scan data to
   315                                                           zenith of the first timestamp. Default is False.
   316                                                       run_check: Option to check for the existence and proper shapes of
   317                                                           required parameters before writing the file. Default is True.
   318                                                       run_check_acceptability: Option to check acceptable range of the values of
   319                                                           required parameters before writing the file. Default is True.
   320                                                   """
   321         1            5      5.0      0.0          if run_check:
   322         1       124998 124998.0      0.4              self.check(run_check_acceptability=run_check_acceptability)
   323                                           
   324         1           14     14.0      0.0          if self.phase_type == 'phased':
   325                                                       pass
   326         1            5      5.0      0.0          elif self.phase_type == 'drift':
   327         1            3      3.0      0.0              if force_phase:
   328         1           53     53.0      0.0                  print('The data are in drift mode and do not have a '
   329                                                                 'defined phase center. Phasing to zenith of the first '
   330                                                                 'timestamp.')
   331         1     15964681 15964681.0     48.9                  self.phase_to_time(self.time_array[0])
   332                                                       else:
   333                                                           raise ValueError('The data are in drift mode. ' +
   334                                                                            'Set force_phase to true to phase the data ' +
   335                                                                            'to zenith of the first timestamp before ' +
   336                                                                            'writing a uvfits file.')
   337                                                   else:
   338                                                       raise ValueError('The phasing type of the data is unknown. '
   339                                                                        'Set the phase_type to drift or phased to '
   340                                                                        'reflect the phasing status of the data')
   341                                           
   342         1           62     62.0      0.0          freq_spacing = self.freq_array[0, 1:] - self.freq_array[0, :-1]
   343         1           64     64.0      0.0          if not np.isclose(np.min(freq_spacing), np.max(freq_spacing),
   344         1          150    150.0      0.0                            rtol=self._freq_array.tols[0], atol=self._freq_array.tols[1]):
   345                                                       raise ValueError('The frequencies are not evenly spaced (probably '
   346                                                                        'because of a select operation). The uvfits format '
   347                                                                        'does not support unevenly spaced frequencies.')
   348         1           18     18.0      0.0          if not np.isclose(np.max(freq_spacing), self.channel_width,
   349         1           83     83.0      0.0                            rtol=self._freq_array.tols[0], atol=self._freq_array.tols[1]):
   350                                                       raise ValueError('The frequencies are separated by more than their '
   351                                                                        'channel width (probably because of a select operation). '
   352                                                                        'The uvfits format does not support frequencies '
   353                                                                        'that are spaced by more than their channel width.')
   354         1           12     12.0      0.0          if self.Npols > 2:
   355                                                       pol_spacing = self.polarization_array[1:] - self.polarization_array[:-1]
   356                                                       if np.min(pol_spacing) < np.max(pol_spacing):
   357                                                           raise ValueError('The polarization values are not evenly spaced (probably '
   358                                                                            'because of a select operation). The uvfits format '
   359                                                                            'does not support unevenly spaced polarizations.')
   360                                           
   361        11         8640    785.5      0.0          for p in self.extra():
   362        10           36      3.6      0.0              param = getattr(self, p)
   363        10           40      4.0      0.0              if param.name in self.uvfits_required_extra:
   364         6           23      3.8      0.0                  if param.value is None:
   365         6           19      3.2      0.0                      if spoof_nonessential:
   366                                                                   # spoof extra keywords required for uvfits
   367         6           25      4.2      0.0                          if isinstance(param, uvp.AntPositionParameter):
   368         1           19     19.0      0.0                              param.apply_spoof(self, 'Nants_telescope')
   369                                                                   else:
   370         5           28      5.6      0.0                              param.apply_spoof()
   371         6           25      4.2      0.0                          setattr(self, p, param)
   372                                                               else:
   373                                                                   raise ValueError('Required attribute {attribute} '
   374                                                                                    'for uvfits not defined. Define or '
   375                                                                                    'set spoof_nonessential to True to '
   376                                                                                    'spoof this attribute.'
   377                                                                                    .format(attribute=p))
   378                                           
   379         1            5      5.0      0.0          weights_array = self.nsample_array * \
   380         1      2142990 2142990.0      6.6              np.where(self.flag_array, -1, 1)
   381         1           39     39.0      0.0          data_array = self.data_array[:, np.newaxis, np.newaxis, :, :, :, np.newaxis]
   382         1            4      4.0      0.0          weights_array = weights_array[:, np.newaxis, np.newaxis, :, :, :,
   383         1            5      5.0      0.0                                        np.newaxis]
   384                                                   # uvfits_array_data shape will be  (Nblts,1,1,[Nspws],Nfreqs,Npols,3)
   385         1            6      6.0      0.0          uvfits_array_data = np.concatenate([data_array.real,
   386         1            4      4.0      0.0                                             data_array.imag,
   387         1      4178199 4178199.0     12.8                                             weights_array], axis=6)
   388                                           
   389         1         5501   5501.0      0.0          uvw_array_sec = self.uvw_array / const.c.to('m/s').value
   390                                                   # jd_midnight = np.floor(self.time_array[0] - 0.5) + 0.5
   391         1           39     39.0      0.0          tzero = np.float32(self.time_array[0])
   392                                           
   393                                                   # uvfits convention is that time_array + relevant PZERO = actual JD
   394                                                   # We are setting PZERO4 = float32(first time of observation)
   395         1         1438   1438.0      0.0          time_array = np.float32(self.time_array - np.float64(tzero))
   396                                           
   397         1          198    198.0      0.0          int_time_array = (np.zeros_like((time_array), dtype=np.float) +
   398         1          456    456.0      0.0                            self.integration_time)
   399                                           
   400         1           10     10.0      0.0          baselines_use = self.antnums_to_baseline(self.ant_1_array,
   401         1            5      5.0      0.0                                                   self.ant_2_array,
   402         1         3807   3807.0      0.0                                                   attempt256=True)
   403                                                   # Set up dictionaries for populating hdu
   404                                                   # Note that uvfits antenna arrays are 1-indexed so we add 1
   405                                                   # to our 0-indexed arrays
   406         1           12     12.0      0.0          group_parameter_dict = {'UU      ': uvw_array_sec[:, 0],
   407         1            4      4.0      0.0                                  'VV      ': uvw_array_sec[:, 1],
   408         1            4      4.0      0.0                                  'WW      ': uvw_array_sec[:, 2],
   409         1            4      4.0      0.0                                  'DATE    ': time_array,
   410         1            4      4.0      0.0                                  'BASELINE': baselines_use,
   411         1          461    461.0      0.0                                  'ANTENNA1': self.ant_1_array + 1,
   412         1          515    515.0      0.0                                  'ANTENNA2': self.ant_2_array + 1,
   413         1          409    409.0      0.0                                  'SUBARRAY': np.ones_like(self.ant_1_array),
   414         1            4      4.0      0.0                                  'INTTIM': int_time_array}
   415         1            5      5.0      0.0          pscal_dict = {'UU      ': 1.0, 'VV      ': 1.0, 'WW      ': 1.0,
   416         1            4      4.0      0.0                        'DATE    ': 1.0, 'BASELINE': 1.0, 'ANTENNA1': 1.0,
   417         1            3      3.0      0.0                        'ANTENNA2': 1.0, 'SUBARRAY': 1.0, 'INTTIM': 1.0}
   418         1            5      5.0      0.0          pzero_dict = {'UU      ': 0.0, 'VV      ': 0.0, 'WW      ': 0.0,
   419         1            4      4.0      0.0                        'DATE    ': tzero, 'BASELINE': 0.0, 'ANTENNA1': 0.0,
   420         1            4      4.0      0.0                        'ANTENNA2': 0.0, 'SUBARRAY': 0.0, 'INTTIM': 0.0}
   421                                           
   422                                                   # list contains arrays of [u,v,w,date,baseline];
   423                                                   # each array has shape (Nblts)
   424         1          515    515.0      0.0          if (np.max(self.ant_1_array) < 255 and
   425         1          458    458.0      0.0                  np.max(self.ant_2_array) < 255):
   426                                                       # if the number of antennas is less than 256 then include both the
   427                                                       # baseline array and the antenna arrays in the group parameters.
   428                                                       # Otherwise just use the antenna arrays
   429         1            4      4.0      0.0              parnames_use = ['UU      ', 'VV      ', 'WW      ',
   430         1            4      4.0      0.0                              'DATE    ', 'BASELINE', 'ANTENNA1',
   431         1            5      5.0      0.0                              'ANTENNA2', 'SUBARRAY', 'INTTIM']
   432                                                   else:
   433                                                       parnames_use = ['UU      ', 'VV      ', 'WW      ', 'DATE    ',
   434                                                                       'ANTENNA1', 'ANTENNA2', 'SUBARRAY', 'INTTIM']
   435                                           
   436         1            4      4.0      0.0          group_parameter_list = [group_parameter_dict[parname] for
   437        10           42      4.2      0.0                                  parname in parnames_use]
   438         1            5      5.0      0.0          hdu = fits.GroupData(uvfits_array_data, parnames=parnames_use,
   439         1      1963159 1963159.0      6.0                               pardata=group_parameter_list, bitpix=-32)
   440         1         9105   9105.0      0.0          hdu = fits.GroupsHDU(hdu)
   441                                           
   442        10           42      4.2      0.0          for i, key in enumerate(parnames_use):
   443         9         1908    212.0      0.0              hdu.header['PSCAL' + str(i + 1) + '  '] = pscal_dict[key]
   444         9         1969    218.8      0.0              hdu.header['PZERO' + str(i + 1) + '  '] = pzero_dict[key]
   445                                           
   446                                                   # ISO string of first time in self.time_array
   447         1           13     13.0      0.0          hdu.header['DATE-OBS'] = Time(self.time_array[0], scale='utc',
   448         1          759    759.0      0.0                                        format='jd').isot
   449                                           
   450         1          225    225.0      0.0          hdu.header['CTYPE2  '] = 'COMPLEX '
   451         1          226    226.0      0.0          hdu.header['CRVAL2  '] = 1.0
   452         1          224    224.0      0.0          hdu.header['CRPIX2  '] = 1.0
   453         1          222    222.0      0.0          hdu.header['CDELT2  '] = 1.0
   454                                           
   455         1          221    221.0      0.0          hdu.header['CTYPE3  '] = 'STOKES  '
   456         1          236    236.0      0.0          hdu.header['CRVAL3  '] = self.polarization_array[0]
   457         1          228    228.0      0.0          hdu.header['CRPIX3  '] = 1.0
   458         1            4      4.0      0.0          try:
   459         1          252    252.0      0.0              hdu.header['CDELT3  '] = np.diff(self.polarization_array)[0]
   460                                                   except(IndexError):
   461                                                       hdu.header['CDELT3  '] = 1.0
   462                                           
   463         1          226    226.0      0.0          hdu.header['CTYPE4  '] = 'FREQ    '
   464         1          247    247.0      0.0          hdu.header['CRVAL4  '] = self.freq_array[0, 0]
   465         1          232    232.0      0.0          hdu.header['CRPIX4  '] = 1.0
   466         1          262    262.0      0.0          hdu.header['CDELT4  '] = np.diff(self.freq_array[0])[0]
   467                                           
   468         1          230    230.0      0.0          hdu.header['CTYPE5  '] = 'IF      '
   469         1          238    238.0      0.0          hdu.header['CRVAL5  '] = 1.0
   470         1          234    234.0      0.0          hdu.header['CRPIX5  '] = 1.0
   471         1          235    235.0      0.0          hdu.header['CDELT5  '] = 1.0
   472                                           
   473         1          232    232.0      0.0          hdu.header['CTYPE6  '] = 'RA'
   474         1          250    250.0      0.0          hdu.header['CRVAL6  '] = self.phase_center_ra_degrees
   475                                           
   476         1          235    235.0      0.0          hdu.header['CTYPE7  '] = 'DEC'
   477         1          243    243.0      0.0          hdu.header['CRVAL7  '] = self.phase_center_dec_degrees
   478                                           
   479         1          237    237.0      0.0          hdu.header['BUNIT   '] = self.vis_units
   480         1          247    247.0      0.0          hdu.header['BSCALE  '] = 1.0
   481         1          243    243.0      0.0          hdu.header['BZERO   '] = 0.0
   482                                           
   483         1          246    246.0      0.0          hdu.header['OBJECT  '] = self.object_name
   484         1          241    241.0      0.0          hdu.header['TELESCOP'] = self.telescope_name
   485         1          362    362.0      0.0          hdu.header['LAT     '] = self.telescope_location_lat_lon_alt_degrees[0]
   486         1          329    329.0      0.0          hdu.header['LON     '] = self.telescope_location_lat_lon_alt_degrees[1]
   487         1          329    329.0      0.0          hdu.header['ALT     '] = self.telescope_location_lat_lon_alt[2]
   488         1          258    258.0      0.0          hdu.header['INSTRUME'] = self.instrument
   489         1          255    255.0      0.0          hdu.header['EPOCH   '] = float(self.phase_center_epoch)
   490                                           
   491         3           16      5.3      0.0          for line in self.history.splitlines():
   492         2         1092    546.0      0.0              hdu.header.add_history(line)
   493                                           
   494                                                   # end standard keywords; begin user-defined keywords
   495         1            8      8.0      0.0          for key, value in self.extra_keywords.iteritems():
   496                                                       # header keywords have to be 8 characters or less
   497                                                       keyword = key[:8].upper()
   498                                                       # print "keyword=-value-", keyword+'='+'-'+str(value)+'-'
   499                                                       if keyword == 'COMMENT':
   500                                                           for line in value.splitlines():
   501                                                               hdu.header.add_comment(line)
   502                                                       else:
   503                                                           hdu.header[keyword] = value
   504                                           
   505                                                   # ADD the ANTENNA table
   506         1           10     10.0      0.0          staxof = np.zeros(self.Nants_telescope)
   507                                           
   508                                                   # 0 specifies alt-az, 6 would specify a phased array
   509         1            7      7.0      0.0          mntsta = np.zeros(self.Nants_telescope)
   510                                           
   511                                                   # beware, X can mean just about anything
   512         1           53     53.0      0.0          poltya = np.full((self.Nants_telescope), 'X', dtype=np.object_)
   513         1           25     25.0      0.0          polaa = [90.0] + np.zeros(self.Nants_telescope)
   514         1           18     18.0      0.0          poltyb = np.full((self.Nants_telescope), 'Y', dtype=np.object_)
   515         1           17     17.0      0.0          polab = [0.0] + np.zeros(self.Nants_telescope)
   516                                           
   517         1            4      4.0      0.0          col1 = fits.Column(name='ANNAME', format='8A',
   518         1          492    492.0      0.0                             array=self.antenna_names)
   519         1            4      4.0      0.0          col2 = fits.Column(name='STABXYZ', format='3D',
   520         1          484    484.0      0.0                             array=self.antenna_positions)
   521                                                   # convert to 1-indexed from 0-indexed indicies
   522         1            5      5.0      0.0          col3 = fits.Column(name='NOSTA', format='1J',
   523         1          422    422.0      0.0                             array=self.antenna_numbers + 1)
   524         1          398    398.0      0.0          col4 = fits.Column(name='MNTSTA', format='1J', array=mntsta)
   525         1          393    393.0      0.0          col5 = fits.Column(name='STAXOF', format='1E', array=staxof)
   526         1          425    425.0      0.0          col6 = fits.Column(name='POLTYA', format='1A', array=poltya)
   527         1          419    419.0      0.0          col7 = fits.Column(name='POLAA', format='1E', array=polaa)
   528                                                   # col8 = fits.Column(name='POLCALA', format='3E', array=polcala)
   529         1          393    393.0      0.0          col9 = fits.Column(name='POLTYB', format='1A', array=poltyb)
   530         1          392    392.0      0.0          col10 = fits.Column(name='POLAB', format='1E', array=polab)
   531                                                   # col11 = fits.Column(name='POLCALB', format='3E', array=polcalb)
   532                                                   # note ORBPARM is technically required, but we didn't put it in
   533                                           
   534         1            4      4.0      0.0          cols = fits.ColDefs([col1, col2, col3, col4, col5, col6, col7, col9,
   535         1         1626   1626.0      0.0                               col10])
   536                                           
   537         1        24300  24300.0      0.1          ant_hdu = fits.BinTableHDU.from_columns(cols)
   538                                           
   539         1          208    208.0      0.0          ant_hdu.header['EXTNAME'] = 'AIPS AN'
   540         1          198    198.0      0.0          ant_hdu.header['EXTVER'] = 1
   541                                           
   542                                                   # write XYZ coordinates if not already defined
   543         1          227    227.0      0.0          ant_hdu.header['ARRAYX'] = self.telescope_location[0]
   544         1          218    218.0      0.0          ant_hdu.header['ARRAYY'] = self.telescope_location[1]
   545         1          216    216.0      0.0          ant_hdu.header['ARRAYZ'] = self.telescope_location[2]
   546         1          205    205.0      0.0          ant_hdu.header['FRAME'] = 'ITRF'
   547         1          216    216.0      0.0          ant_hdu.header['GSTIA0'] = self.gst0
   548         1          220    220.0      0.0          ant_hdu.header['FREQ'] = self.freq_array[0, 0]
   549         1          205    205.0      0.0          ant_hdu.header['RDATE'] = self.rdate
   550         1          214    214.0      0.0          ant_hdu.header['UT1UTC'] = self.dut1
   551                                           
   552         1          211    211.0      0.0          ant_hdu.header['TIMSYS'] = self.timesys
   553         1            5      5.0      0.0          if self.timesys == 'IAT':
   554                                                       warnings.warn('This file has an "IAT" time system. Files of '
   555                                                                     'this type are not properly supported')
   556         1          215    215.0      0.0          ant_hdu.header['ARRNAM'] = self.telescope_name
   557         1          207    207.0      0.0          ant_hdu.header['NO_IF'] = self.Nspws
   558         1          218    218.0      0.0          ant_hdu.header['DEGPDY'] = self.earth_omega
   559                                                   # ant_hdu.header['IATUTC'] = 35.
   560                                           
   561                                                   # set mandatory parameters which are not supported by this object
   562                                                   # (or that we just don't understand)
   563         1          207    207.0      0.0          ant_hdu.header['NUMORB'] = 0
   564                                           
   565                                                   # note: Bart had this set to 3. We've set it 0 after aips 117. -jph
   566         1          211    211.0      0.0          ant_hdu.header['NOPCAL'] = 0
   567                                           
   568         1          214    214.0      0.0          ant_hdu.header['POLTYPE'] = 'X-Y LIN'
   569                                           
   570                                                   # note: we do not support the concept of "frequency setups"
   571                                                   # -- lists of spws given in a SU table.
   572         1          209    209.0      0.0          ant_hdu.header['FREQID'] = -1
   573                                           
   574                                                   # if there are offsets in images, this could be the culprit
   575         1          284    284.0      0.0          ant_hdu.header['POLARX'] = 0.0
   576         1          228    228.0      0.0          ant_hdu.header['POLARY'] = 0.0
   577                                           
   578         1          216    216.0      0.0          ant_hdu.header['DATUTC'] = 0  # ONLY UTC SUPPORTED
   579                                           
   580                                                   # we always output right handed coordinates
   581         1          219    219.0      0.0          ant_hdu.header['XYZHAND'] = 'RIGHT'
   582                                           
   583                                                   # ADD the FQ table
   584                                                   # skipping for now and limiting to a single spw
   585                                           
   586                                                   # write the file
   587         1           92     92.0      0.0          hdulist = fits.HDUList(hdus=[hdu, ant_hdu])
   588         1            9      9.0      0.0          if float(astropy.__version__[0:3]) < 1.3:
   589         1      8180674 8180674.0     25.1              hdulist.writeto(filename, clobber=True)
   590                                                   else:
   591                                                       hdulist.writeto(filename, overwrite=True)

back to top