Skip to main content
  • Home
  • Development
  • Documentation
  • Donate
  • Operational login
  • Browse the archive

swh logo
SoftwareHeritage
Software
Heritage
Archive
Features
  • Search

  • Downloads

  • Save code now

  • Add forge now

  • Help

https://github.com/RadioAstronomySoftwareGroup/pyuvdata
05 July 2025, 16:03:44 UTC
  • Code
  • Branches (68)
  • Releases (1)
  • Visits
    • Branches
    • Releases
    • HEAD
    • refs/heads/add_bitshuffle
    • refs/heads/add_feko_read
    • refs/heads/allow-extrap-in-beam
    • refs/heads/allow_exclude_ants
    • refs/heads/allow_power_units
    • refs/heads/bitshuffle_v2
    • refs/heads/combine_uvdata_func
    • refs/heads/double_prec_uvws
    • refs/heads/faster-interp
    • refs/heads/faster-uvh5-indexing
    • refs/heads/fits_speedup
    • refs/heads/fix-cst-beam-read
    • refs/heads/frame_attr
    • refs/heads/indexIng_tools
    • refs/heads/main
    • refs/heads/miriad_tweaks_v3
    • refs/heads/mwa_van_vleck_2
    • refs/heads/ovro-lwa
    • refs/heads/py03
    • refs/heads/simple-set-antdiam
    • refs/heads/sma_dev
    • refs/heads/snap_convert
    • refs/heads/use_gh_cache
    • refs/tags/V2.1.5
    • refs/tags/v1.1
    • refs/tags/v1.2
    • refs/tags/v1.3
    • refs/tags/v1.4
    • refs/tags/v1.5
    • refs/tags/v2.0.0
    • refs/tags/v2.0.1
    • refs/tags/v2.0.2
    • refs/tags/v2.1.0
    • refs/tags/v2.1.1
    • refs/tags/v2.1.2
    • refs/tags/v2.1.3
    • refs/tags/v2.1.4
    • refs/tags/v2.2.0
    • refs/tags/v2.2.1
    • refs/tags/v2.2.10
    • refs/tags/v2.2.11
    • refs/tags/v2.2.12
    • refs/tags/v2.2.2
    • refs/tags/v2.2.3
    • refs/tags/v2.2.4
    • refs/tags/v2.2.5
    • refs/tags/v2.2.6
    • refs/tags/v2.2.7
    • refs/tags/v2.2.8
    • refs/tags/v2.2.9
    • refs/tags/v2.3.0
    • refs/tags/v2.3.1
    • refs/tags/v2.3.2
    • refs/tags/v2.3.3
    • refs/tags/v2.4.0
    • refs/tags/v2.4.1
    • refs/tags/v2.4.2
    • refs/tags/v2.4.3
    • refs/tags/v2.4.4
    • refs/tags/v2.4.5
    • refs/tags/v3.0.0
    • refs/tags/v3.1.0
    • refs/tags/v3.1.1
    • refs/tags/v3.1.2
    • refs/tags/v3.1.3
    • refs/tags/v3.2.0
    • refs/tags/v3.2.1
    • refs/tags/v3.2.2
    • v1.0
  • aacb583
  • /
  • pyuvdata
  • /
  • tests
  • /
  • Line_Profiler_full.out
Raw File Download
Take a new snapshot of a software origin

If the archived software origin currently browsed is not synchronized with its upstream version (for instance when new commits have been issued), you can explicitly request Software Heritage to take a new snapshot of it.

Use the form below to proceed. Once a request has been submitted and accepted, it will be processed as soon as possible. You can then check its processing state by visiting this dedicated page.
swh spinner

Processing "take a new snapshot" request ...

Permalinks

To reference or cite the objects present in the Software Heritage archive, permalinks based on SoftWare Hash IDentifiers (SWHIDs) must be used.
Select below a type of object currently browsed in order to display its associated SWHID and permalink.

  • content
  • directory
  • revision
  • snapshot
origin badgecontent badge Iframe embedding
swh:1:cnt:90f7bf37023e6ef0584f5af0c376e347151fb859
origin badgedirectory badge Iframe embedding
swh:1:dir:c2b4186780122327d7d2d231e2556cd6ebdd52b6
origin badgerevision badge
swh:1:rev:c5a1b463f29b8164cda381e97e3d93a549be4179
origin badgesnapshot badge
swh:1:snp:bcfc5cc16e5c167943740f8c7459e8b553b12133
Citations

This interface enables to generate software citations, provided that the root directory of browsed objects contains a citation.cff or codemeta.json file.
Select below a type of object currently browsed in order to generate citations for them.

  • content
  • directory
  • revision
  • snapshot
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
Generate software citation in BibTex format (requires biblatex-software package)
Generating citation ...
Tip revision: c5a1b463f29b8164cda381e97e3d93a549be4179 authored by Bryna Hazelton on 29 April 2020, 15:53:31 UTC
Update changelog for new version, minor changelog correction.
Tip revision: c5a1b46
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

Software Heritage — Copyright (C) 2015–2025, The Software Heritage developers. License: GNU AGPLv3+.
The source code of Software Heritage itself is available on our development forge.
The source code files archived by Software Heritage are available under their own copyright and licenses.
Terms of use: Archive access, API— Contact— JavaScript license information— Web API