""" @brief Second Series of skymap functions for GALPROP to sim EGRET datons @author A. Connors """ # @file second_simpler_skymap_datons.py # # #from numarray import * ## Changed to numpy Aug 2007 from numpy import * import pyfits ###import random class simpler_direct_interpol: def __init__(self, InBigSkyModel_hdu, TargetCoordsHeader): # Rather than doing anything at all that might resemble smoothing, # this module: # - takes in a bigskymap assumed in "GALPROP" Gamma-Ray units, # ( that is, MeV**2 /(cm**2 sr s MeV) ) # BOTH data and Input Coordinate system (hdu.data, hdu.header) # - TRANSFORMS to number-count of photons/input-bin [or something] # - takes in the coordinate system (from hdu.header) # from the wished-for model-map # ? Maybe it should check that energy-bins are consistent? # - THEN it makes a TEMPORARY 4xNbin binning of Target HDU-map. # - For each TempBin: # + it calculates the InputCoordinates of the bin corners # + it looks up the InputBigSkyMap value for each corner # + it does a quadrilateral sum, taking into account the # respective solid angles if Input and Target maps. # + It sums them up and puts them in the TargetMap Channel-bin. # - The output is an hdu in the target coordinate system. # # 1/ Initialize data format, based on input bigskymap format: # # 1.1/ Initialize header: # 1.1.1/ Transform Old Axis Descriptions Into New: # Warning! Doesn't always transform "cdelt" and "crval" properly!! ** old_naxis1 = bigskymap_hdu.header['NAXIS1'] naxis1 = old_naxis1/want_rebinfact[0] # old_cdelt1 = bigskymap_hdu.header['CDELT1'] # cdelt1 = old_cdelt1*want_rebinfact[0] # old_crval1 = bigskymap_hdu.header['CRVAL1'] # crval1 = old_crval1 + (want_rebinfact[0]*old_cdelt1/2.0000) # Use bin median (middle) as new crval old_naxis2 = bigskymap_hdu.header['NAXIS2'] naxis2 = old_naxis2/want_rebinfact[1] # old_cdelt2 = bigskymap_hdu.header['CDELT2'] # cdelt2 = old_cdelt2*want_rebinfact[1] # old_crval2 = bigskymap_hdu.header['CRVAL2'] # crval2 = old_crval2 + (want_rebinfact[1]*old_cdelt2/2.0000) # Use bin median (middle) as new crval # naxis3 = old_naxis3 = bigskymap_hdu.header['NAXIS3'] # crval3 = old_crval3 = bigskymap_hdu.header['CRVAL3'] # cdelt3 = old_cdelt3 = bigskymap_hdu.header['CDELT3'] naxis3 = 1 old_naxis3 = 1 data0 = zeros((naxis3,naxis2,naxis1),dtype=float ) # self = pyfits.PrimaryHDU(data0) ### PROBLEM WITH PASSING HDU FORMAT - KLUGE BELOW: self.data = data0 self.header = bigskymap_hdu.header # # 2/ Transform data from bigskymap units (flux) by summing bins print 'Debug: naxis1, naxis2, old,new naxis3:', naxis1, naxis2, old_naxis3, naxis3 print 'old,new:', bigskymap_hdu.data.shape, self.data.shape for l in range(old_naxis1): new_l = l/want_rebinfact[0] for b in range(old_naxis2): new_b = b/want_rebinfact[1] # print 'i: ',i,' l, new_l', l, new_l, ' b:',b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l], # print 'gp.data[which3rdbin][i][b][l]:',bigskymap_hdu.data[which3rdbin][i][b][l] self.data[0][new_b][new_l] = self.data[0][new_b][new_l] + bigskymap_hdu.data[which3rdbin][b][l] # # end loops over all image bins # # 3/ And finally, trying to update/clean out just the output header: # self.header.update('NAXIS1', naxis1) # self.header.update('CRVAL1', crval1) # self.header.update('CDELT1', cdelt1) self.header.update('NAXIS2', naxis2) # self.header.update('CRVAL2', crval2) # self.header.update('CDELT2', cdelt2) self.header.add_history('Rebinned.') # #------------------------------------------------------------------------------------------ # class simpler_average_rebin: def __init__(self, bigskymap_hdu,want_rebinfact,which_3rdbin): # # 1/ Initialize data format, based on input bigskymap format: # # 1.1/ Initialize header: # 1.1.1/ Transform Old Axis Descriptions Into New: # Warning! Doesn't always transform "cdelt" and "crval" properly!! ** old_naxis1 = bigskymap_hdu.header['NAXIS1'] naxis1 = old_naxis1/want_rebinfact[0] old_cdelt1 = bigskymap_hdu.header['CDELT1'] cdelt1 = old_cdelt1*want_rebinfact[0] old_crval1 = bigskymap_hdu.header['CRVAL1'] crval1 = old_crval1 + (want_rebinfact[0]*old_cdelt1/2.0000) # Use bin median (middle) as new crval old_naxis2 = bigskymap_hdu.header['NAXIS2'] naxis2 = old_naxis2/want_rebinfact[1] old_cdelt2 = bigskymap_hdu.header['CDELT2'] cdelt2 = old_cdelt2*want_rebinfact[1] old_crval2 = bigskymap_hdu.header['CRVAL2'] crval2 = old_crval2 + (want_rebinfact[1]*old_cdelt2/2.0000) # Use bin median (middle) as new crval naxis3 = 1 crval3 = 1. cdelt3 = 1. # 1.2/ Initialze Data data0 = zeros((naxis3,naxis2,naxis1),dtype=float ) # 1.3/ Turn it into fits "PrimaryHDU" format! # self = pyfits.PrimaryHDU(data0) ### PROBLEM WITH PASSING HDU FORMAT - KLUGE BELOW: self.data = data0 # 1.4/ Put in correct initial header: self.header = bigskymap_hdu.header # # 2/ Transform data from bigskymap units (flux) by summing bins print 'Debug: naxis1, naxis2, new naxis3:', naxis1, naxis2, naxis3 print 'old,new:', bigskymap_hdu.data.shape, self.data.shape newsize = 1.000000*(want_rebinfact[0]*want_rebinfact[1]) i = which_3rdbin for l in range(old_naxis1): new_l = l/want_rebinfact[0] for b in range(old_naxis2): new_b = b/want_rebinfact[1] # print 'i: ',i,' l, new_l', l, new_l, ' b:',b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l], # print 'gp.data[which3rdbin][i][b][l]:',bigskymap_hdu.data[i][b][l] self.data[0][new_b][new_l] = self.data[0][new_b][new_l] + (bigskymap_hdu.data[which_3rdbin][b][l]/newsize) # # end loops over all image bins # # 3/ And finally, trying to update/clean out just the output header: # self.header.update('NAXIS1', naxis1) self.header.update('CRVAL1', crval1) self.header.update('CDELT1', cdelt1) self.header.update('NAXIS2', naxis2) self.header.update('CRVAL2', crval2) self.header.update('CDELT2', cdelt2) self.header.update('NAXIS3', naxis3) self.header.add_history('Averageed.') # #------------------------------------------------------------------------------------------ # class simpler_outer_products: def __init__(self, modelflux_hdu, exposr_rebinned_hdu): # # 1/ Initialize data format, based on input modelflux format: # 1.1.0/ Get Axis Descriptions: naxis1 = modelflux_hdu.header['NAXIS1'] # cdelt1 = modelflux_hdu.header['CDELT1'] # crval1 = modelflux_hdu.header['CRVAL1'] naxis2 = modelflux_hdu.header['NAXIS2'] # cdelt2 = modelflux_hdu.header['CDELT2'] # crval2 = modelflux_hdu.header['CRVAL2'] exp_naxis1 = exposr_rebinned_hdu.header['NAXIS1'] # exp_cdelt1 = exposr_rebinned_hdu.header['CDELT1'] # exp_crval1 = modelflux_hdu.header['CRVAL1'] exp_naxis2 = exposr_rebinned_hdu.header['NAXIS2'] # exp_cdelt2 = exposr_rebinned_hdu.header['CDELT2'] # exp_crval2 = exposr_rebinned_hdu.header['CRVAL2'] naxis3 = 1 # crval3 = 1. # cdelt3 = 1. # # 1.1/ Do both have same naxis, cdelt? ASSUMED SO FOR NOW epsilon = 1.e-5 if naxis1 != exp_naxis1: print 'Error: Array Dims mis-match: data naxis1 ',naxis1,' is not equal to exposure naxis1:', exp_naxis1 return if naxis2 != exp_naxis2: print 'Error: Array Dims mis-match: data naxis2 ',naxis2,' is not equal to exposure naxis2:', exp_naxis2 return # if abs(cdelt1 - exp_cdelt1) > epsilon: # print 'Error: Array Dims mis-match: data cdelt1 ',cdelt1,' is not equal to exposure cdelt1:', exp_cdelt1 # return # if abs(cdelt2 - exp_cdelt2) > epsilon: # print 'Error: Array Dims mis-match: data cdelt2 ',cdelt2,' is not equal to exposure cdelt2:', exp_cdelt2 # return # if abs(crval1 - exp_crval1) > epsilon: # print 'Error: Array Dims mis-match: data crval1 ',crval1,' is not equal to exposure crval1:', exp_crval1 # return # if abs(crval2 - exp_crval2) > epsilon: # print 'Error: Array Dims mis-match: data crval2 ',crval2,' is not equal to exposure crval2:', exp_crval2 # return # # 1.2/ Initialize output data-model: data0 = zeros((naxis3,naxis2,naxis1),dtype=float ) self.data = data0 # # 1.3/ Put in fits "PrimaryHDU" format! # self = pyfits.PrimaryHDU(data0) ### PROBLEM WITH PASSING HDU FORMAT - KLUGE BELOW: # # 1.4/ Initialize output header: self.header = modelflux_hdu.header # # 2/ Transform data from modelflux units (flux) to model_data by outer-product: print 'Debug: naxis1, naxis2, new naxis3:', naxis1, naxis2, naxis3 print 'old,new:', modelflux_hdu.data.shape, self.data.shape print 'middle of old, new:',modelflux_hdu.data[0][naxis2/2][naxis1/2], print exposr_rebinned_hdu.data[0][naxis2/2][naxis1/2] print self.data[0][naxis2/2][naxis1/2] loprintlim, hiprintlim = (naxis2/2-10), (naxis2/2+10) for l in range(naxis1): for b in range(naxis2): exposure = exposr_rebinned_hdu.data[0][b][l] self.data[0][b][l] = modelflux_hdu.data[0][b][l]*exposure #if ( loprintlim < b and b < hiprintlim): # print 'l, b:', l,b # print ' modelflux, exposure:',modelflux_hdu.data[0][b][l],exposure # print ' self: ',self.data[0][b][l] #### self.data[0][b][l] = modelflux.data[0][b][l] # # end loops over all image bins # # 3/ And finally, trying to update/clean out just the output header: # self.header.add_history('Exposure*FluxModel*IntTime') # #------------------------------------------------------------------------------------------ # class simpler_twosum: def __init__(self, one_hdu, two_hdu): # # 1/ Initialize data format, based on input modelflux format: # 1.1.0/ Get Axis Descriptions: one_naxis1 = one_hdu.header['NAXIS1'] # one_cdelt1 = one_hdu.header['CDELT1'] # one_crval1 = one_hdu.header['CRVAL1'] one_naxis2 = one_hdu.header['NAXIS2'] # one_cdelt2 = one_hdu.header['CDELT2'] # one_crval2 = one_hdu.header['CRVAL2'] two_naxis1 = two_hdu.header['NAXIS1'] # two_cdelt1 = two_hdu.header['CDELT1'] # two_crval1 = two_hdu.header['CRVAL1'] two_naxis2 = two_hdu.header['NAXIS2'] # two_cdelt2 = two_hdu.header['CDELT2'] # two_crval2 = two_hdu.header['CRVAL2'] one_naxis3 = 1 two_naxis3 = 1 # crval3 = 1. # cdelt3 = 1. # # 1.1/ Do both have same naxis, cdelt? ASSUMED SO FOR NOW epsilon = 1.e-5 if one_naxis1 != two_naxis1: print 'Error: Array Dims mis-match: data naxis1 ',one_naxis1,' is not equal to exposure naxis1:', two_naxis1 return if one_naxis2 != two_naxis2: print 'Error: Array Dims mis-match: data naxis2 ',one_naxis2,' is not equal to exposure naxis2:', two_naxis2 return # if abs(cdelt1 - two_cdelt1) > epsilon: # print 'Error: Array Dims mis-match: data cdelt1 ',one_cdelt1,' is not equal to exposure cdelt1:', two_cdelt1 # return # if abs(cdelt2 - two_cdelt2) > epsilon: # print 'Error: Array Dims mis-match: data cdelt2 ',one_cdelt2,' is not equal to exposure cdelt2:', two_cdelt2 # return # if abs(crval1 - two_crval1) > epsilon: # print 'Error: Array Dims mis-match: data crval1 ',one_crval1,' is not equal to exposure crval1:', two_crval1 # return # if abs(crval2 - two_crval2) > epsilon: # print 'Error: Array Dims mis-match: data crval2 ',one_crval2,' is not equal to exposure crval2:', two_crval2 # return # # 1.2/ Initialize output data-model: data0 = zeros((one_naxis3,one_naxis2,one_naxis1),dtype=float ) self.data = data0 # # 1.3/ Put in fits "PrimaryHDU" format! # self = pyfits.PrimaryHDU(data0) ### PROBLEM WITH PASSING HDU FORMAT - KLUGE BELOW: # # 1.4/ Initialize output header: self.header = one_hdu.header print 'TwoSums check: input self header:', self.header # # 2/ Transform data print 'one_naxis1, one_naxis2, one_naxis3:', one_naxis1, one_naxis2, one_naxis3 print 'old,new:', one_hdu.data.shape, self.data.shape print 'middle of old, new:',one_hdu.data[0][one_naxis2/2][one_naxis1/2], print two_hdu.data[0][one_naxis2/2][one_naxis1/2] print self.data[0][one_naxis2/2][one_naxis1/2] loprintlim, hiprintlim = (one_naxis2/2-10), (one_naxis2/2+10) for l in range(one_naxis1): for b in range(one_naxis2): two = two_hdu.data[0][b][l] self.data[0][b][l] = one_hdu.data[0][b][l] + two #if ( loprintlim < b and b < hiprintlim): # print 'l, b:', l,b # print ' modelflux, exposure:',one_hdu.data[0][b][l],exposure # print ' self: ',self.data[0][b][l] #### self.data[0][b][l] = modelflux.data[0][b][l] # # end loops over all image bins # # 3/ And finally, trying to update/clean out just the output header: # self.header.add_history('Exposure*FluxModel*IntTime') # #------------------------------------------------------------------------------------------ # class simpler_poisson_datons: def __init__(self, modeldata_hdu): # # 1/ Initialize data format, based on input modeldata format: # 1.0/ Do both have same naxis, cdelt? ASSUMED SO FOR NOW # 1.1/ Transform Old Axis Descriptions Into New: naxis1 = modeldata_hdu.header['NAXIS1'] # cdelt1 = modeldata_hdu.header['CDELT1'] # crval1 = modeldata_hdu.header['CRVAL1'] naxis2 = modeldata_hdu.header['NAXIS2'] # cdelt2 = modeldata_hdu.header['CDELT2'] # crval2 = modeldata_hdu.header['CRVAL2'] naxis3 = 1 # crval3 = 1. # cdelt3 = 1. # # 1.2/ Initialize data data0 = zeros((naxis3,naxis2,naxis1),dtype=float ) # # 1.3/ Fits format: Primary HDU! # self = pyfits.PrimaryHDU(data0) ### PROBLEM WITH PASSING HDU FORMAT - KLUGE BELOW: self.data = data0 # # 1.4/ Initialize header: self.header = modeldata_hdu.header # # 2/ Transform data from modeldata units (flux) to model_data by outer-product: print 'Debug: naxis1, naxis2, new naxis3:', naxis1, naxis2, naxis3 print 'old,new:', modeldata_hdu.data.shape, self.data.shape for l in range(naxis1): for b in range(naxis2): # print ' l', l, ' b:',b # print 'self.data[0][b][l]:',self.data[0][b][l], # print 'modeldata_hdu.data[0][b][l]:',modeldata_hdu.data[0][b][l] mu = modeldata_hdu.data[0][b][l] if mu > 0: self.data[0][b][l] = random.poisson(mu) else: self.data[0][b][l] = 0. # # end loops over all image bins # # 3/ And finally, trying to update/clean out just the output header: # self.header.add_history('GammaVariate of DataModel') # #------------------------------------------------------------------------------------------ # class simpler_sphereskymap_pad_to_square_power2: def __init__(self, modeldata_hdu): # PowersOf2 = [1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,32768] # **Needs exception for data bigger than len[]!! # 1/ Initialize data format, based on input modeldata format: # # 1.1/ Transform Old Axis Descriptions Into New: naxis1 = modeldata_hdu.header['NAXIS1'] naxis2 = modeldata_hdu.header['NAXIS2'] naxis3 = 1 crval3 = 1. cdelt3 = 1. naxis_big = max(naxis1,naxis2) for i2 in PowersOf2: if i2/2 <= naxis_big < i2 : print naxis_big, i2 naxissq = i2 # 1.2/ Initialize data ---- with new "naxs1"="naxis2"=naxissq: data0 = zeros((naxis3,naxissq,naxissq),dtype=float) # 1.3/ Put into Fits format: Primary HDU! # self = pyfits.PrimaryHDU(data0) ### PROBLEM WITH PASSING HDU FORMAT - KLUGE BELOW: self.data = data0 # 1.4/ Initialize header: self.header = modeldata_hdu.header print 'Original naxis1, naxis2=',naxis1,naxis2,' to be padded to new square naxissq:', naxissq # # 2/ Transform data from modeldata units (flux) to model_data by outer-product: # # 2.1/ Find the new padding amounts on Gal-l minus; gal-l plus; gal-b bottom(-); gal-b top(+): delpad_l_m = (naxissq - naxis1)/2 delpad_l_p = naxissq - naxis1 - delpad_l_m delpad_b_m = (naxissq - naxis2)/2 delpad_b_p = naxissq - naxis2 - delpad_b_m print 'new padding on Gal-l minus,plus:',delpad_l_m, delpad_l_p, print ' Gal-b minus, plus:',delpad_b_m, delpad_b_p, print 'old,new:', modeldata_hdu.data.shape, self.data.shape # 2.2/ Fill in the middle "square" with the original skymap: for l in range(naxis1): new_l = l + delpad_l_m for b in range(naxis2): new_b = b + delpad_b_m # print ' l, new_l', l, ' b, new_b:',b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l], # print 'modeldata_hdu.data[0][b][l]:',modeldata_hdu.data[0][b][l] self.data[0][new_b][new_l] = modeldata_hdu.data[0][b][l] # 2.3/ Fill in the horizontal (i.e. Gal-l) part 1st, as it has simpler "cycling" rule: # 2.3.1/ Fill in the "minus" side of Gal-l: for new_l in range(delpad_l_m): old_l = naxis1 - delpad_l_m + new_l for b in range(naxis2): new_b = b + delpad_b_m # print ' old_l, new_l', old_l, ' b, new_b:',b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l], # print 'modeldata_hdu.data[0][b][old_l]:',modeldata_hdu.data[0][b][old_l] self.data[0][new_b][new_l] = modeldata_hdu.data[0][b][old_l] # 2.3.2/ Fill in the "plus" side of Gal-l: for del_l in range(delpad_l_p): new_l = del_l + naxis1 + delpad_l_m old_l = del_l for b in range(naxis2): new_b = b + delpad_b_m # print ' old_l, new_l', old_l, new_l,' b, new_b:',b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l], # print 'modeldata_hdu.data[0][b][old_l]:',modeldata_hdu.data[0][b][old_l] self.data[0][new_b][new_l] = modeldata_hdu.data[0][b][old_l] # 2.4/ Fill in the Vertical (i.e. Gal-b) part 2nd, as it has more complex "cycling" rule: # # | /---------\ # | / /---\ \ | # a18 | a11 a12 a13 a14 a15 a16 a17 a18 | a11 # a28 | a21 a22 a23 a24 a25 a26 a27 a28 | a21 # a38 | a31 a32 a33 a34 a35 a36 a37 a38 | a31 # a48 | a41 a42 a43 a44 a45 a46 a47 a48 | a41 # | \ \---/ / | # | \---------/ | # # 2.4.1/ Fill in the "minus" side of Gal-b: for old_l in range(naxissq): new_l = naxissq - old_l - 1 for new_b in range(delpad_b_m): old_b = 2*delpad_b_m - new_b - 1 # print ' new_l', new_l ' old_b, new_b:', old_b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l] self.data[0][new_b][new_l] = self.data[0][old_b][old_l] # # 2.4.2/ Fill in the "plus" side of Gal-b: for old_l in range(naxissq): new_l = naxissq - old_l - 1 for del_b in range(delpad_b_p): new_b = delpad_b_m + naxis2 + del_b old_b = delpad_b_m + naxis2 - 1 - del_b # print ' new_l', new_l ' old_b, new_b:', old_b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l] self.data[0][new_b][new_l] = self.data[0][old_b][old_l] # end loops over all image bins # # 3/ And finally, trying to update/clean out just the output header: # self.header.update('NAXIS1', naxissq) self.header.update('NAXIS2', naxissq) self.header.add_history('Sphere-projection skymap padded to square power of 2') # #------------------------------------------------------------------------------------------ # class simpler_zerozero_pad_to_square_power2: def __init__(self, modeldata_hdu): # PowersOf2 = [1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,32768] # **Needs exception for data bigger than len[]!! # 1/ Initialize data format, based on input modeldata format: # # 1.1/ Transform Old Axis Descriptions Into New: naxis1 = modeldata_hdu.header['NAXIS1'] naxis2 = modeldata_hdu.header['NAXIS2'] naxis3 = 1 crval3 = 1. cdelt3 = 1. naxis_big = max(naxis1,naxis2) for i2 in PowersOf2: if i2/2 <= naxis_big < i2 : print naxis_big, i2 naxissq = i2 # 1.2/ Initialize data ---- with new "naxs1"="naxis2"=naxissq: data0 = zeros((naxis3,naxissq,naxissq),dtype=float) # 1.3/ Put into Fits format: Primary HDU! # self = pyfits.PrimaryHDU(data0) ### PROBLEM WITH PASSING HDU FORMAT - KLUGE BELOW: self.data = data0 # 1.4/ Initialize header: self.header = modeldata_hdu.header print 'Original naxis1, naxis2=',naxis1,naxis2,' to be padded to new square naxissq:', naxissq # # 2/ Transform data from modeldata units (flux) to model_data by outer-product: # # 2.1/ Find the new padding amounts on Gal-l minus; gal-l plus; gal-b bottom(-); gal-b top(+): delpad_l_m = (naxissq - naxis1)/2 delpad_l_p = naxissq - naxis1 - delpad_l_m delpad_b_m = (naxissq - naxis2)/2 delpad_b_p = naxissq - naxis2 - delpad_b_m print 'new padding on Gal-l minus,plus:',delpad_l_m, delpad_l_p, print ' Gal-b minus, plus:',delpad_b_m, delpad_b_p, print 'old,new:', modeldata_hdu.data.shape, self.data.shape # 2.2/ Fill in the middle "square" with the original skymap: for l in range(naxis1): new_l = l + delpad_l_m for b in range(naxis2): new_b = b + delpad_b_m # print ' l, new_l', l, ' b, new_b:',b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l], # print 'modeldata_hdu.data[0][b][l]:',modeldata_hdu.data[0][b][l] self.data[0][new_b][new_l] = modeldata_hdu.data[0][b][l] # 2.3/ Fill in the horizontal (i.e. Gal-l) part 1st, as it has simpler "cycling" rule: # 2.3.1/ Fill in the "minus" side of Gal-l: for new_l in range(delpad_l_m): old_l = naxis1 - delpad_l_m + new_l for b in range(naxis2): new_b = b + delpad_b_m # print ' old_l, new_l', old_l, ' b, new_b:',b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l], # print 'modeldata_hdu.data[0][b][old_l]:',modeldata_hdu.data[0][b][old_l] self.data[0][new_b][new_l] = modeldata_hdu.data[0][b][old_l] # 2.3.2/ Fill in the "plus" side of Gal-l: for del_l in range(delpad_l_p): new_l = del_l + naxis1 + delpad_l_m old_l = del_l for b in range(naxis2): new_b = b + delpad_b_m # print ' old_l, new_l', old_l, new_l,' b, new_b:',b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l], # print 'modeldata_hdu.data[0][b][old_l]:',modeldata_hdu.data[0][b][old_l] self.data[0][new_b][new_l] = modeldata_hdu.data[0][b][old_l] # 2.4/ Fill in the Vertical (i.e. Gal-b) part with zeros or min, for when cycling isn't useful # # 2.4.1/ Fill in the "minus" side of Gal-b: for old_l in range(naxissq): new_l = naxissq - old_l - 1 for new_b in range(delpad_b_m): old_b = 2*delpad_b_m - new_b - 1 # print ' new_l', new_l ' old_b, new_b:', old_b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l] self.data[0][new_b][new_l] = 1.e-26 # # 2.4.2/ Fill in the "plus" side of Gal-b: for old_l in range(naxissq): new_l = naxissq - old_l - 1 for del_b in range(delpad_b_p): new_b = delpad_b_m + naxis2 + del_b old_b = delpad_b_m + naxis2 - 1 - del_b # print ' new_l', new_l ' old_b, new_b:', old_b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l] self.data[0][new_b][new_l] = 1.e-26 # end loops over all image bins # # 3/ And finally, trying to update/clean out just the output header: # self.header.update('NAXIS1', naxissq) self.header.update('NAXIS2', naxissq) self.header.add_history('Sphere-projection skymap padded to square power of 2') # #------------------------------------------------------------------------------------------ # class simpler_zero_pad_to_square_power2: def __init__(self, modeldata_hdu): # PowersOf2 = [1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,32768] # **Needs exception for data bigger than len[]!! # 1/ Initialize data format, based on input modeldata format: # # 1.1/ Transform Old Axis Descriptions Into New: naxis1 = modeldata_hdu.header['NAXIS1'] naxis2 = modeldata_hdu.header['NAXIS2'] naxis3 = 1 crval3 = 1. cdelt3 = 1. naxis_big = max(naxis1,naxis2) for i2 in PowersOf2: if i2/2 <= naxis_big < i2 : print naxis_big, i2 naxissq = i2 # 1.2/ Initialize data ---- with new "naxs1"="naxis2"=naxissq: data0 = zeros((naxis3,naxissq,naxissq),dtype=float) # 1.3/ Put into Fits format: Primary HDU! # self = pyfits.PrimaryHDU(data0) ### PROBLEM WITH PASSING HDU FORMAT - KLUGE BELOW: self.data = data0 # 1.4/ Initialize header: self.header = modeldata_hdu.header print 'Original naxis1, naxis2=',naxis1,naxis2,' to be padded to new square naxissq:', naxissq # # 2/ Transform data from modeldata units (flux) to model_data by outer-product: # # 2.1/ Find the new padding amounts on Gal-l minus; gal-l plus; gal-b bottom(-); gal-b top(+): delpad_l_m = (naxissq - naxis1)/2 delpad_l_p = naxissq - naxis1 - delpad_l_m delpad_b_m = (naxissq - naxis2)/2 delpad_b_p = naxissq - naxis2 - delpad_b_m print 'new padding on Gal-l minus,plus:',delpad_l_m, delpad_l_p, print ' Gal-b minus, plus:',delpad_b_m, delpad_b_p, print 'old,new:', modeldata_hdu.data.shape, self.data.shape # 2.2/ Fill in the middle "square" with the original skymap: for l in range(naxis1): new_l = l + delpad_l_m for b in range(naxis2): new_b = b + delpad_b_m # print ' l, new_l', l, ' b, new_b:',b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l], # print 'modeldata_hdu.data[0][b][l]:',modeldata_hdu.data[0][b][l] self.data[0][new_b][new_l] = modeldata_hdu.data[0][b][l] # 2.3/ Fill in the horizontal (i.e. Gal-l) part 1st, as it has simpler "cycling" rule: # 2.3.1/ Fill in the "minus" side of Gal-l: for new_l in range(delpad_l_m): old_l = naxis1 - delpad_l_m + new_l for b in range(naxis2): new_b = b + delpad_b_m # print ' old_l, new_l', old_l, ' b, new_b:',b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l], # print 'modeldata_hdu.data[0][b][old_l]:',modeldata_hdu.data[0][b][old_l] self.data[0][new_b][new_l] = modeldata_hdu.data[0][b][old_l] # 2.3.2/ Fill in the "plus" side of Gal-l: for del_l in range(delpad_l_p): new_l = del_l + naxis1 + delpad_l_m old_l = del_l for b in range(naxis2): new_b = b + delpad_b_m # print ' old_l, new_l', old_l, new_l,' b, new_b:',b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l], # print 'modeldata_hdu.data[0][b][old_l]:',modeldata_hdu.data[0][b][old_l] self.data[0][new_b][new_l] = modeldata_hdu.data[0][b][old_l] # 2.4/ Fill in the Vertical (i.e. Gal-b) part with zeros or min, for when cycling isn't useful # # 2.4.1/ Fill in the "minus" side of Gal-b: for old_l in range(naxissq): new_l = naxissq - old_l - 1 for new_b in range(delpad_b_m): old_b = 2*delpad_b_m - new_b - 1 # print ' new_l', new_l ' old_b, new_b:', old_b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l] self.data[0][new_b][new_l] = 1.e-13 # # 2.4.2/ Fill in the "plus" side of Gal-b: for old_l in range(naxissq): new_l = naxissq - old_l - 1 for del_b in range(delpad_b_p): new_b = delpad_b_m + naxis2 + del_b old_b = delpad_b_m + naxis2 - 1 - del_b # print ' new_l', new_l ' old_b, new_b:', old_b, new_b # print 'self.data[0][new_b][new_l]:',self.data[0][new_b][new_l] self.data[0][new_b][new_l] = 1.e-13 # end loops over all image bins # # 3/ And finally, trying to update/clean out just the output header: # self.header.update('NAXIS1', naxissq) self.header.update('NAXIS2', naxissq) self.header.add_history('Sphere-projection skymap padded to square power of 2') # #------------------------------------------------------------------------------------------ #