Small Bugs

[jaym878]

I have a couple of bugs in my code right now and I was wondering if you could have a look.

1. I implemented your plotting function on the cuts to get the grid and the axes using ra and dec. The bug is that the selected coordinates earlier are not the same as the centre of the cuts. I cant really tell why its doing this.

2. I'm trying to reference the hdu from the image files to save the cut and stacked images as new files and to pass the data in the required format to do the aperture photometry.

   fluxes = []
   im_name = []
   for i in my_image_files:
   with fits.open(i) as hdulist:
       header = hdulist[0].header
       data = hdulist[0].data
       hdu = datasets.load_spitzer_image()
       wcs = WCS(header)
       position = SkyCoord(test_src_coord[0] * u.degree, test_src_coord[1] * u.degree, frame='icrs')
       apertures = SkyCircularAperture(position, r=Rad[SNR] * u.arcsec)
       phot_table = aperture_photometry(hdu, apertures)
   
       converted_aperture_sum = (phot_table['aperture_sum'] * arc_pix)
       fluxes.append(converted_aperture_sum)
       im_name.append(my_image_file(i))
   
       print(phot_table)

[PaddyKavanagh]

For the first one, the WCS object you are passing to plot_image is not the same as the WCS object in the cutout. The orginal WCS object is trimmed down by cutout2D and is in the .wcs attribute. Try:

plot_image(cutout.data, cutout.wcs, coords=test_src_coord)

For the second, is it that you are trying to recreate the format of hdu like in the example? I'm not sure of the exact format (what is the output of: print(hdulist.info) ) but you could try passing either hdulist or hdulist[0] in place if hdu.

[jaym878]

The aperture photometry code I have is working, but its returning the aperture sum as a NaN value. The required input is the primary HDU and the aperture shape. It does return the correct position and aperture shape so the problem must be with the handling of the data from the HDU. I have attached the code below.

for i in my_image_files:
    with fits.open(i) as hdulist:
        r = Rad[SNR] * u.arcsec

        position = SkyCoord(test_src_coord[0] * u.degree, test_src_coord[1] * u.degree, frame='icrs')
        apertures = SkyCircularAperture(position, r=Rad[SNR] * u.arcsec)

        phot_table = aperture_photometry(hdulist[0], apertures)

        print(phot_table) 
print(fluxes)

The returned data is in the format:

 id      xcenter       ... aperture_sum
           pix         ...   MJy / sr  
--- ------------------ ... ------------
  1 1561.7809690335544 ...          nan
[<Quantity [nan] MJy / sr>]

[PaddyKavanagh]

There may be some NaNs in the data array. I downloaded a test image and found that aperture_photometry was also returned nan for me. I changed the NaNs to zeroes and it works fine:

for i in my_image_files: with fits.open(i) as hdulist:

    my_hdu = hdulist[0]
    my_hdu.data = np.nan_to_num(my_hdu.data)

    r = Rad[SNR] * u.arcsec

    position = SkyCoord(test_src_coord[0] * u.degree, test_src_coord[1] * u.degree, frame='icrs')
    apertures = SkyCircularAperture(position, r=r)

    phot_table = aperture_photometry(my_hdu, apertures)

    print(phot_table) 

[PaddyKavanagh]

Sorry, the syntax highlighting is a bit off in the last comment:

for i in my_image_files:
    with fits.open(i) as hdulist:

        my_hdu = hdulist[0]
        my_hdu.data = np.nan_to_num(my_hdu.data)

        r = Rad[SNR] * u.arcsec

        position = SkyCoord(test_src_coord[0] * u.degree, test_src_coord[1] * u.degree, frame='icrs')
        apertures = SkyCircularAperture(position, r=r)

        phot_table = aperture_photometry(my_hdu, apertures)

        print(phot_table) 

[jaym878]

I've attached the flux calculation here. It might make it easier to read.

