<commit_before>import pyspeckit
import astropy.io.fits as fits
import numpy as np
import os
from spectral_cube import SpectralCube
import signal_id
from radio_beam import Beam
import astropy.constants as con
import astropy.units as u
OneOneIntegrated = 'NGC1333_W11.fits'
OneOneFile = 'NGC1333_ABCDEFGH_NH3_11_all_blsub.rg.fits'
TwoTwoFile = 'NGC1333_ABCDEFGH_NH3_22_all_blsub.fits'
ThreeThreeFile = 'NGC1333_ABCDEFGH_NH3_33_all_blsub.fits'


beam11 = Beam.from_fits_header(fits.getheader(OneOneFile))
cube11sc = SpectralCube.read(OneOneFile)
cube22sc = SpectralCube.read(TwoTwoFile)

noisy = signal_id.noise.Noise(cube11sc)
noisy.estimate_noise(spectral_flat=True,niter=2)
errmap11 = (noisy.spatial_norm*noisy.scale)
rms = np.nanmedian(errmap11)
mask = (noisy.snr>5)*(noisy.spatial_norm<2.5)
maskcube = cube11sc.with_mask(mask)
maskcube = maskcube.with_spectral_unit(u.km/u.s,velocity_convention='radio')
slab = maskcube#.spectral_slab(5*u.km/u.s,15*u.km/u.s)
w11=slab.moment(0,0).value
peakloc = np.nanargmax(w11)
xmax,ymax = np.unravel_index(peakloc,w11.shape)
moment1 = slab.moment(1,0).value
moment2 = (slab.moment(2,0).value)**0.5
moment2[np.isnan(moment2)]=0.2
moment2[moment2<0.2]=0.2

maskmap = w11>0.5
cube11 = pyspeckit.Cube(OneOneFile,maskmap=maskmap)
cube11.units="K"
cube22 = pyspeckit.Cube(TwoTwoFile,maskmap=maskmap)
cube22.units="K"
cube33 = pyspeckit.Cube(ThreeThreeFile,maskmap=maskmap)
cube33.units="K"
cubes = pyspeckit.CubeStack([cube11,cube22,cube33],maskmap=maskmap)
cubes.units="K"

guesses = np.zeros((6,)+cubes.cube.shape[1:])
guesses[0,:,:] = 12                    # Kinetic temperature 
guesses[1,:,:] = 3                     # Excitation  Temp
guesses[2,:,:] = 14.5                  # log(column)
guesses[3,:,:] = moment2  # Line width / 5 (the NH3 moment overestimates linewidth)               guesses[4,:,:] = moment1  # Line centroid              
guesses[5,:,:] = 0.5                   # F(ortho) - ortho NH3 fraction (fixed)

F=False
T=True


cubes.fiteach(fittype='ammonia', multifit=True, guesses=guesses,
              integral=False, verbose_level=3, fixed=[F,F,F,F,F,T], signal_cut=3,
              limitedmax=[F,F,F,F,T,T],
              maxpars=[0,0,0,0,15,1],
              limitedmin=[T,T,F,T,T,T],
              minpars=[2.73,2.73,0,0,-5,0],
              use_nearest_as_guess=True, start_from_point=(xmax,ymax),
              multicore=1,
              errmap=errmap11)

fitcubefile = fits.PrimaryHDU(data=np.concatenate([cubes.parcube,cubes.errcube]), header=cubes.header)
fitcubefile.header.update('PLANE1','TKIN')
fitcubefile.header.update('PLANE2','TEX')
fitcubefile.header.update('PLANE3','COLUMN')
fitcubefile.header.update('PLANE4','SIGMA')
fitcubefile.header.update('PLANE5','VELOCITY')
fitcubefile.header.update('PLANE6','FORTHO')
fitcubefile.header.update('PLANE7','eTKIN')
fitcubefile.header.update('PLANE8','eTEX')
fitcubefile.header.update('PLANE9','eCOLUMN')
fitcubefile.header.update('PLANE10','eSIGMA')
fitcubefile.header.update('PLANE11','eVELOCITY')
fitcubefile.header.update('PLANE12','eFORTHO')
fitcubefile.header.update('CDELT3',1)
fitcubefile.header.update('CTYPE3','FITPAR')
fitcubefile.header.update('CRVAL3',0)
fitcubefile.header.update('CRPIX3',1)
fitcubefile.writeto("parameter_maps.fits",clobber=True)<commit_msg>Revised property mapping code. Requires changes to pyspeckit<commit_after>import pyspeckit
import astropy.io.fits as fits
import numpy as np
import os
from spectral_cube import SpectralCube
import signal_id
from radio_beam import Beam
import astropy.constants as con
import astropy.units as u
from skimage.morphology import remove_small_objects,closing,disk
OneOneIntegrated = 'NGC1333/NGC1333_NH3_11_mom0.fits'
OneOneFile = 'NGC1333/NGC1333_NH3_11_base1.fits'
RMSFile = 'NGC1333/NGC1333_NH3_11_base1_rms.fits'
TwoTwoFile = 'NGC1333/NGC1333_NH3_22_base1.fits'
ThreeThreeFile = 'NGC1333/NGC1333_NH3_33_base1.fits'

