Refactor of grid code

[low-sky]

Hi @rfriesen @jpinedaf

First, don't merge quite yet. I wanted to get this out for discussion.

My fix to the spikes in the (3,3) window is to completely rewrite the gridding code. Or rather, what I've done is to use some code that has been evolving and refactored into a stand-alone package called gbtpipe that includes a lot of the common framework used across spectral line surveys. The code hasn't changed too much, but it does have some additional flagging that I had already built for KEYSTONE, namely flagSpike which flags spikes in VEGAS. I'm running a full-end-to-end right now for L1688 but it seems to be working as intended.

The gbtpipe repo is at https://github.com/low-sky/gbtpipe/ with relevant changes to the gridding code here.

Check out the files in /lustre/pipeline/scratch/GAS/images/L1688 as they come in over the next day.

[rfriesen]

Hey Erik, Thanks for working on this! I've taken a look at the L1688 data. It looks like the new gridding code is getting rid of the spikes, at least as far as I've checked in the (3,3) cube. It also seems to be getting rid of good signal, though. Here's a comparison of the same velocity channel in the DR1 vs. new cube for NH3 (1,1): (Same colour range, DR1 is on the right)

[low-sky]

Thanks; good check. The defaults must be too aggressive.

[low-sky]

The issue is that some feeds were producing a lot of scans flagged by rms flagging. Namely, if the rms in a spectrum is > 1.25 times the radiometer expectation, it was getting flagged. This is also new in this refactor. Some files lost 90% of their data, which appear to be okay otherwise. Maybe the Tsys data were bad? Regardless, I've upped the default GAS rejection threshold to 1.5 and this cleans up the issues with missing data.

I'm running some new versions now.

[low-sky]

Here's the comparable image using the new code.

I need to wait for the (3,3) to finish to make sure it's working as intended.

[rfriesen]

@low-sky Running through the gridding on the GB machines, and I ended up with an 'UnboundLocalError' in the gbtpipe calls re: the eulerFlag variable being referenced before assignment. It looks like I've got gbtpipe up-to-date on the GB machines, so am I missing something?

[low-sky]

Hi @rfriesen -- Okay; I pushed an untested fix to main (the default for the flag needs to be set to False). I'll try to test it ASAP, but if it clears things up; let me know.

[low-sky]

I squashed a couple more bugs in GBTpipe (pesky commas) and it looks to be running again. thanks for the heads-up.

[rfriesen]

Seems to be running now!

[rfriesen]

Coming up with an error again while trying to grid, this time with the header of the new cube:

IOError Traceback (most recent call last)

