SkybackgroundMap.py 11.3 KB
Newer Older
1
2
from ObservationSim.MockObject.SpecDisperser import SpecDisperser
from ObservationSim.MockObject.SpecDisperser import rotate90
Fang Yuedong's avatar
Fang Yuedong committed
3
4
5
6
7
8
9

import galsim
import numpy as np
from astropy.table import Table
from scipy import interpolate

import galsim
10
import astropy.constants as cons
Fang Yuedong's avatar
Fang Yuedong committed
11
12
13

import os

Zhang Xin's avatar
Zhang Xin committed
14
15
import time

Fang Yuedong's avatar
Fang Yuedong committed
16
17
18
19
20
21
try:
    import importlib.resources as pkg_resources
except ImportError:
    # Try backported to PY<37 'importlib_resources'
    import importlib_resources as pkg_resources

Fang Yuedong's avatar
Fang Yuedong committed
22
23
24

###calculate sky map by sky SED

Fang Yuedong's avatar
Fang Yuedong committed
25
def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='sky_emiss_hubble_50_50_A.dat', conf=[''], pixelSize=0.074, isAlongY=0,
26
                            split_pos=3685, flat_cube = None, zoldial_spec = None):
27
28
29
30
    # skyMap = np.ones([yLen, xLen], dtype='float32')
    #
    # if isAlongY == 1:
    #     skyMap = np.ones([xLen, yLen], dtype='float32')
Fang Yuedong's avatar
Fang Yuedong committed
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

    # for i in range(len(conf)):
    #     conf[i] = os.path.join(SLSSIM_PATH, conf[i])
    conf1 = conf[0]
    conf2 = conf[0]
    if np.size(conf) == 2:
        conf2 = conf[1]

    skyImg = galsim.Image(skyMap, xmin=0, ymin=0)

    tbstart = blueLimit
    tbend = redLimit

    fimg = np.zeros_like(skyMap)

    fImg = galsim.Image(fimg)
47
    try:
48
        with pkg_resources.files('ObservationSim.Straylight.data.sky').joinpath(skyfn) as data_path:
49
50
            skySpec = np.loadtxt(data_path)
    except AttributeError:
51
        with pkg_resources.path('ObservationSim.Straylight.data.sky', skyfn) as data_path:
52
            skySpec = np.loadtxt(data_path)
Fang Yuedong's avatar
Fang Yuedong committed
53
    # skySpec = np.loadtxt(skyfn)
Fang Yuedong's avatar
Fang Yuedong committed
54
55
    spec = Table(np.array([skySpec[:, 0], skySpec[:, 1]]).T, names=('WAVELENGTH', 'FLUX'))

56
57
58
59
60
61
62
63
64
65
    if zoldial_spec is not None:
        deltL = 0.5
        lamb = np.arange(2000, 11000, deltL)

        speci = interpolate.interp1d(zoldial_spec['WAVELENGTH'], zoldial_spec['FLUX'])

        y = speci(lamb)
        # erg/s/cm2/A --> photo/s/m2/A
        s_flux = y * lamb / (cons.h.value * cons.c.value) * 1e-13
        spec = Table(np.array([lamb, s_flux]).T, names=('WAVELENGTH', 'FLUX'))
Fang Yuedong's avatar
Fang Yuedong committed
66
67
68
69
70
71
72
73
74
75
76
77
78
79
    if isAlongY == 0:
        directParm = 0
    if isAlongY ==1:
        directParm = 1

    if split_pos >= skyImg.array.shape[directParm]:
        skyImg1 = galsim.Image(skyImg.array)
        origin1 = [0, 0]
        # sdp = specDisperser.specDisperser(orig_img=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y,
        #                                   full_img=fimg, tar_spec=spec, band_start=tbstart, band_end=tbend,
        #                                   origin=origin1,
        #                                   conf=conf1)
        # sdp.compute_spec_orders()

80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
        y_len = skyMap.shape[0]
        x_len = skyMap.shape[1]
        delt_x = 100
        delt_y = 100

        sub_y_start_arr = np.arange(0, y_len, delt_y)
        sub_y_end_arr = sub_y_start_arr + delt_y
        sub_y_end_arr[-1] = min(sub_y_end_arr[-1], y_len)

        sub_x_start_arr = np.arange(0, x_len, delt_x)
        sub_x_end_arr = sub_x_start_arr + delt_x
        sub_x_end_arr[-1] = min(sub_x_end_arr[-1], x_len)

        for i,k1 in enumerate(sub_y_start_arr):
            sub_y_s = k1
            sub_y_e = sub_y_end_arr[i]

            sub_y_center = (sub_y_s+sub_y_e)/2.

            for j,k2 in enumerate(sub_x_start_arr):
                sub_x_s = k2
                sub_x_e = sub_x_end_arr[j]

                skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e])
                origin_sub = [sub_y_s, sub_x_s]
                sub_x_center = (sub_x_s + sub_x_e) / 2.

                sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center, origin=origin_sub,
                                tar_spec=spec,
                                band_start=tbstart, band_end=tbend,
                                conf=conf2,