beam11 = Beam.from_fits_header(fits.getheader(OneOneFile))
cube11sc = SpectralCube.read(OneOneFile)
cube22sc = SpectralCube.read(TwoTwoFile)
errmap11 = fits.getdata(RMSFile)
snr = cube11sc.filled_data[:].value/errmap11
peaksnr = np.max(snr,axis=0)
rms = np.nanmedian(errmap11)
planemask = (peaksnr>3.5)*(errmap11 < 0.2)
planemask = remove_small_objects(planemask,min_size=20)
planemask = closing(planemask,disk(1))
#planemask = (peaksnr>20) * (errmap11 < 0.2)

mask = (snr>3)*planemask
maskcube = cube11sc.with_mask(mask.astype(bool))
maskcube = maskcube.with_spectral_unit(u.km/u.s,velocity_convention='radio')
slab = maskcube.spectral_slab(15*u.km/u.s,5*u.km/u.s)
w11=slab.moment(0,0).value
peakloc = np.nanargmax(w11)
ymax,xmax = np.unravel_index(peakloc,w11.shape)
moment1 = slab.moment(1,0).value
moment2 = (slab.moment(2,0).value)**0.5
moment2[np.isnan(moment2)]=0.2
moment2[moment2<0.2]=0.2
maskmap = w11>0.5
cube11 = pyspeckit.Cube(OneOneFile,maskmap=planemask)
cube11.units="K"
cube22 = pyspeckit.Cube(TwoTwoFile,maskmap=planemask)
cube22.units="K"
cube33 = pyspeckit.Cube(ThreeThreeFile,maskmap=planemask)
cube33.units="K"
cubes = pyspeckit.CubeStack([cube11,cube22,cube33],maskmap=planemask)
cubes.units="K"
guesses = np.zeros((6,)+cubes.cube.shape[1:])
guesses[0,:,:] = 12                    # Kinetic temperature 
guesses[1,:,:] = 3                     # Excitation  Temp
guesses[2,:,:] = 14.5                  # log(column)
guesses[3,:,:] = moment2  # Line width / 5 (the NH3 moment overestimates linewidth)               
guesses[4,:,:] = moment1  # Line centroid              
guesses[5,:,:] = 0.5                   # F(ortho) - ortho NH3 fraction (fixed)
F=False
T=True
print('start fit')
cubes.fiteach(fittype='ammonia',  guesses=guesses,
              integral=False, verbose_level=3, 
              fixed=[F,F,F,F,F,T], signal_cut=2,
              limitedmax=[F,F,F,F,T,T],
              maxpars=[0,0,0,0,15,1],
              limitedmin=[T,T,F,T,T,T],
              minpars=[5,2.8,0,0,5,0],
              start_from_point=(xmax,ymax),
              use_neighbor_as_guess=True, 
              position_order = 1/peaksnr,
              errmap=errmap11)

fitcubefile = fits.PrimaryHDU(data=np.concatenate([cubes.parcube,cubes.errcube]), header=cubes.header)
fitcubefile.header.update('PLANE1','TKIN')
fitcubefile.header.update('PLANE2','TEX')
fitcubefile.header.update('PLANE3','COLUMN')
fitcubefile.header.update('PLANE4','SIGMA')
fitcubefile.header.update('PLANE5','VELOCITY')
fitcubefile.header.update('PLANE6','FORTHO')
fitcubefile.header.update('PLANE7','eTKIN')
fitcubefile.header.update('PLANE8','eTEX')
fitcubefile.header.update('PLANE9','eCOLUMN')
fitcubefile.header.update('PLANE10','eSIGMA')
fitcubefile.header.update('PLANE11','eVELOCITY')
fitcubefile.header.update('PLANE12','eFORTHO')
fitcubefile.header.update('CDELT3',1)
fitcubefile.header.update('CTYPE3','FITPAR')
fitcubefile.header.update('CRVAL3',0)
fitcubefile.header.update('CRPIX3',1)
fitcubefile.writeto("parameter_maps.fits",clobber=True)
{'fraction_non_alphanumeric': 0.15178184966094357, 'fraction_numerical': 0.05295051219160294, 'mean_word_length': 5.161777777777778, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3262288', 'n_tokens_mistral': 3143, 'n_tokens_neox': 2858, 'n_words': 336}