astrodendro

import aplpy
from astropy.io import fits
import matplotlib.pyplot as plt
import numpy as np
import os
from astrodendro import Dendrogram
from astrodendro.analysis import PPStatistic
from astrodendro import Dendrogram, pp_catalog
from astropy import units as u
import numpy as np
from matplotlib.patches import Ellipse

inputfits='1.3mm_910klambda_subfits.fits'
outputfits='1.3mm_910klambda_subfits_aplpy.fits'
hdu = fits.open(inputfits)
hdr = hdu[0].header
hdr

def make_new_fits():
    #change header:
    f = fits.open(inputfits)
    hdr = f[0].header
    hdr.set('NAXIS',2)
    hdr.remove('NAXIS3')                                             
    hdr.remove('NAXIS4')                                                                          
    hdr.remove('PC1_1')                                             
    hdr.remove('PC2_1')                                             
    hdr.remove('PC3_1')                                             
    hdr.remove('PC4_1')                                             
    hdr.remove('PC1_2')                                             
    hdr.remove('PC2_2')                                             
    hdr.remove('PC3_2')                                             
    hdr.remove('PC4_2')                                             
    hdr.remove('PC1_3')                                             
    hdr.remove('PC2_3')                                             
    hdr.remove('PC3_3')                                             
    hdr.remove('PC4_3')                                             
    hdr.remove('PC1_4')                                             
    hdr.remove('PC2_4')                                             
    hdr.remove('PC3_4')                                             
    hdr.remove('PC4_4')                                                                                  
    hdr.remove('CTYPE3')                                             
    hdr.remove('CRVAL3')                                             
    hdr.remove('CDELT3')                                             
    hdr.remove('CRPIX3')                                             
    hdr.remove('CUNIT3')                                             
    hdr.remove('CTYPE4')                                             
    hdr.remove('CRVAL4')                                             
    hdr.remove('CDELT4')                                             
    hdr.remove('CRPIX4')                                             
    hdr.remove('CUNIT4')   
    #make new fits:
    hdu = fits.PrimaryHDU(f[0].data[0][0],header=hdr)
    path_exists = os.path.exists(outputfits)
    #print(path_exists)
    if path_exists==True:
        os.remove(outputfits)
        
    hdu.writeto(outputfits)
    
make_new_fits()

hdu = fits.open(outputfits)[0]
d = Dendrogram.compute(hdu.data, min_value=0.0009,min_delta = 0.00013,min_npix=15)
p = d.plotter()

fig = plt.figure(figsize=(60,20))
ax = fig.add_subplot(1,1,1)
ax.imshow(hdu.data,origin = 'lower',interpolation='nearest',cmap=plt.cm.jet,vmax=0.02)
for leaf in d.leaves:
    p.plot_contour(ax,astructure=leaf,lw=2,colors='red')

    s = PPStatistic(leaf)
    ellipse = s.to_mpl_ellipse(edgecolor = 'black',facecolor='none')

    ax.add_patch(ellipse)
plt.savefig('test.pdf')

d.leaves
metadata = {}
metadata['data_unit'] = u.Jy / u.beam
metadata['spatial_scale'] =  0.005 * u.arcsec
metadata['beam_major'] =  0.027 * u.arcsec
metadata['beam_minor'] =  0.022 * u.arcsec

cat = pp_catalog(d.leaves, metadata)
cat.pprint(show_unit=True, show_name=True, max_lines=70)
cat.write('table.csv', format='ascii.csv')