xin's avatar
xin committed
111
                                flat_cube=flat_cube, ignoreBeam=['D','E'])
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

                spec_orders = sdp.compute_spec_orders()

                for k, v in spec_orders.items():
                    img_s = v[0]
                    origin_order_x = v[1]
                    origin_order_y = v[2]
                    ssImg = galsim.ImageF(img_s)
                    ssImg.setOrigin(origin_order_x, origin_order_y)
                    bounds = ssImg.bounds & fImg.bounds
                    if bounds.area() == 0:
                        continue
                    fImg[bounds] = fImg[bounds] + ssImg[bounds]

        # sdp = SpecDisperser(orig_img=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y, origin=origin1,
        #                 tar_spec=spec,
        #                 band_start=tbstart, band_end=tbend,
        #                 conf=conf2,
        #                 flat_cube=flat_cube, ignoreBeam=['D','E'])
        #
        # spec_orders = sdp.compute_spec_orders()
        #
        # for k, v in spec_orders.items():
        #     img_s = v[0]
        #     origin_order_x = v[1]
        #     origin_order_y = v[2]
        #     ssImg = galsim.ImageF(img_s)
        #     ssImg.setOrigin(origin_order_x, origin_order_y)
        #     bounds = ssImg.bounds & fImg.bounds
        #     if bounds.area() == 0:
        #         continue
        #     fImg[bounds] = fImg[bounds] + ssImg[bounds]
Fang Yuedong's avatar
Fang Yuedong committed
144
145
146
147
148
        

        
    else:

149
150
151
152
        # skyImg1 = galsim.Image(skyImg.array[:, 0:split_pos])
        # origin1 = [0, 0]
        # skyImg2 = galsim.Image(skyImg.array[:, split_pos:])
        # origin2 = [0, split_pos]
Fang Yuedong's avatar
Fang Yuedong committed
153
154
155
156
157
158
159

        # sdp = specDisperser.specDisperser(orig_img=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y,
        #                                   full_img=fimg, tar_spec=spec, band_start=tbstart, band_end=tbend,
        #                                   origin=origin1,
        #                                   conf=conf1)

        # sdp.compute_spec_orders()
160
161
162
        y_len = skyMap.shape[0]
        x_len = skyMap.shape[1]
        delt_x = 500
Zhang Xin's avatar
Zhang Xin committed
163
        delt_y = y_len
164
165
166
167

        sub_y_start_arr = np.arange(0, y_len, delt_y)
        sub_y_end_arr = sub_y_start_arr + delt_y
        sub_y_end_arr[-1] = min(sub_y_end_arr[-1], y_len)
Zhang Xin's avatar
Zhang Xin committed
168
169
        
        delt_x = split_pos-0
170
171
172
173
174
175
176
177
178
179
180
181
182
        sub_x_start_arr = np.arange(0, split_pos, delt_x)
        sub_x_end_arr = sub_x_start_arr + delt_x
        sub_x_end_arr[-1] = min(sub_x_end_arr[-1], split_pos)

        for i,k1 in enumerate(sub_y_start_arr):
            sub_y_s = k1
            sub_y_e = sub_y_end_arr[i]

            sub_y_center = (sub_y_s+sub_y_e)/2.

            for j,k2 in enumerate(sub_x_start_arr):
                sub_x_s = k2
                sub_x_e = sub_x_end_arr[j]
xin's avatar
xin committed
183
                # print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e)
Zhang Xin's avatar
Zhang Xin committed
184
                T1 = time.time()
185
186
187
188
189
190
191
192
                skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e])
                origin_sub = [sub_y_s, sub_x_s]
                sub_x_center = (sub_x_s + sub_x_e) / 2.

                sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center, origin=origin_sub,
                                tar_spec=spec,
                                band_start=tbstart, band_end=tbend,
                                conf=conf1,
xin's avatar
xin committed
193
                                flat_cube=flat_cube)
194
195
196
197
198
199
200
201
202
203
204
205
206

                spec_orders = sdp.compute_spec_orders()

                for k, v in spec_orders.items():
                    img_s = v[0]
                    origin_order_x = v[1]
                    origin_order_y = v[2]
                    ssImg = galsim.ImageF(img_s)
                    ssImg.setOrigin(origin_order_x, origin_order_y)
                    bounds = ssImg.bounds & fImg.bounds
                    if bounds.area() == 0:
                        continue
                    fImg[bounds] = fImg[bounds] + ssImg[bounds]
Zhang Xin's avatar
Zhang Xin committed
207
208
209
210
        
                T2 = time.time()

                print('time: %s ms'% ((T2 - T1)*1000))
