Python WLanalysis.coords2grid方法代码示例

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import coords2grid [as 别名]
def partialdata2grid (icount):
    '''for a small portion of the data, put into a grid'''
    print 'icount',icount
    idata = load(mask_dir+'smaller/weight0_cat_W%i_step%i_start%i.npy'%(Wx,step, icount))
    print 'loaded',icount
    ix, iy=np.indices(idata.shape)
    iy+=step*icount
    #radeclist = (array(w.wcs_pix2world(ix, iy, 0)).reshape(2,-1)).T ////jia changed on 12/9, since the coordinates seems to be off..
    #radeclist = (array(w.wcs_pix2world(iy, ix, 0)).reshape(2,-1)).T 
    #y, x = f_Wx (radeclist)
    y, x = f_Wx ((array(w.wcs_pix2world(iy, ix, 0)).reshape(2,-1)).T )
    print icount,'done f_wx %s'%(icount)#,time.strftime("%Y-%m-%d %H:%M")

    if Wx!=1:
        ipix_mask, ipix = WLanalysis.coords2grid(x,y, idata.flatten().reshape(1,-1), size=sizes[Wx-1])
        
    else:
        istep=ceil(len(y)/10.0)
        ipix_mask = zeros(shape=(sizes[Wx-1], sizes[Wx-1]))
        ipix = ipix_mask.copy()
        jj=0    
        while jj< len(y):
            print 'icount, jj',icount,jj
            iipix_mask,iipix = WLanalysis.coords2grid(x[jj:jj+istep], y[jj:jj+istep], idata.flatten().reshape(1,-1)[:,jj:jj+istep], size=sizes[Wx-1])
            ipix_mask += iipix_mask
            ipix += iipix
            jj+=istep
    print icount,'W%i done coords2grid %s'%(Wx,icount)#,time.strftime("%Y-%m-%d %H:%M")    
    
    save(mask_dir+'smaller/weight0_W%i_%i_numpix'%(Wx,icount), ipix)
    save(mask_dir+'smaller/weight0_W%i_%i_nummask'%(Wx,icount), ipix_mask)
    #ipix is the num. of pixels fall in that big pix, ipix_mask is the mask
    return ipix, ipix_mask

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import coords2grid [as 别名]
def kmapPk_1sim (r):
	print r,'1sim'
	k, s1, s2 = kappaGen_1sim(r)[:3]
	A, galn = WLanalysis.coords2grid(x, y, array([s1,s2]))
	Ms1,Ms2 = A
	s1_smooth = WLanalysis.weighted_smooth(Ms1, galn, PPA=PPA512, sigmaG=sigmaG)
	s2_smooth = WLanalysis.weighted_smooth(Ms1, galn, PPA=PPA512, sigmaG=sigmaG)
	kmap = WLanalysis.KSvw(s1_smooth, s2_smooth)
	pk = WLanalysis.peaks_mask_hist(kmap, mask, bins=25, kmin = -0.04, kmax = 0.12)
	return pk

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import coords2grid [as 别名]
def kmapPs (r):
	print r
	k, s1, s2 = kappaGen(r)[:3]
	A, galn = WLanalysis.coords2grid(x, y, array([s1,s2 ]))
	Ms1,Ms2 = A
	s1_smooth = WLanalysis.weighted_smooth(Ms1, galn, PPA=PPA512, sigmaG=sigmaG)
	s2_smooth = WLanalysis.weighted_smooth(Ms1, galn, PPA=PPA512, sigmaG=sigmaG)
	kmap = WLanalysis.KSvw(s1_smooth, s2_smooth)
	kmap*=mask
	ps = WLanalysis.PowerSpectrum(kmap,sizedeg=12.0)[-1]
	return ps

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import coords2grid [as 别名]
def fileGen(i, R, cosmo):
	'''
	Put catalogue to grid, with (1+m)w correction. Mw is already done.
	also add randomly rotated noise
	Input:
	i: subfield range from (1, 2..13)
	R: realization range from (1..1000)
	cosmo: one of the 100 cosmos
	Return:
	Me1 = e1*w
	Me2 = e2*w
	
	'''
	#y, x, e1, e2, w, m = yxewm_arr[i-1].T
	s1, s2 = (WLanalysis.readFits(SIMfn(i,cosmo,R)).T)[[1,2]]
	A, galn = WLanalysis.coords2grid(x, y, array([s1*w, s2*w]))
	Ms1, Ms2 = A
	return Ms1, Ms2

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import coords2grid [as 别名]
def createBadFieldMask (sf):
	sf_splitfiles = os.listdir(sf_dir(sf))
	genfromtxtA = lambda fn: genfromtxt(sf_dir(sf)+fn)
	datas = map(genfromtxtA,sf_splitfiles)#3 columns: RA, DEC, GB
	datas = concatenate(datas,axis=0)
	idx = where(datas[:,-1]==1)[0]
	datas = datas[idx]
	y, x, k = datas.T
	k, galn = WLanalysis.coords2grid(x, y, array([k,]))
	for sigmaG in sigmaG_arr:
		print 'createBadFieldMask sf, sigmaG:', sf, sigmaG
		Allmask = WLanalysis.readFits(mask_fcn(sigmaG, sf))#mask for all field
		badmask_fn = badmask_fcn(sigmaG, sf)#file name for bad pointing mask, which is 75% area of Allmask
		galn_smooth = snd.filters.gaussian_filter(galn.astype(float),sigmaG*PPA512, mode='constant')
		#smooth the galn grid
		Mmask = ones(shape=galn.shape)#create mask grid
		Mmask[where(galn_smooth < ngal_cut)]=0#find the low density region in galn_smooth
		Mmask = adHocFix(Mmask,sf)
		Mmask *= Allmask#since I didn't do redshift cut in badmask, so here it takes care of it, since ALl mask has redshift cuts
		WLanalysis.writeFits(Mmask, badmask_fn)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import coords2grid [as 别名]
