Bug fixes for AMUSING processing

.gitignore

@@ -46,3 +46,4 @@ demo_notebk/wong2024/bpt_sncut3.pdf
 demo_notebk/wong2024/bpt_sncut3.png
 demo_notebk/wong2024/bpt_sncut3_zoh.pdf
 demo_notebk/wong2024/bpt_sncut3_zoh.png
+edge_pydb/img_comom/aquest_comom_2p5

bin/edge_db_gen_ifu_large

@@ -55,7 +55,7 @@ do_pipe3d(gallist=amusing, outfile=outfile,
           ssptable=pipe3d_dir+'gsd01_156.fits', nsm=2,
           ortpar='/Users/tonywong/Scratch2/amusing/tables/amusing_leda.csv', 
           distpar='/Users/tonywong/Scratch2/amusing/tables/amusing_global.csv', 
-          coln_dmpc='DistP3d', prob=False, 
+          coln_dmpc='DistP3d', prob=False, blankval=-1,
           prodtype=['SFH', 'SSP', 'indices', 'flux_elines'],
           leadstr=['', '', 'indices.CS.', 'flux_elines.'],
           tailstr=['.SFH','.SSP','',''])

bin/edge_db_gen_ifugrid

@@ -58,3 +58,25 @@ do_comom(gallist=aca60, outfile=outfile,
          alphaco=6.6, fitsdir=comom_dir+'aca12', allpix=True, append=True,
          p3d_dir=pipe3d_dir+'fits_natv_aca')
 
+# ALMaQUEST galaxies, MaNGA native resolution processing smoothed to 2.5 kpc CO
+outfile = output_dir + 'almaquest_allmatch.2p5kpc.hdf5'
+if os.path.isfile(outfile): os.remove(outfile)
+aquest = [os.path.basename(file).split('.')[0] for file in sorted(
+          glob.glob(comom_dir+'aquest_comom_2p5/*_dil.snrpk.fits.gz'))]
+aquest2 = [os.path.basename(file).split('.')[0].split('-',1)[1] for file in sorted(
+          glob.glob(pipe3d_dir+'fits_natv_aq/*.Pipe3D.cube.fits.gz'))]
+do_pipe3d(gallist=aquest2, outfile=outfile,
+          fitsdir=pipe3d_dir+'fits_natv_aq', p3dstruct='manga', packed=True,
+          matchres=True, ext='', comomdir=comom_dir+'aquest_comom_2p5', 
+          cotempl='manga_GNAME.co.kpc2p5_dil.snrpk.fits.gz',
+          ssptable='/Users/tonywong/Work/projects/MaNGA/mastar_237.fits',
+          ortpar='/Users/tonywong/Work/projects/MaNGA/MaNGA_props_pipe3d.csv', 
+          distpar='/Users/tonywong/Work/projects/MaNGA/MaNGA_props_pipe3d.csv', 
+          coln_ra='objra', coln_dc='objdec', coln_pa='nsa_sersic_phi',
+          coln_inc='nsa_inclination', coln_dmpc='nsa_z_dMpc', allpix=True)
+do_comom(gallist=aquest, outfile=outfile, fitsdir=comom_dir+'aquest_comom_2p5', 
+         linelbl=['co'], lines=['12'], seq='kpc2p5',
+         ortpar='/Users/tonywong/Work/projects/MaNGA/MaNGA_props_pipe3d.csv', 
+         ortlabel='p3d', coln_ra='objra', coln_dc='objdec', coln_pa='nsa_sersic_phi', 
+         coln_inc='nsa_inclination', append=True, manganame=True, allpix=True,
+         p3d_dir=pipe3d_dir+'fits_natv_aq', p3dtempl='GNAME.Pipe3D.cube.fits.gz')

edge_pydb/fitsextract.py

@@ -297,14 +297,14 @@ def getlabels(product, p3dstruct='califa'):
         # We select the bright lines that are also in ELINES, plus [OI]6300
         has_errors = True
         flux  = [0, 26, 27, 28, 41, 45, 46, 47, 49, 50]
-        nline = len(flux)
         if p3dstruct == 'califa':
             nz = 408
         elif p3dstruct == 'manga':
             nz = 456
         elif p3dstruct == 'amusing':
             nz = 240
             flux  = [1, 2, 3, 19, 20, 21, 22, 24, 25]
+        nline = len(flux)
         nfelines = nz // 8
         vel   = list(np.array(flux)+nfelines)
         disp  = list(np.array(flux)+nfelines*2)

edge_pydb/img_califa/do_pipe3d.py