211

Zhang Xin's avatar
Zhang Xin committed
212
        delt_x = x_len-split_pos
213
214
215
216
217
218
219
220
221
222
223
224
225
        sub_x_start_arr = np.arange(split_pos, x_len, delt_x)
        sub_x_end_arr = sub_x_start_arr + delt_x
        sub_x_end_arr[-1] = min(sub_x_end_arr[-1], x_len)

        for i, k1 in enumerate(sub_y_start_arr):
            sub_y_s = k1
            sub_y_e = sub_y_end_arr[i]

            sub_y_center = (sub_y_s + sub_y_e) / 2.

            for j, k2 in enumerate(sub_x_start_arr):
                sub_x_s = k2
                sub_x_e = sub_x_end_arr[j]
xin's avatar
xin committed
226
                # print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e)
Zhang Xin's avatar
Zhang Xin committed
227
228
                
                T1 = time.time()
229
230
231
232
233
234
235
236
237

                skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e])
                origin_sub = [sub_y_s, sub_x_s]
                sub_x_center = (sub_x_s + sub_x_e) / 2.

                sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center, origin=origin_sub,
                                tar_spec=spec,
                                band_start=tbstart, band_end=tbend,
                                conf=conf2,
xin's avatar
xin committed
238
                                flat_cube=flat_cube)
239
240
241
242
243
244
245
246
247
248
249
250
251

                spec_orders = sdp.compute_spec_orders()

                for k, v in spec_orders.items():
                    img_s = v[0]
                    origin_order_x = v[1]
                    origin_order_y = v[2]
                    ssImg = galsim.ImageF(img_s)
                    ssImg.setOrigin(origin_order_x, origin_order_y)
                    bounds = ssImg.bounds & fImg.bounds
                    if bounds.area() == 0:
                        continue
                    fImg[bounds] = fImg[bounds] + ssImg[bounds]
Zhang Xin's avatar
Zhang Xin committed
252
253
254
                T2 = time.time()

                print('time: %s ms'% ((T2 - T1)*1000))
Fang Yuedong's avatar
Fang Yuedong committed
255
256
257
258
259
260
261
262
263
264
265

    if isAlongY == 1:
        fimg, tmx, tmy = rotate90(array_orig=fImg.array, xc=0, yc=0, isClockwise=0)
    else:
        fimg = fImg.array

    fimg = fimg * pixelSize * pixelSize

    return fimg
    
def calculateSkyMap(xLen=9232, yLen=9126, blueLimit=4200, redLimit=6500,
Fang Yuedong's avatar
Fang Yuedong committed
266
                    skyfn='sky_emiss_hubble_50_50_A.dat', conf='', pixelSize=0.074, isAlongY=0):
Fang Yuedong's avatar
Fang Yuedong committed
267
268
269
270
271
272
273
274
275
276
277
278
    skyMap = np.ones([yLen, xLen], dtype='float32')

    if isAlongY == 1:
        skyMap = np.ones([xLen, yLen], dtype='float32')

    skyImg = galsim.Image(skyMap)

    tbstart = blueLimit
    tbend = redLimit

    fimg = np.zeros_like(skyMap)
    fImg = galsim.Image(fimg)
279
    try:
280
        with pkg_resources.files('ObservationSim.Straylight.data.sky').joinpath(skyfn) as data_path:
281
282
            skySpec = np.loadtxt(data_path)
    except AttributeError:
283
        with pkg_resources.path('ObservationSim.Straylight.data.sky', skyfn) as data_path:
284
            skySpec = np.loadtxt(data_path)
Fang Yuedong's avatar
Fang Yuedong committed
285
286
    # skySpec = np.loadtxt(skyfn)
    
Fang Yuedong's avatar
Fang Yuedong committed
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
    spec = Table(np.array([skySpec[:, 0], skySpec[:, 1]]).T, names=('WAVELENGTH', 'FLUX'))

    sdp = SpecDisperser(orig_img=skyImg, xcenter=skyImg.center.x, ycenter=skyImg.center.y, origin=[1, 1],
                        tar_spec=spec,
                        band_start=tbstart, band_end=tbend,
                        conf=conf)

    spec_orders = sdp.compute_spec_orders()

    for k, v in spec_orders.items():
        img_s = v[0]
        origin_order_x = v[1]
        origin_order_y = v[2]
        ssImg = galsim.ImageF(img_s)
        ssImg.setOrigin(origin_order_x, origin_order_y)
        bounds = ssImg.bounds & fImg.bounds
        fImg[bounds] = fImg[bounds] + ssImg[bounds]

    if isAlongY == 1:
        fimg, tmx, tmy = rotate90(array_orig=fImg.array, xc=0, yc=0, isClockwise=0)
    else:
        fimg = fImg.array
        
    fimg = fimg * pixelSize * pixelSize

    return fimg