def fileGen(i):
	'''
	Input:
	i range from (1, 2..13)
	Return:
	Me1 = e1*w
	Me2 = (e2-c2)*w
	Mw = (1+m)*w
	galn = number of galaxies per pixel
	'''
	Me1_fn = KS_dir+'CFHT_subfield%02d_Me1.fits'%(i)
	Me2_fn = KS_dir+'CFHT_subfield%02d_Me2.fits'%(i)
	Mw_fn = KS_dir+'CFHT_subfield%02d_Mw.fits'%(i)
	galn_fn = KS_dir+'CFHT_subfield%02d_galn.fits'%(i)
	
	print 'fileGen', i
	if WLanalysis.TestComplete((Me1_fn,Me2_fn,Mw_fn,galn_fn),rm = True):
		Me1 = WLanalysis.readFits(Me1_fn)
		Me2 = WLanalysis.readFits(Me2_fn)
		Mw =  WLanalysis.readFits(Mw_fn)
		galn =WLanalysis.readFits(galn_fn)
	else:
		ifile = np.genfromtxt(full_dir+'full_subfield'+str(i) ,usecols=[0, 1, 2, 3, 4, 9, 10, 11, 16, 17])
		# cols: y, x, z_peak, z_rnd1, z_rnd2, e1, e2, w, m, c2

		#redshift cut 0.2< z <1.3
		zs = ifile[:,[2,3,4]]
		print 'zs'
		idx = np.where((amax(zs,axis=1) <= zmax) & (amin(zs,axis=1) >= zmin))[0]
		
		y, x, z_peak, z_rnd1, z_rnd2, e1, e2, w, m, c2 = ifile[idx].T

		k = array([e1*w, (e2-c2)*w, (1+m)*w])
		Ms, galn = WLanalysis.coords2grid(x, y, k)
		print 'coords2grid'
		Me1, Me2, Mw = Ms
		WLanalysis.writeFits(Me1,Me1_fn)
		WLanalysis.writeFits(Me2,Me2_fn)
		WLanalysis.writeFits(Mw,Mw_fn)
		WLanalysis.writeFits(galn,galn_fn)
	return Me1, Me2, Mw, galn

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import coords2grid [as 别名]
	for i in arange(1,14):
		sG_arr = array([1, 3.5, 5.3, 8.9])
		for sG in sG_arr:
			KS_45_fn = test_dir+'KS_45_%i_%i.fit'%(i,sG)
			KS_rand_fn = test_dir+'KS_rand_%i_%i.fit'%(i,sG)
			print i, sG
			if os.path.isfile(KS_45_fn) and os.path.isfile(KS_rand_fn):
				KS_45=WLanalysis.readFits(KS_45_fn)
				KS_rand=WLanalysis.readFits(KS_rand_fn)
			
			else:
				print 'generating KS'
				y, x, e1, e2, w = WLanalysis.readFits(test_dir+'yxew_subfield%i_zcut0213.fit'%(i)).T
				e1_45, e2_45 = WLanalysis.rndrot(e1, e2, deg=45)
				e1_rand, e2_rand = WLanalysis.rndrot(e1, e2, iseed=0)
				mat_e1_45,mat_e2_45,mat_e1_rand,mat_e2_rand, Mw = WLanalysis.coords2grid(x, y, array([e1_45*w, e2_45*w, e1_rand*w, e2_rand*w, w]) )[0]

				mat_e1_45_smoothed  = WLanalysis.weighted_smooth(mat_e1_45  , Mw, sigmaG=sG)
				mat_e2_45_smoothed  = WLanalysis.weighted_smooth(mat_e2_45  , Mw, sigmaG=sG)
				mat_e1_rand_smoothed= WLanalysis.weighted_smooth(mat_e1_rand, Mw, sigmaG=sG)
				mat_e2_rand_smoothed= WLanalysis.weighted_smooth(mat_e2_rand, Mw, sigmaG=sG)
				KS_45 =WLanalysis.KSvw(mat_e1_45_smoothed,mat_e2_45_smoothed)
				KS_rand=WLanalysis.KSvw(mat_e1_rand_smoothed,mat_e2_rand_smoothed)
				
				
				WLanalysis.writeFits(KS_45,KS_45_fn)
				WLanalysis.writeFits(KS_rand,KS_rand_fn)
			
			kappa_45 = KS_45.flatten()
			std_45=std(kappa_45)
			#hist_45,binedges45 = histogram(kappa_45,bins=100,range=(-5*std_45,5*std_45))

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import coords2grid [as 别名]
		
		icat = cat_gen(Wx).T #
		ra, dec, redshift, weight, MAGi, Mhalo, Rvir, DL = icat
		##########next 3 lines, cut at z=0.4 for source galaxies
		k = np.load(obsPK_dir+'kappa_predict_W%i.npy'%(Wx))
		idx_k = where(k>0)[0]
		k = k[idx_k]
		icat = icat[:,idx_k]
		
		#####################################
		f_Wx = WLanalysis.gnom_fun(center)	
		y, x = array(f_Wx(icat[:2]))
		weight = icat[3]
		#k = np.load(obsPK_dir+'kappa_predict_W%i.npy'%(Wx))#kappa_predict_Mmax2e15_W%i.npy
		A, galn = WLanalysis.coords2grid(x, y, array([k*weight, weight, k]), size=isize)
		Mkw, Mw, Mk = A
		###########################################
		
		for sigmaG in  (0.5, 1.0, 1.8, 3.5, 5.3, 8.9):
			print Wx, sigmaG
			
			mask0 = maskGen(Wx, 0.5, sigmaG)
			mask = WLanalysis.smooth(mask0, 5.0)
			################ make maps ######################
			kmap_predict = WLanalysis.weighted_smooth(Mkw, Mw, PPA=PPA512, sigmaG=sigmaG)
			kmap_predict*=mask
			np.save(obsPK_dir+'maps/kmap_W%i_predict_sigmaG%02d.npy'%(Wx, sigmaG*10), kmap_predict)
			###########################################