in () ----> 1 GAS.run_grid_regions.grid_SerAqu() /lustre/pipeline/scratch/GAS/python/GAS/GAS/run_grid_regions.pyc in grid_SerAqu(release) 92 startChannel = startChannel, endChannel = endChannel, 93 templateHeader=hd_temp, ---> 94 Sessions=mySessions, file_extension=file_extension) 95 96 /lustre/pipeline/scratch/GAS/python/GAS/GAS/gridregion.pyc in griddata(templateHeader, gridFunction, rootdir, region, dirname, startChannel, endChannel, doBaseline, baselineRegion, blorder, Sessions, file_extension, flagRMS, flagSpike, **kwargs) 191 outdir=outdir, outname=outname, 192 VlsrByCoord=VlsrByCoord, --> 193 flagSpike=flagSpike, **kwargs) 194 # Convolve the beam size up by 10% in size 195 gbtpipe.Gridding.postConvolve(outdir+outname) /lustre/pipeline/scratch/GAS/python/anaconda/lib/python2.7/site-packages/gbtpipe-0.1.3.dev34-py2.7.egg/gbtpipe/Gridding.pyc in griddata(filelist, pixPerBeam, templateHeader, gridFunction, startChannel, endChannel, doBaseline, baselineRegion, blorder, rebase, rebaseorder, beamSize, OnlineDoppler, flagRMS, flagRipple, flagSpike, rmsThresh, spikeThresh, projection, outdir, outname, dtype, **kwargs) 398 if templateHeader is None: 399 wcsdict = autoHeader(filelist, beamSize=beamSize, --> 400 pixPerBeam=pixPerBeam, projection=projection) 401 w.wcs.crpix = [wcsdict['CRPIX1'], wcsdict['CRPIX2'], crpix3] 402 w.wcs.cdelt = np.array([wcsdict['CDELT1'], wcsdict['CDELT2'], cdelt3]) /lustre/pipeline/scratch/GAS/python/anaconda/lib/python2.7/site-packages/gbtpipe-0.1.3.dev34-py2.7.egg/gbtpipe/Gridding.pyc in autoHeader(filelist, beamSize, pixPerBeam, projection) 175 DEClist = [] 176 for thisfile in filelist: --> 177 s = fits.getdata(thisfile) 178 try: 179 RAlist = RAlist + [s['CRVAL2']] /lustre/pipeline/scratch/GAS/python/anaconda/lib/python2.7/site-packages/astropy/io/fits/convenience.pyc in getdata(filename, *args, **kwargs) 205 if data is None and extidx == 0: 206 try: --> 207 hdu = hdulist[1] 208 data = hdu.data 209 except IndexError: /lustre/pipeline/scratch/GAS/python/anaconda/lib/python2.7/site-packages/astropy/io/fits/hdu/hdulist.pyc in __getitem__(self, key) 312 # instead 313 return self._try_while_unread_hdus(super(HDUList, self).__getitem__, --> 314 self._positive_index_of(key)) 315 316 def __contains__(self, item): /lustre/pipeline/scratch/GAS/python/anaconda/lib/python2.7/site-packages/astropy/io/fits/hdu/hdulist.pyc in _try_while_unread_hdus(self, func, index, *args, **kwargs) 1095 return func(index, *args, **kwargs) 1096 except Exception: -> 1097 if self._read_next_hdu(): 1098 continue 1099 else: /lustre/pipeline/scratch/GAS/python/anaconda/lib/python2.7/site-packages/astropy/io/fits/hdu/hdulist.pyc in _read_next_hdu(self) 1133 fileobj.seek(offset, os.SEEK_SET) 1134 -> 1135 hdu = _BaseHDU.readfrom(fileobj, **kwargs) 1136 except EOFError: 1137 self._read_all = True /lustre/pipeline/scratch/GAS/python/anaconda/lib/python2.7/site-packages/astropy/io/fits/hdu/base.pyc in readfrom(cls, fileobj, checksum, ignore_missing_end, **kwargs) 327 hdu = cls._readfrom_internal(fileobj, checksum=checksum, 328 ignore_missing_end=ignore_missing_end, --> 329 **kwargs) 330 331 # If the checksum had to be checked the data may have already been read /lustre/pipeline/scratch/GAS/python/anaconda/lib/python2.7/site-packages/astropy/io/fits/hdu/base.pyc in _readfrom_internal(cls, data, header, checksum, ignore_missing_end, **kwargs) 392 if header is None: 393 header_offset = data.tell() --> 394 header = Header.fromfile(data, endcard=not ignore_missing_end) 395 hdu_fileobj = data 396 data_offset = data.tell() # *after* reading the header /lustre/pipeline/scratch/GAS/python/anaconda/lib/python2.7/site-packages/astropy/io/fits/header.pyc in fromfile(cls, fileobj, sep, endcard, padding) 448 449 return cls._from_blocks(block_iter, is_binary, sep, endcard, --> 450 padding)[1] 451 finally: 452 if close_file: /lustre/pipeline/scratch/GAS/python/anaconda/lib/python2.7/site-packages/astropy/io/fits/header.pyc in _from_blocks(cls, block_iter, is_binary, sep, endcard, padding) 517 # TODO: Pass this error to validation framework as an ERROR, 518 # rather than raising an exception --> 519 raise IOError('Header missing END card.') 520 521 header_str = ''.join(read_blocks) IOError: Header missing END card.

[low-sky]

That looks like a readfits failure. Probably triggered by trying to read an invalid file at some point. I will try to reproduce on the GBT machines to isolate what's causing it.