unit = flux * 23.5045
ap_ar = np.pi * Rad[SNR]**2
jan = unit * ap_ar
erg = jan * 10**-17
band_list = [3.6, 4.5, 5.8, 8.0, 24, 70, 160]
if wavelength == 3.6:
    band = 3.9-3.1
if wavelength == 4.5:
    band = 5.0-3.9
if wavelength == 5.8:
    band = 6.2-4.9
if wavelength == 8.0:
    band = 9.3-6.2
if wavelength == 24:
    band = 28-20
if wavelength == 70:
    band = 90-50
if wavelength == 160:
    band = 190-130
f = 300000000/band
fluxIR = erg * f
FluxIR = fluxIR[0]
print(FluxIR)

[PaddyKavanagh]

You should convert the band list from micron to Hz before subtracting one from the other.

Also, in the current form, there is a factor of 1e-6 missing in: f = 300000000/band to convert micron to metre.

[jaym878]

Current bug error

FileNotFoundError                         Traceback (most recent call last)
<ipython-input-1-4e176392aef4> in <module>()
     81 b = 8.0
     82 my_test = ImagePull(MCSNR, RA, DE, Rad, sensor, b)
---> 83 my_test.run()
     84 #for SNR in range(len(np.genfromtxt(open('SNR_list.csv', "r"), names=True, delimiter=',', dtype=None))):
     85 

/mnt/home_cr/murphyj/Desktop/Coding/astropull_pk.py in run(self)
    330 
    331         # Ap Phot Function
--> 332         flux = self.ap_phot(my_image_file, self.radius, test_src_coord, self.wavelength)
    333 
    334         print("Flux determined for {} in {}-{}:".format(self.name, self.sensor, self.wavelength))

/mnt/home_cr/murphyj/Desktop/Coding/astropull_pk.py in ap_phot(self, my_image_files, Rad, test_src_coord, wavelength)
    244 
    245             # load the file data, header, and wcs
--> 246             with fits.open(i) as hdulist:
    247                 my_hdu = hdulist[0]
    248                 my_hdu.data = np.nan_to_num(my_hdu.data)

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/hdu/hdulist.py in fitsopen(name, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)
    162 
    163     return HDUList.fromfile(name, mode, memmap, save_backup, cache,
--> 164                             lazy_load_hdus, **kwargs)
    165 
    166 

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/hdu/hdulist.py in fromfile(cls, fileobj, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)
    398         return cls._readfrom(fileobj=fileobj, mode=mode, memmap=memmap,
    399                              save_backup=save_backup, cache=cache,
--> 400                              lazy_load_hdus=lazy_load_hdus, **kwargs)
    401 
    402     @classmethod

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/hdu/hdulist.py in _readfrom(cls, fileobj, data, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)
    985             if not isinstance(fileobj, _File):
    986                 # instantiate a FITS file object (ffo)
--> 987                 fileobj = _File(fileobj, mode=mode, memmap=memmap, cache=cache)
    988             # The Astropy mode is determined by the _File initializer if the
    989             # supplied mode was None

~/anaconda3/lib/python3.6/site-packages/astropy/utils/decorators.py in wrapper(*args, **kwargs)
    486                         # one with the name of the new argument to the function
    487                         kwargs[new_name[i]] = value
--> 488             return function(*args, **kwargs)
    489 
    490         return wrapper

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/file.py in __init__(self, fileobj, mode, memmap, overwrite, cache)
    173             self._open_fileobj(fileobj, mode, overwrite)
    174         elif isinstance(fileobj, str):
--> 175             self._open_filename(fileobj, mode, overwrite)
    176         else:
    177             self._open_filelike(fileobj, mode, overwrite)

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/file.py in _open_filename(self, filename, mode, overwrite)
    526 
    527         if not self._try_read_compressed(self.name, magic, mode, ext=ext):
