Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
# NAME:
# indexFits_hdf5
# PURPOSE:
# Return a healpix indexed catalog from a set of Fits
# CALLING:
# write_StampsIndex(dir_cat=dir_temp)
# INPUTS:
# dir_cat - the directory of Fits catalog, "<dir_cat>/stampCats/*.fits"
# OUTPUTS:
# <dir_cat>/stampCatsIndex.hdf5
# OPTIONAL:
# test_fits(nfits=100, dir_cat=None) - generate a set of Fits by galsim gaussian
# HISTORY:
# Written by Chengliang Wei, 13 Apr. 2023
# Included by csst-simulation, C.W. 25 Apr. 2023
#
#
import os
import numpy as np
import astropy.io.fits as fitsio
import h5py
import healpy
import galsim
def test_fits(nfits=100, dir_cat=None):
for ifits in range(nfits):
gal = galsim.Gaussian(sigma=np.random.uniform(0.2, 0.3)).shear(g1=np.random.uniform(-0.5, 0.5), g2=np.random.uniform(-0.5, 0.5))
arr = gal.drawImage(nx=64, ny=64, scale=0.074).array
hdu = fitsio.PrimaryHDU()
hdu.data = arr
hdu.header.set('index', ifits)
hdu.header.set('ra', 60.+np.random.uniform(-0.2, 0.2))
hdu.header.set('dec', -40.+np.random.uniform(-0.2, 0.2))
hdu.header.set('mag_g', 22+np.random.uniform(-1, 1))
hdu.header.set('pixScale', 0.074)
fout=dir_cat+"stampCats/testStamp_{:}.fits".format(ifits)
if os.path.exists(fout):
os.remove(fout)
hdu.writeto(fout)
def write_StampsIndex(dir_cat=None, DEBUG=False):
MAXNUMBERINDEX = 10000
NSIDE = 128
fp = h5py.File(dir_cat+'stampCatsIndex.hdf5', 'w')
grp1 = fp.create_group('Stamps')
dataSet_Size = np.zeros(healpy.nside2npix(NSIDE), dtype=np.int64)
fitsList = os.listdir(dir_cat+'stampCats/') #获取fits文件列表
for istamp in range(len(fitsList)):
print(istamp, ': ', fitsList[istamp], end='\r')
hdu=fitsio.open(dir_cat+"stampCats/"+fitsList[istamp])
tra = hdu[0].header['RA']
tdec= hdu[0].header['DEC']
healpixID= healpy.ang2pix(NSIDE, tra, tdec, nest=False, lonlat=True)
if not(str(healpixID) in grp1):
grp2 = grp1.create_group(str(healpixID))
else:
grp2 = grp1[str(healpixID)]
if not('ra' in grp2):
dset_ra = grp2.create_dataset('ra', (0,), dtype='f16' , maxshape=(MAXNUMBERINDEX, ))
dset_dec= grp2.create_dataset('dec', (0,), dtype='f16', maxshape=(MAXNUMBERINDEX, ))
dt = h5py.special_dtype(vlen=str)
dset_fn = grp2.create_dataset('filename', (0,), dtype=dt, maxshape=(MAXNUMBERINDEX, ))
else:
dset_ra = grp2['ra']
dset_dec = grp2['dec']
dset_fn = grp2['filename']
dataSet_Size[healpixID] = dataSet_Size[healpixID]+1
grp2['ra'].resize((dataSet_Size[healpixID],))
grp2['dec'].resize((dataSet_Size[healpixID],))
grp2['filename'].resize((dataSet_Size[healpixID],))
dset_ra[dataSet_Size[healpixID]-1] = tra
dset_dec[dataSet_Size[healpixID]-1]= tdec
dset_fn[dataSet_Size[healpixID]-1]= fitsList[istamp]
fp.close()
if DEBUG:
print('\n')
ff = h5py.File(dir_cat+"stampCatsIndex.hdf5","r")
ss = 0
for kk in ff['Stamps'].keys():
print(kk, ff['Stamps'][kk]['ra'].size)
ss = ss+ff['Stamps'][kk]['ra'].size
print(ss)
if __name__ == '__main__':
dir_temp = "./Catalog_test/"
#test_fits(dir_cat=dir_temp)
write_StampsIndex(dir_cat=dir_temp)