@@ -32,7 +32,7 @@ def do_pipe3d(outfile='NGC4047.pipe3d.hdf5', gallist=['NGC4047'], fitsdir=None,
               distpar='edge_califa.csv', ortlabel='LEDA', coln_ra='ledaRA', 
               coln_dc='ledaDE', coln_pa='ledaPA', coln_inc='ledaAxIncl',
               coln_dmpc='caDistP3d', hexgrid=False, allpix=False, debug=False, 
-              keepnan=True, blankzero=True, prob=True, discard_cdmatrix=False, 
+              keepnan=True, blankval=0, prob=True, discard_cdmatrix=False, 
               overwrite=True, matchres=False,
               prodtype=['ELINES', 'SFH', 'SSP', 'indices', 'flux_elines'],
               leadstr=['', '', '', 'indices.CS.', 'flux_elines.'],
@@ -109,9 +109,11 @@ def do_pipe3d(outfile='NGC4047.pipe3d.hdf5', gallist=['NGC4047'], fitsdir=None,
     debug : boolean
         True to generate some additional output
     keepnan : boolean
-        True to replace zeroes in the FITS file with NaN.  This is the default.
-    blankzero : boolean
-        True to replace zeroes in the FITS file with NaN.  This is the default.
+        If False, the output table drops rows which are all-NaN.
+        Default is True (keep the NaNs).
+    blankval : int
+        Integer value in the FITS file to replace with NaN.  Default is 0.
+        Set blankval=None to omit any additional blanking.
     prob : boolean
         True to add BPT probability column to flux_elines table.
     discard_cdmatrix : boolean
@@ -251,9 +253,9 @@ def do_pipe3d(outfile='NGC4047.pipe3d.hdf5', gallist=['NGC4047'], fitsdir=None,
                                             y_size=copsf.shape[0])
                 elif p3dstruct == 'manga':
                     if 'RFWHM' in p3dhd.keys():
-                        target_fwhm = hdr['RFWHM']/pixsca.value  # convert to pixels
+                        fwhm = p3dhd['RFWHM']/pixsca.value  # convert to pixels
                     else:
-                        target_fwhm = 2.54/pixsca.value  # MaNGA sample median, (Yan+16)
+                        fwhm = 2.54/pixsca.value  # MaNGA sample median, (Yan+16)
                     print('Gaussian profile parameters: fwhm={}'.format(fwhm))
                     ifupsf = Gaussian2DKernel(fwhm/np.sqrt(8*np.log(2)),
                                      x_size=copsf.shape[1], y_size=copsf.shape[0])
@@ -317,8 +319,8 @@ def do_pipe3d(outfile='NGC4047.pipe3d.hdf5', gallist=['NGC4047'], fitsdir=None,
             desckeys = [key for key in list(cahd.keys()) if key.startswith('DESC')]
             for key in desckeys:
                 w_cahd[key] = cahd[key]
-            if blankzero:
-                cadat[cadat==0] = np.nan
+            if blankval is not None:
+                cadat[cadat==blankval] = np.nan
 
             # Regrid to the CO template
             if not matchres and comomdir is not None:

edge_pydb/img_comom/do_comom.py

@@ -134,12 +134,24 @@ def do_comom(outfile='NGC4047.2d_smo7.hdf5', gallist=['NGC4047'], seq='smo7',
             newhd = hdul[0].header.copy()
             if p3d_dir is not None:
                 p3dfile = os.path.join(p3d_dir,p3dtempl.replace('GNAME',gal))
+                if not os.path.exists(p3dfile):
+                    # Fudge for almaquest file naming convention
+                    gal2 = gal.replace('_','-')
+                    p3dfile = os.path.join(p3d_dir,p3dtempl.replace('GNAME',gal2))
+                    if not os.path.exists(p3dfile):
+                        print('####### Cannot find',p3dfile)
+                        continue
                 p3dhd = fits.getheader(p3dfile)
                 hd2d = WCS(p3dhd).celestial.to_header()
                 for key in hd2d.keys():
                     newhd[key] = hd2d[key]
+                # For packed files NAXIS1 and NAXIS2 may be missing
                 for key in ['NAXIS1', 'NAXIS2']:
-                    newhd[key] = p3dhd[key]
+                    if key in p3dhd.keys():
+                        newhd[key] = p3dhd[key]
+                    else:
+                        p3dhd1 = fits.getheader(p3dfile, 1)
+                        newhd[key] = p3dhd1[key]
                 newim = reproject_interp(hdul[0], newhd, order=interp_order,
                                          return_footprint=False)
 #                 if debug:

edge_pydb/plotting.py

@@ -343,7 +343,7 @@ def gridplot(edgetab=None, gallist=None, columnlist=None,
         # Plot all galaxies by default
         if gallist is None:
             gallist = list(np.unique(edgetab['Name']))
-        print('Plotting column',columnlist[0],'for',len(gallist),'galaxies')
+        print('\nPlotting column',columnlist[0],'for',len(gallist),'galaxies')
         pagelist = gallist
         if allnorm:
             if not np.isnan(edgetab[columnlist[0]]).all():