if plot_predict_maps:
	def plot_predict_maps_fcn(WxsigmaG):

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import coords2grid [as 别名]
def SumSplitFile2Grid(Wx):
	'''For Wx field, read in each split file, 
	and create e1, e2 grid for mass construction.
	Input: Wx=1,2,3,4
	Output: (Me1, Me2, Mw, galn) split in each redshift bins'''
	isize = sizes[Wx-1]
	ishape = (len(zbins), isize, isize)
	ishape_hi = (len(zbins)-1, isize, isize)#no need to do hi for zcut=1.3 since it's everything
	Me1_hi = zeros(shape=ishape_hi)
	Me2_hi = zeros(shape=ishape_hi)#hi is for higher redshift bins, lo is lower redshift
	Mw_hi = zeros(shape=ishape_hi)
	#Mk_hi = zeros(shape=ishape_hi)
	galn_hi = zeros(shape=ishape_hi)

	Me1_lo = zeros(shape=ishape)
	Me2_lo = zeros(shape=ishape)
	Mw_lo = zeros(shape=ishape)
	#Mk_lo = zeros(shape=ishape)
	galn_lo = zeros(shape=ishape)
	
	Wfiles = os.listdir(W_dir(Wx))#get the list of split file for Wx
	for iW in Wfiles:
		datas = WLanalysis.readFits(W_dir(Wx)+iW)
		#cols: x, y, ra, dec, e1, e2, w, r, snr, m, c2, mag, z_peak, z_rand1, z_rand2
		z = datas.T[-3]#z_peak, -2 is z_rand1, -1 is z_rand2
		i = 0 #zbin count
		for zcut in zbins:
			idx0 = where(z<zcut)[0]
			idx1 = where(z>=zcut)[0]
			for idx in [idx0,idx1]:
				y, x, e1, e2, w, m = (datas[idx].T)[[0,1,4,5,6,9]]#note x, y is reversed in python
				k = array([e1*w, e2*w, (1+m)*w])
				x = radians(x)
				y = radians(y)
				print 'W'+str(Wx), iW, 'coords2grid, zbin =',zbins[i]
				A, galn = WLanalysis.coords2grid(x, y, k, size=isize)
				if len(idx)==0:#no need to calculate hi bin for zcut=1.3
					continue
				elif idx[0] == idx0[0]:
					Me1_lo[i] += A[0]
					Me2_lo[i] += A[1]
					Mw_lo[i] += A[2]
					galn_lo[i] += galn
				
				else:
					Me1_hi[i] += A[0]
					Me2_hi[i] += A[1]
					Mw_hi[i] += A[2]
					galn_hi[i] += galn
			i+=1
	print 'Done collecting small fields for W'+str(Wx)
	
	for i in range(len(zbins)):
		for hl in ('lo','hi'):
			Me1_fn = cat_dir+'Me_Mw_galn/W%i_Me1w_%s_%s.fit'%(Wx, zbins[i],hl)
			Me2_fn = cat_dir+'Me_Mw_galn/W%i_Me2w_%s_%s.fit'%(Wx, zbins[i],hl)
			Mw_fn = cat_dir+'Me_Mw_galn/W%i_Mwm_%s_%s.fit'%(Wx, zbins[i],hl)
			galn_fn = cat_dir+'Me_Mw_galn/W%i_galn_%s_%s.fit'%(Wx, zbins[i],hl)
			if hl=='hi' and i==len(zbins)-1:
				continue
			elif hl=='lo':
				WLanalysis.writeFits(Me1_lo[i],Me1_fn, rewrite = True)
				WLanalysis.writeFits(Me2_lo[i],Me2_fn, rewrite = True)
				WLanalysis.writeFits(Mw_lo[i],Mw_fn, rewrite = True)
				WLanalysis.writeFits(galn_lo[i],galn_fn, rewrite = True)
			else:
				WLanalysis.writeFits(Me1_hi[i],Me1_fn, rewrite = True)
				WLanalysis.writeFits(Me2_hi[i],Me2_fn, rewrite = True)
				WLanalysis.writeFits(Mw_hi[i],Mw_fn, rewrite = True)
				WLanalysis.writeFits(galn_hi[i],galn_fn, rewrite = True)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import coords2grid [as 别名]
