from ObservationSim.MockObject.SpecDisperser import SpecDisperser from ObservationSim.MockObject.SpecDisperser import rotate90 import galsim import numpy as np from astropy.table import Table from scipy import interpolate import galsim import os import time try: import importlib.resources as pkg_resources except ImportError: # Try backported to PY<37 'importlib_resources' import importlib_resources as pkg_resources ###calculate sky map by sky SED def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='sky_emiss_hubble_50_50_A.dat', conf=[''], pixelSize=0.074, isAlongY=0, split_pos=3685, flat_cube = None): # skyMap = np.ones([yLen, xLen], dtype='float32') # # if isAlongY == 1: # skyMap = np.ones([xLen, yLen], dtype='float32') # 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) try: with pkg_resources.files('ObservationSim.MockObject.data').joinpath(skyfn) as data_path: skySpec = np.loadtxt(data_path) except AttributeError: with pkg_resources.path('ObservationSim.MockObject.data', skyfn) as data_path: skySpec = np.loadtxt(data_path) # skySpec = np.loadtxt(skyfn) spec = Table(np.array([skySpec[:, 0], skySpec[:, 1]]).T, names=('WAVELENGTH', 'FLUX')) 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() 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, 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] # 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] else: # skyImg1 = galsim.Image(skyImg.array[:, 0:split_pos]) # origin1 = [0, 0] # skyImg2 = galsim.Image(skyImg.array[:, split_pos:]) # origin2 = [0, split_pos] # 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() y_len = skyMap.shape[0] x_len = skyMap.shape[1] delt_x = 500 delt_y = y_len 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) delt_x = split_pos-0 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] # print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e) T1 = time.time() 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, flat_cube=flat_cube) 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] T2 = time.time() print('time: %s ms'% ((T2 - T1)*1000)) delt_x = x_len-split_pos 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] # print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e) T1 = time.time() 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, flat_cube=flat_cube) 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] T2 = time.time() print('time: %s ms'% ((T2 - T1)*1000)) 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, skyfn='sky_emiss_hubble_50_50_A.dat', conf='', pixelSize=0.074, isAlongY=0): 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) try: with pkg_resources.files('ObservationSim.MockObject.data').joinpath(skyfn) as data_path: skySpec = np.loadtxt(data_path) except AttributeError: with pkg_resources.path('ObservationSim.MockObject.data', skyfn) as data_path: skySpec = np.loadtxt(data_path) # skySpec = np.loadtxt(skyfn) 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