add_sky_background.py 6.56 KB
Newer Older
1
2
3
import numpy as np
import galsim
from ObservationSim.Straylight import calculateSkyMap_split_g
Zhang Xin's avatar
Zhang Xin committed
4
from ObservationSim.Instrument import FilterParam
5

6
7
8
from astropy.time import Time
from datetime import datetime, timezone

9
def add_sky_background_sci(self, chip, filt, tel, pointing, catalog, obs_param):
10
11
12
13
14
15
16
    
    # Get exposure time
    if (obs_param) and ("exptime" in obs_param) and (obs_param["exptime"] is not None):
        exptime = obs_param["exptime"]
    else:
        exptime = pointing.exp_time
    
17
18
19
20
21
22
    flat_normal = np.ones_like(chip.img.array)
    if obs_param["flat_fielding"] == True:
        flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
    if obs_param["shutter_effect"] == True:
        flat_normal = flat_normal * chip.shutter_img
        flat_normal = np.array(flat_normal, dtype='float32')
23
24
        self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT'], values = [True])
    else:
25
        self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN1','SHTCLOS0'], values = [True,self.h_ext['SHTCLOS1'],self.h_ext['SHTOPEN0']])
26
27
28
29
30
31
32
33
34
35
36
37
    
    if obs_param["enable_straylight_model"]:
        # Filter.sky_background, Filter.zodical_spec will be updated
        filt.setFilterStrayLightPixel(
            jtime = pointing.jdt,
            sat_pos = np.array([pointing.sat_x, pointing.sat_y, pointing.sat_z]),
            pointing_radec = np.array([pointing.ra,pointing.dec]),
            sun_pos = np.array([pointing.sun_x, pointing.sun_y, pointing.sun_z]))
        self.chip_output.Log_info("================================================")
        self.chip_output.Log_info("sky background + stray light pixel flux value: %.5f"%(filt.sky_background))
    
    if chip.survey_type == "photometric":
38
        sky_map = filt.getSkyNoise(exptime = exptime)
39
40
41
42
43
44
45
46
47
48
49
50
51
        sky_map = sky_map * np.ones_like(chip.img.array) * flat_normal
        sky_map = galsim.Image(array=sky_map)
    else:
        # chip.loadSLSFLATCUBE(flat_fn='flat_cube.fits')
        sky_map = calculateSkyMap_split_g(
                skyMap=flat_normal,
                blueLimit=filt.blue_limit,
                redLimit=filt.red_limit,
                conf=chip.sls_conf,
                pixelSize=chip.pix_scale,
                isAlongY=0,
                flat_cube=chip.flat_cube, 
                zoldial_spec = filt.zodical_spec)
52
        sky_map = (sky_map + filt.sky_background)*exptime
53
54
55
56
57
58
59
60
61
62
63
    
    # sky_map = sky_map * tel.pupil_area * obs_param["exptime"]
    chip.img += sky_map
    return chip, filt, tel, pointing

def add_sky_flat_calibration(self, chip, filt, tel, pointing, catalog, obs_param):

    if not hasattr(self, 'h_ext'):
        _, _ = self.prepare_headers(chip=chip, pointing=pointing)
    chip_wcs = galsim.FitsWCS(header = self.h_ext)

64
65
66
67
68
69
    # Get exposure time
    if (obs_param) and ("exptime" in obs_param) and (obs_param["exptime"] is not None):
        exptime = obs_param["exptime"]
    else:
        exptime = pointing.exp_time
    
70
71
72
73
    skyback_level = obs_param["flat_level"]

    filter_param = FilterParam()
    sky_level_filt = obs_param["flat_level_filt"]
74
    norm_scaler = skyback_level/exptime / filter_param.param[sky_level_filt][5]
75
76
77
78
79
80
81

    flat_normal = np.ones_like(chip.img.array)
    if obs_param["flat_fielding"] == True:
        flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
    if obs_param["shutter_effect"] == True:
        flat_normal = flat_normal * chip.shutter_img
        flat_normal = np.array(flat_normal, dtype='float32')
82
83
        if self.overall_config["output_setting"]["shutter_output"] == True:    # output 16-bit shutter effect image with pixel value <=65535
            shutt_gsimg = galsim.ImageUS(chip.shutter_img*6E4)
84
            shutt_gsimg.write("%s/ShutterEffect_%s_1.fits" % (self.chip_output.subdir, str(chip.chipID).rjust(2, '0')))
85
86
        self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT'], values = [True])
    else:
87
        self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN1','SHTCLOS0'], values = [True,self.h_ext['SHTCLOS1'],self.h_ext['SHTOPEN0']])
88
89
90
    

    if chip.survey_type == "photometric":
91
        sky_map = flat_normal * np.ones_like(chip.img.array) * norm_scaler * filter_param.param[chip.filter_type][5] / tel.pupil_area * exptime
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
    elif chip.survey_type == "spectroscopic":
        # flat_normal = np.ones_like(chip.img.array)
        if obs_param["flat_fielding"] == True:
            
            flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
        if obs_param["shutter_effect"] == True:
            
            flat_normal = flat_normal * chip.shutter_img
            flat_normal = np.array(flat_normal, dtype='float32')
        sky_map = calculateSkyMap_split_g(
            skyMap=flat_normal,
            blueLimit=filt.blue_limit,
            redLimit=filt.red_limit,
            conf=chip.sls_conf,
            pixelSize=chip.pix_scale,
            isAlongY=0,
            flat_cube=chip.flat_cube)
109
        sky_map = sky_map * norm_scaler * exptime
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
    
    chip.img += sky_map
    return chip, filt, tel, pointing

def add_sky_background(self, chip, filt, tel, pointing, catalog, obs_param):
    if not hasattr(self, 'h_ext'):
        _, _ = self.prepare_headers(chip=chip, pointing=pointing)
    chip_wcs = galsim.FitsWCS(header = self.h_ext)

    if "flat_level" not in obs_param or "flat_level_filt" not in obs_param:
        chip, filt, tel, pointing = self.add_sky_background_sci(chip, filt, tel, pointing, catalog, obs_param)
    else:
        if obs_param.get('flat_level') is None or obs_param.get('flat_level_filt')is None:
            chip, filt, tel, pointing = self.add_sky_background_sci(chip, filt, tel, pointing, catalog, obs_param)
        else:
            chip, filt, tel, pointing = self.add_sky_flat_calibration(chip, filt, tel, pointing, catalog, obs_param)
126
127
128
129
130
131
132
133
134
135
136
137
138

    # renew header info
    datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
    datetime_obs = datetime_obs.replace(tzinfo=timezone.utc)
    t_obs = Time(datetime_obs)
    
    ##ccd刷新2s,等待0.5s,开始曝光
    t_obs_renew = Time(t_obs.mjd - (2.+0.5) / 86400., format="mjd")

    t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
    self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])

    #dark time : 曝光时间+刷新后等带时间0.5s+关闭快门时间1.5s+管快门后读出前等待0.5s
Zhang Xin's avatar
Zhang Xin committed
139
    self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.+0.0+0.0+pointing.exp_time])
140
141


142
    return chip, filt, tel, pointing