def galn_gen(i):
	print i
	y, x, z = WLanalysis.readFits(emucat_dir+'emulator_subfield%i_zcut0213.fit'%(i)).T
	Mz, galn = WLanalysis.coords2grid(x, y, np.array([z,]))
	WLanalysis.writeFits(galn, galn_fn)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import coords2grid [as 别名]
	z = WLanalysis.readFits('/Users/jia/weaklensing/CFHTLenS/catalogue/emulator_galpos_zcut0213/emulator_subfield1_zcut0213.fit').T[-1]
	
	znoise1 = eobs_analytical (z)
	znoise2 = eobs_analytical (z)
	znoise27 = np.random.normal(0, 0.29,size=len(z))
	
	s1 = s1o*(1+m)
	s2 = s2o*(1+m)
	
	eint1, eint2 = WLanalysis.rndrot(e1, e2, iseed=1000)
	eint1_45, eint2_45 = WLanalysis.rndrot(e1, e2, iseed=1000,deg=45.0)

	e1_reduce, e2_reduce = eobs_fun(s1, s2, kappa, eint1, eint2)
	e1_add, e2_add = s1+eint1, s2+eint2
	
	A, galn = WLanalysis.coords2grid(x, y, array([e1_reduce*w, e2_reduce*w, w*(1+m), s1o, s2o, kappa, e1_add*w, e2_add*w, s1+znoise1, s2+znoise2, kappa+znoise1]))
	#Mk, Ms1, Ms2 = A
	Me1, Me2, Mw, Ms1, Ms2, Mk, Me1add, Me2add, Ms1n, Ms2n, Mkn = A
	
	B, galn = WLanalysis.coords2grid(x, y, array([znoise1, eint1, eint2, eint1_45, eint2_45, znoise27]))
	Mn, Men1, Men2, Men1b, Men2b, Mn27 = B
	### pure analytical noise ######
	Mn_smooth =  WLanalysis.weighted_smooth(Mn, galn, sigmaG=0.5)
	Mn27_smooth =  WLanalysis.weighted_smooth(Mn27, galn, sigmaG=0.5)
	Men1_smooth =  WLanalysis.weighted_smooth(Men1, galn, sigmaG=0.5)
	Men2_smooth =  WLanalysis.weighted_smooth(Men2, galn, sigmaG=0.5)
	Men1b_smooth =  WLanalysis.weighted_smooth(Men1b, galn, sigmaG=0.5)
	Men2b_smooth =  WLanalysis.weighted_smooth(Men2b, galn, sigmaG=0.5)
	Nmap = WLanalysis.KSvw(Men1_smooth, Men2_smooth)
	Nmap45 = WLanalysis.KSvw(Men1b_smooth, Men2b_smooth)