--> 528             self._file = fileobj_open(self.name, IO_FITS_MODES[mode])
    529             self.close_on_error = True
    530 

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/util.py in fileobj_open(filename, mode)
    386     """
    387 
--> 388     return open(filename, mode, buffering=0)
    389 
    390 

FileNotFoundError: [Errno 2] No such file or directory: 'J'

[jaym878]

I have an error when applying the rectangular aperture.

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-1-7d0ff050a681> in <module>()
    102             # Run image pull as object
    103             my_test = ImagePull(MCSNR, RA, DE, Rad, sensor, b)
--> 104             ret = my_test.run()
    105             FluxIR = ret[0]
    106             flux_lim = ret[1]

/mnt/home_cr/murphyj/Desktop/Coding/astropull.py in run(self)
    467         flux = self.ap_phot(my_image_file, self.radius, test_src_coord, self.wavelength)
    468 
--> 469         flux_lim = self.flux_arc(my_image_file, self.radius, test_src_coord, self.wavelength)
    470 
    471         return flux, flux_lim

/mnt/home_cr/murphyj/Desktop/Coding/astropull.py in flux_arc(self, my_image_files, Rad, test_src_coord, wavelength)
    399                 t4_pos = SkyCoord((test_src_coord[0] - y) * u.degree, (test_src_coord[1] - y) * u.degree, frame='icrs')
    400 
--> 401                 ap1 = SkyRectangularAperture((t1_pos, 60 * u.arcsec, 120 * u.arcsec, 0))
    402                 ap2 = SkyRectangularAperture((t1_pos, 60 * u.arcsec, 120 * u.arcsec, 0))
    403                 ap3 = SkyRectangularAperture((t1_pos, 60 * u.arcsec, 120 * u.arcsec, 0))

TypeError: __init__() missing 3 required positional arguments: 'w', 'h', and 'theta'

I define these values when calling the SkyRectangularAperture function but its asking to have them defined in the init() method. The code is uploaded on github and the error is in the flux arc method.

[PaddyKavanagh]

I think a pair of brackets need to be removed. Try:

ap1 = SkyRectangularAperture(t1_pos, 60 * u.arcsec, 120 * u.arcsec, theta=0.)

[jaym878]

I ran the program for all SNRs and it worked for the first 12 then crashed and gave this error. I cant workout why it crashed on that SNR

ValueError                                Traceback (most recent call last)
<ipython-input-1-9967dbbdaba9> in <module>()
     99             # Run image pull as object
    100                 my_test = ImagePull(MCSNR, RA, DE, Rad, sensor, b)
--> 101                 ret = my_test.run()
    102                 FluxIR = ret[0]
    103                 flux_high = ret[1]

/mnt/home_cr/murphyj/Desktop/Coding/astropull.py in run(self)
    601 
    602         # Query Function
--> 603         self.query(test_src_coord, self.sensor, self.wavelength)
    604 
    605         # Cuts Function

/mnt/home_cr/murphyj/Desktop/Coding/astropull.py in query(self, test_src_coord, sensor, wavelength)
    167                 distance_check[n] = sep.value
    168 
--> 169             closest = np.argmin(distance_check)
    170             my_final_row = tbl[closest]
    171             url = my_final_row['accessUrl'].strip()

~/anaconda3/lib/python3.6/site-packages/numpy/core/fromnumeric.py in argmin(a, axis, out)
   1066 
   1067     """
-> 1068     return _wrapfunc(a, 'argmin', axis=axis, out=out)
   1069 
   1070 

~/anaconda3/lib/python3.6/site-packages/numpy/core/fromnumeric.py in _wrapfunc(obj, method, *args, **kwds)
     50 def _wrapfunc(obj, method, *args, **kwds):
     51     try:
---> 52         return getattr(obj, method)(*args, **kwds)
     53 
     54     # An AttributeError occurs if the object does not have

ValueError: attempt to get argmin of an empty sequence

[PaddyKavanagh]

Looks to me like there is just no Spitzer data at that location which is causing the code to crash. Two options:

1. delete that SNR line from the catalogue

2. write in some check in astropull.query to skip the SNR if there is no data

Since we are a bit time constrained for your poster, I'd say just delete the row in the catalogue (and the first 12 successful rows) and run the rest.