forked from spacetelescope/ramp_simulator
-
Notifications
You must be signed in to change notification settings - Fork 0
/
ramp_simulator.py
executable file
·6349 lines (5176 loc) · 310 KB
/
ramp_simulator.py
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
104
105
106
107
108
109
110
111
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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
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
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#! /usr/bin/env python
'''
Rewritten version of Kevin Volk's rampsim.py. Updated to be more general
and flexible such that other instruments can use the code
'''
__version__ = "1.1" # Fairly major update to remove dependence on
# JWST pipeline if requested. Also reworked
# the linearized dark ramp input to work.
import argparse,sys, glob, os
import yaml
from copy import copy
import scipy.ndimage.interpolation as interpolation
import scipy.signal as s1
import numpy as np
import math,random
from astropy.io import fits,ascii
from astropy.table import Table
from astropy.time import Time, TimeDelta
from itertools import izip
import datetime
from time import sleep
from astropy.modeling.models import Sersic2D
import rotations #this is Colin's set of rotation functions
from asdf import AsdfFile
import polynomial #more of Colin's functions
from astropy.convolution import convolve
import read_siaf_table
import moving_targets
import read_fits
import set_telescope_pointing_separated as set_telescope_pointing
#make sure you query the build 7.1 engineering database
os.environ['ENG_RESTFUL_URL'] = 'http://iwjwdmsbemwebv.stsci.edu/JWDMSEngFqAccB71/TlmMnemonicDataSrv.svc/'
inst_list = ['nircam']
modes = {'nircam':['imaging','moving_target','wfss','coron']}
inst_abbrev = {'nircam':'NRC'}
pixelScale = {'nircam':{'sw':0.031,'lw':0.063}}
full_array_size = {'nircam':2048}
allowedOutputFormats = ['DMS']
class RampSim():
def __init__(self):
#if a grism signal rate image is requested, expand
#the width and height of the signal rate image by this
#factor, so that the grism simulation software can
#track sources that are outside the requested subarray
#in order to calculate contamination.
self.grism_direct_factor = np.sqrt(2.)
#self.coord_adjust contains the factor by which the nominal output array size
#needs to be increased (used for TSO and WFSS modes), as well as the coordinate
#offset between the nominal output array coordinates, and those of the expanded
#array. These are needed mostly for TSO observations, where the nominal output
#array will not sit centered in the expanded output image.
self.coord_adjust = {'x':1.,'xoffset':0.,'y':1.,'yoffset':0.}
def run(self):
#NO LONGER SUPPORTED
print("WARNING WARNING!!")
print("This code is no longer supported!")
print("Updated, easier to use code is located at:")
print("https://github.com/spacetelescope/ramp_simulator")
print("We recommend you download the updated version")
print("of the simulator before running simulations.")
print("WARNING WARNING!!")
print('')
print('')
print('')
print('')
print('')
print('')
sleep(5)
#read in the parameter file
self.readParameterFile()
print('CHEAT CHEAT CHEAT!! manually setting output directory')
print('to the working directory so that I can finish WFSS')
print('example dataset!!!')
self.params['Output']['directory'] = './'
# expand locations in the parameter
# file to be full paths
self.fullPaths()
#if the user requests that a linearized dark be used
#rather than a raw dark, then set the use_JWST_pipeline
#option to false
#if self.params['Reffiles']['linearized_darkfile'] != 'None':
# self.params['Inst']['use_JWST_pipeline'] = False
#check the consistency of the input readpattern and nframe/nskip
self.readpatternCheck()
#from parameters, generate other useful variables
self.checkParams()
#read in the subarray definition file
self.readSubarrayDefinitionFile()
#define the subarray bounds from param file
self.getSubarrayBounds()
#read in the input dark current frame
if self.runStep['linearized_darkfile'] == False:
self.getBaseDark()
self.linDark = None
else:
print(self.params['Reffiles']['linearized_darkfile'])
self.readLinearDark()
self.dark = self.linDark
print("As read in:")
print(self.dark.data[0,0:3,400,400])
#print("self.dark as read in: {} {} {}".format(self.dark.data.shape,self.dark.zeroframe.shape,self.dark.sbAndRefpix.shape))
#print("self.linDark as read in: {} {} {}".format(self.linDark.data.shape,self.linDark.zeroframe.shape,self.linDark.sbAndRefpix.shape))
#make sure there is enough data (frames/groups) in the input integration to produce
#the proposed output integration
self.dataVolumeCheck(self.dark)
#print("After volume check:")
#print(self.dark.data[0,0:3,400,400])
#print("After volume check:")
#compare the requested number of integrations to the number
#of integrations in the input dark
print("Dark shape as read in: {}".format(self.dark.data.shape))
self.darkints()
print("Dark shape after copying integrations to match request: {}".format(self.dark.data.shape))
print("After darkints:")
print(self.dark.data[0,0:3,400,400])
# Put the input dark (or linearized dark) into the requested
# readout pattern
#self.dark,darkzeroframe,sbzeroframe = self.reorderDark(self.dark)
self.dark,sbzeroframe = self.reorderDark(self.dark)
print('DARK has been reordered to {} to match the input readpattern of {}'.format(self.dark.data.shape,self.dark.header['READPATT']))
#If a raw dark was read in, create linearized version here
#using the SSB pipeline. Better to do this
#on the full dark before cropping, so that reference pixels can be
#used in the processing.
if ((self.params['Inst']['use_JWST_pipeline']) & (self.runStep['linearized_darkfile'] == False)):
#Linearize the dark ramp via the SSB pipeline. Also save a diff image of
#the original dark minus the superbias and refpix subtracted dark, to use later.
#in order to linearize the dark, the JWST pipeline must be present, and self.dark
#will have to be translated back into a RampModel instance
#self.linDark, self.sbAndRefpixEffects = self.linearizeDark(self.dark)
self.linDark = self.linearizeDark(self.dark)
print("Linearized dark shape: {}".format(self.linDark.data.shape))
if self.params['Readout']['readpatt'].upper() == 'RAPID':
print("output is rapid, grabbing zero frame from linearized dark")
zeroModel = read_fits.Read_fits()
zeroModel.data = self.linDark.data[:,0,:,:] #can't grab all zeroframes (for multiple ints)'
zeroModel.sbAndRefpix = self.linDark.sbAndRefpix[:,0,:,:]
elif ((self.params['Readout']['readpatt'].upper() != 'RAPID') & (self.dark.zeroframe is not None)):
print("Now we need to linearize the zeroframe because the")
print("output readpattern is not RAPID")
#Now we need to linearize the zeroframe. Place it into a RampModel instance
#before running the pipeline steps
zeroModel = read_fits.Read_fits()
zeroModel.data = self.dark.zeroframe
zeroModel.header = self.linDark.header
zeroModel.header['NGROUPS'] = 1
zeroModel = self.linearizeDark(zeroModel)
else:
zeroModel = None
#Now crop self.linDark, self.dark, and zeroModel to requested subarray
self.dark = self.cropDark(self.dark)
self.linDark= self.cropDark(self.linDark)
if zeroModel is not None:
zeroModel = self.cropDark(zeroModel)
else:
#if no pipeline is run
zeroModel = read_fits.Read_fits()
zeroModel.data = self.dark.zeroframe
zeroModel.sbAndRefpix = sbzeroframe
#Crop the linearized dark to the requested
#subarray size
#THIS WILL CROP self.dark AS WELL SINCE
#self.linDark IS JUST A REFERENCE IN THE NON
#PIPELINE CASE!!
self.linDark = self.cropDark(self.linDark)
if zeroModel.data is not None:
zeroModel = self.cropDark(zeroModel)
print('cropped #2 the linearized dark current zero frame. {}'.format(zeroModel.data.shape))
#try:
# print('After Pipeline raw dark, data shape: {}'.format(self.dark.data.shape))
#except:
# print('After Pipeline raw dark, data extension not present!!')
#try:
# print('After Pipeline raw dark, zeroframe shape: {}'.format(self.dark.zeroframe.shape))
#except:
# print('After Pipeline raw dark, zeroframe extension not present!!')
#try:
# print('After Pipeline raw dark, sbAndRefpix shape: {}'.format(self.dark.sbAndRefpix.shape))
#except:
# print('After Pipeline raw dark, sbAndRefpix extension not present!!')
#try:
# print('After Pipeline linear dark, data shape: {}'.format(self.linDark.data.shape))
#except:
# print('After Pipeline linear dark, data extension not present!!')
#try:
# print('After Pipeline linear dark, zeroframe shape: {}'.format(self.linDark.zeroframe.shape))
#except:
# print('After Pipeline linear dark, zeroframe extension not present!!')
#try:
# print('After Pipeline linear dark, sbAndRefpix shape: {}'.format(self.linDark.sbAndRefpix.shape))
#except:
# print('After Pipeline linear dark, sbAndRefpix extension not present!!')
#try:
# print('After Pipeline zeromodel, data shape: {}'.format(zeroModel.data.shape))
#except:
# print('After Pipeline zeromodel, data extension not present!!')
#try:
# print('After Pipeline zeromodel, zeroframe shape: {}'.format(zeroModel.zeroframe.shape))
#except:
# print('After Pipeline zeromodel, zeroframe extension not present!!')
#try:
# print('After Pipeline zeromodel, sbAndRefpix shape: {}'.format(zeroModel.sbAndRefpix.shape))
#except:
# print('After Pipeline zeromodel, sbAndRefpix extension not present!!')
# print("self.dark after pipeline: {} {} {}".format(self.dark.data.shape,self.dark.zeroframe.shape,self.dark.sbAndRefpix.shape))
# try:
# print("self.linDark after pipline: {} {} {}".format(self.linDark.data.shape,self.linDark.zeroframe.shape,self.linDark.sbAndRefpix.shape))
# except:
# pass
#save the linearized dark for testing
if self.params['Output']['save_intermediates']:
h0=fits.PrimaryHDU()
h1 = fits.ImageHDU(self.linDark.data)
h2 = fits.ImageHDU(zeroModel.data)
hl=fits.HDUList([h0,h1,h2])
#hl.writeto(self.params['Output']['file'][0:-5] + '_linearizedDark.fits',overwrite=True)
hl.writeto(os.path.join(self.params['Output']['directory'],self.params['Output']['file'][0:-5] + '_linearizedDark.fits'),overwrite=True)
#save the cropped dark for testing
if self.params['Output']['save_intermediates']:
h0=fits.PrimaryHDU()
h1 = fits.ImageHDU(self.dark.data)
hl=fits.HDUList([h0,h1])
hl.writeto(os.path.join(self.params['Output']['directory'],self.params['Output']['file'][0:-5] + '_croppedDark.fits'),overwrite=True)
print("After cropping to subarrays, the shapes of self.dark, self.linDark are: {},{}"
.format(self.dark.data.shape,self.linDark.data.shape))
print("Size of the sbandrefpix extension: {}".format(self.linDark.sbAndRefpix.shape))
#Find the difference between the cropped original dark and the cropped linearized dark
#This will be used later to re-add the superbias and refpix-associated signal to the final
#output ramp.
#if self.params['newRamp']['combine_method'].upper() in ['HIGHSIG','PROPER']:
# if self.sbAndRefpixEffects is not None:
# self.sbAndRefpixEffects = self.sbAndRefpixEffects[:,:,self.subarray_bounds[1]:self.subarray_bounds[3]+1,self.subarray_bounds[0]:self.subarray_bounds[2]+1]
#calculate the exposure time of a single frame, based on the size of the subarray
self.calcFrameTime()
#calculate the rate of cosmic ray hits expected per frame
self.getCRrate()
#using the input instrument name, load appropriate
#instrument-specific dictionaries
self.instrument_specific_dicts(self.params['Inst']['instrument'].lower())
#If the output is a TSO ramp, use the slew rate and angle (and whether a grism
#direct image is requested) to determine how much larger than nominal the
#signal rate image should be
if self.params['Inst']['mode'] == 'moving_target' or self.params['Output']['grism_source_image']:
self.calcCoordAdjust()
#image dimensions
self.nominal_dims = np.array([self.subarray_bounds[3]-self.subarray_bounds[1]+1,self.subarray_bounds[2]-self.subarray_bounds[0]+1])
self.output_dims = (self.nominal_dims * np.array([self.coord_adjust['y'],self.coord_adjust['x']])).astype(np.int)
#generate the signal image from the input reference files
#The output contains signal in ADU that accumulate
#in one frametime.
frameimage = None
if self.params['Inst']['mode'].lower() != 'moving_target':
print('Creating signal rate image of synthetic inputs.')
frameimage = self.addedSignals()
#If moving target mode is used (i.e. tracking a non-sidereal target), then
#everything in the point source, galaxy, and extended source lists needs to
#be treated as a moving target.
if self.params['Inst']['mode'].lower() == 'moving_target':
print("need to adjust moving target work for multiple integrations! everything above has been modified")
#create the count rate image for the non-sidereal target(s)
nonsidereal_countrate,ra_vel,dec_vel,vel_flag = self.nonsidereal_CRImage(self.params['simSignals']['movingTargetToTrack'])
#expand into a RAPID ramp and convert from signal rate to signals
ns_yd,ns_xd = nonsidereal_countrate.shape
totframes = self.params['Readout']['ngroup'] * (self.params['Readout']['nframe']+self.params['Readout']['nskip'])
tmptimes = self.frametime * np.arange(1,totframes+1)
non_sidereal_ramp = np.zeros((totframes,ns_yd,ns_xd))
for i in xrange(totframes):
non_sidereal_ramp[i,:,:] = nonsidereal_countrate * tmptimes[i]
non_sidereal_zero = non_sidereal_ramp[0,:,:]
#Now we need to collect all the other sources (point sources, galaxies, extended)
#in the other input files, and treat them as targets which will move across
#the field of view during the exposure.
mtt_data_list = []
mtt_zero_list = []
if self.runStep['pointsource']:
#now ptsrc is a list, which we need to get into movingTargetInputs...
mtt_ptsrc = self.movingTargetInputs(self.params['simSignals']['pointsource'],'pointSource',MT_tracking=True,tracking_ra_vel=ra_vel,tracking_dec_vel=dec_vel,trackingPixVelFlag=vel_flag)
mtt_ptsrc_zero = mtt_ptsrc[0,:,:]
mtt_data_list.append(mtt_ptsrc)
mtt_zero_list.append(mtt_ptsrc_zero)
print("Done with creating moving targets from {}".format(self.params['simSignals']['pointsource']))
if self.runStep['galaxies']:
mtt_galaxies = self.movingTargetInputs(self.params['simSignals']['galaxyListFile'],'galaxies',MT_tracking=True,tracking_ra_vel=ra_vel,tracking_dec_vel=dec_vel,trackingPixVelFlag=vel_flag)
mtt_galaxies_zero = mtt_galaxies[0,:,:]
mtt_data_list.append(mtt_galaxies)
mtt_zero_list.append(mtt_galaxies_zero)
print("Done with creating moving targets from {}".format(self.params['simSignals']['galaxyListFile']))
#hh00 = fits.PrimaryHDU(mtt_galaxies)
#hhllist = fits.HDUList([hh00])
#hhllist.writeto('junk.fits',overwrite=True)
#print('Moving target ramp with galaxies only written to junk.fits')
if self.runStep['extendedsource']:
mtt_ext = self.movingTargetInputs(self.params['simSignals']['extended'],'extended',MT_tracking=True,tracking_ra_vel=ra_vel,tracking_dec_vel=dec_vel,trackingPixVelFlag=vel_flag)
mtt_ext_zero = mtt_ext[0,:,:]
mtt_data_list.append(mtt_ext)
mtt_zero_list.append(mtt_ext_zero)
print("Done with creating moving targets from {}".format(self.params['simSignals']['extended']))
print("Number of moving target ramps to combine {}".format(len(mtt_data_list)))
#add in the other objects which are not being tracked on
#(i.e. the sidereal targets)
if len(mtt_data_list) > 0:
for i in xrange(len(mtt_data_list)):
non_sidereal_ramp += mtt_data_list[i]
non_sidereal_zero += mtt_zero_list[i]
#if moving targets are requested (KBOs, asteroids, etc, NOT moving_target mode
#where the telescope slews), then create a RAPID integration which
#includes those targets
mov_targs_ramps = []
if self.runStep['movingTargets']:
print('Starting moving targets for point sources!')
mov_targs_ptsrc = self.movingTargetInputs(self.params['simSignals']['movingTargetList'],'pointSource',MT_tracking=False)
mov_targs_ramps.append(mov_targs_ptsrc)
#moving target using a sersic object
if self.runStep['movingTargetsSersic']:
print("Moving targets, sersic!")
mov_targs_sersic = self.movingTargetInputs(self.params['simSignals']['movingTargetSersic'],'galaxies',MT_tracking=False)
mov_targs_ramps.append(mov_targs_sersic)
#moving target using an extended object
if self.runStep['movingTargetsExtended']:
print("Extended moving targets!!!")
mov_targs_ext = self.movingTargetInputs(self.params['simSignals']['movingTargetExtended'],'extended',MT_tracking=False)
mov_targs_ramps.append(mov_targs_ext)
mov_targs_integration = None
if self.runStep['movingTargets'] or self.runStep['movingTargetsSersic'] or self.runStep['movingTargetsExtended']:
#Combine the ramps of the moving targets if there is more than one type
mov_targs_integration = mov_targs_ramps[0]
if len(mov_targs_ramps) > 1:
for i in xrange(1,len(mov_targs_ramps)):
mov_targs_integration += mov_targs_ramps[0]
#Combine the signal rate image and the moving target ramp. If signalramp is
#a ramp, then rearrange into the requested readout pattern
if mov_targs_integration is not None:
#if self.params['Inst']['mode'].lower() == 'imaging':
if self.params['Inst']['mode'].lower() != 'moving_target':
signalramp = self.combineSimulatedDataSources('countrate',frameimage,None,mov_targs_integration)
elif self.params['Inst']['mode'].lower() == 'moving_target':
signalramp = self.combineSimulatedDataSources('ramp',non_sidereal_ramp,non_sidereal_zero,mov_targs_integration)
#save the combined static+moving targets ramp
if self.params['Output']['save_intermediates']:
h0=fits.PrimaryHDU()
h1 = fits.ImageHDU(signalramp)
hl=fits.HDUList([h0,h1])
hl.writeto(os.path.join(self.params['Output']['directory'],self.params['Output']['file'][0:-5] + '_noiseless_static_plus_movingtargs_ramp.fits'),overwrite=True)
print('Ramp of static plus moving target simulated sources saved to {}'.format(self.params['Output']['file'][0:-5] + '_noiseless_static_plus_movingtargs_ramp.fits'))
else:
#if self.params['Inst']['mode'].lower() in ['imaging','wfss']:
if self.params['Inst']['mode'].lower() != 'moving_target':
num_frames = self.params['Readout']['ngroup'] * (self.params['Readout']['nframe'] + self.params['Readout']['nskip'])
yd,xd = frameimage.shape
signalramp = np.zeros((num_frames,yd,xd))
for i in xrange(num_frames):
signalramp[i,:,:] = frameimage * self.frametime * (i+1)
#signalzero = frameimage * self.frametime
elif self.params['Inst']['mode'].lower() == 'moving_target':
signalramp = non_sidereal_ramp
#signalzero = non_sidereal_zero
#save the combined static+moving targets ramp
if self.params['Output']['save_intermediates']:
h0=fits.PrimaryHDU()
h1 = fits.ImageHDU(signalramp)
hl=fits.HDUList([h0,h1])
hl.writeto(os.path.join(self.params['Output']['directory'],self.params['Output']['file'][0:-5] + '_noiseless_static_targets_ramp.fits'),overwrite=True)
#ILLUMINATION FLAT
if self.runStep['illuminationflat']:
illuminationflat,illuminationflatheader = self.readCalFile(self.params['Reffiles']['illumflat'])
signalramp *= illuminationflat
#signalzero *= illuminationflat
#PIXEL FLAT
if self.runStep['pixelflat']:
pixelflat,pixelflatheader = self.readCalFile(self.params['Reffiles']['pixelflat'])
signalramp *= pixelflat
#signalzero *= pixelflat
#IPC EFFECTS
if self.runStep['ipc']:
#ipcimage,ipcimageheader = self.readCalFile(self.params['Reffiles']['ipc'])
ipcimage = fits.getdata(self.params['Reffiles']['ipc'])
#Assume that the IPC kernel is designed for the removal of IPC, which means we need to
#invert it.
if self.params['Reffiles']['invertIPC']:
print("Iverting IPC kernel prior to convolving with image")
yk,xk = ipcimage.shape
newkernel = 0. - ipcimage
newkernel[(yk-1)/2,(xk-1)/2] = 1. - (ipcimage[1,1]-np.sum(ipcimage))
ipcimage = newkernel
g,y,x = signalramp.shape
for group in xrange(g):
signalramp[group,:,:] = s1.fftconvolve(signalramp[group,:,:],ipcimage,mode="same")
#signalzero = s1.fftconvolve(signalzero,ipcimage,mode='same')
#CROSSTALK
if self.runStep['crosstalk']:
if self.params['Readout']['namp'] == 4:
xdet = self.detector[3:5].upper()
if xdet[1] == 'L':
xdet = xdet[0] + '5'
xtcoeffs = self.readCrossTalkFile(self.params['Reffiles']['crosstalk'],xdet)
#Only sources on the detector will create crosstalk. If signalimage is larger than full frame
#because we are creating a grism image, then extract the pixels corresponding to the actual
#detector, and only create crosstalk values for those.
gd,yd,xd = signalramp.shape
xs = 0
xe = xd
ys = 0
ye = yd
if self.params['Output']['grism_source_image']:
xs,xe,ys,ye = self.extractFromGrismImage(signalramp[0,:,:])
for group in xrange(g):
xtinput = signalramp[group,ys:ye,xs:xe]
xtimage = self.crossTalkImage(xtinput,xtcoeffs)
#Now add the crosstalk image to the signalimage
signalramp[group,ys:ye,xs:xe] += xtimage
#signalzero[ys:ye,xs:xe] += self.crossTalkImage(signalzero[ys:ye,xs:xe],xtcoeffs)
else:
print("Crosstalk calculation requested, but the chosen subarray is read out using only 1 amplifier.")
print("Therefore there will be no crosstalk. Skipping this step.")
#Pixel Area Map
if self.runStep['pixelAreaMap']:
pixAreaMap = self.simpleGetImage(self.params['Reffiles']['pixelAreaMap'])
#If we are making a grism direct image, we need to embed the true pixel area
#map in an array of the appropriate dimension, where any pixels outside the
#actual aperture are set to 1.0
if self.params['Output']['grism_source_image']:
mapshape = pixAreaMap.shape
g,yd,xd = signalramp.shape
pam = np.ones((yd,xd))
ys = self.coord_adjust['yoffset']
xs = self.coord_adjust['xoffset']
pam[ys:ys+mapshape[0],xs:xs+mapshape[1]] = np.copy(pixAreaMap)
pixAreaMap = pam
signalramp *= pixAreaMap
#signalzero *= pixAreaMap
#save the combined static+moving targets ramp
if self.params['Output']['save_intermediates']:
h0=fits.PrimaryHDU()
h1 = fits.ImageHDU(signalramp)
hl=fits.HDUList([h0,h1])
hl.writeto(os.path.join(self.params['Output']['directory'],self.params['Output']['file'][0:-5] + '_noiseless_sources_plus_det_effects_ramp.fits'),overwrite=True)
#If grism direct data is requested, then save the signal
if self.params['Output']['grism_source_image']:
if ((self.params['Inst']['mode'].lower() == 'moving_target') or (mov_targs_integration is not None)):
#rearrange the RAPID mov_targs integration to the requested readout pattern
gd_integration = signalramp
if self.params['Readout']['readpatt'].upper() != 'RAPID':
gd_integration = self.changeReadPattern(gd_integration,self.params['Readout']['nframe'],self.params['Readout']['nskip'],self.params['Readout']['ngroup'])
self.createGrismDirectData(gd_integration)
else:
self.createGrismDirectData(frameimage)
#If we continue after saving the grism direct image, then
#we need to chop down signalimage to the nominal aperture
#size.
if not self.params['Output']['grism_input_only']:
xs,xe,ys,ye = self.extractFromGrismImage(signalramp[0,:,:])
signalramp = signalramp[:,ys:ye,xs:xe]
#signalzero = signalzero[ys:ye,xs:xe]
print('Simulated signal data cropped back down to nominal shape ({}), and proceeding.'.format(signalramp.shape))
else:
print("grism_input_only set to True in {}. Quitting.".format(self.paramfile))
return 0
#above here, signalramp is always RAPID, and always the full ramp, even if there
#are no moving targets.
#read in cosmic ray library files if CRs are to be added to the data later
if self.runStep['cosmicray']:
self.readCRFiles()
#read in saturation map
if self.runStep['saturation_lin_limit']:
try:
self.satmap,self.satheader = self.readCalFile(self.params['Reffiles']['saturation'])
bad = ~np.isfinite(self.satmap)
self.satmap[bad] = 1.e6
except:
print('WARNING: unable to open saturation file {}.'.format(self.params['Reffiles']['saturation']))
print("Please provide a valid file, or place 'none' in the saturation entry in the parameter file,")
print("in which case the nonlin limit value in the parameter file ({}) will be used for all pixels.".format(self.params['nonlin']['limit']))
sys.exit()
else:
print('CAUTION: no saturation map provided. Using {} for all pixels.'.format(self.params['nonlin']['limit']))
dy,dx = self.dark.data.shape[2:]
self.satmap = np.zeros((dy,dx)) + self.params['nonlin']['limit']
# Combine the dark current ramp and the synthetic signal ramp
print(signalramp.shape)
print(self.linDark.data.shape)
print(self.linDark.sbAndRefpix.shape)
print(zeroModel.data.shape)
print(zeroModel.sbAndRefpix.shape)
#print("stopping at doProperCombine. keep an eye on shapes of zeromodel. keep as 2d always for now? (assumes no multiple integrations. make sure this is true both for pipeline and non-pipeline situations")
proper_ramp,proper_zero = self.doProperCombine(signalramp,self.linDark.sbAndRefpix,zeroModel.data,zeroModel.sbAndRefpix)
final_ramp = proper_ramp
final_zero = proper_zero
#Set any pixels with signals above 65535 to 65535. We want to prevent crazy values from
#popping up when the ramp is saved. When being saved in DMS format, the ramp is converted
#to a 16-bit integer, so we don't want to feed in really large values to that conversion.
#Similarly, set any large negative values to zero. These values most likely
#come from doNonLin, where a bad initial guess or bad linearity coefficients
#caused the solution to fail to converge
if final_zero is not None:
bad = final_zero > 65535
final_zero[bad] = 65535
bad2 = final_zero < 0
final_zero[bad2] = 0
badhigh = final_ramp > 65535
badlow = final_ramp < 0
final_ramp[badhigh] = 65535
final_ramp[badlow] = 0
#save the integration
#if self.params['Output']['format'].upper() == 'DMS':
if self.params['Inst']['use_JWST_pipeline']:
print("Saving exposure in DMS format using datamodels.")
self.saveDMS(final_ramp,final_zero,self.params['Output']['file'])
else:
print("Saving exposure in DMS format using astropy.")
self.savefits(final_ramp,final_zero,os.path.join(self.params['Output']['directory'],self.params['Output']['file']))
#query the engineering database in order to popularte basic WCS info
set_telescope_pointing.add_wcs(os.path.join(self.params['Output']['directory'],self.params['Output']['file']),roll=self.params['Telescope']['rotation'])
def fullPaths(self):
# Expand all input paths to be full paths
# This is to allow for easier Condor-ization of
# many runs
pathdict = {'Reffiles':['dark','linearized_darkfile','hotpixmask','superbias',
'subarray_defs','readpattdefs','linearity',
'saturation','gain','phot','pixelflat',
'illumflat','astrometric','distortion_coeffs','ipc',
'crosstalk','occult','filtpupilcombo','pixelAreaMap',
'flux_cal'],
'cosmicRay':['path'],
'simSignals':['pointsource','psfpath','galaxyListFile','extended',
'movingTargetList','movingTargetSersic',
'movingTargetExtended','movingTargetToTrack'],
'newRamp':['dq_configfile','sat_configfile','superbias_configfile',
'refpix_configfile','linear_configfile'],
'Output':['file','directory']}
for key1 in pathdict:
for key2 in pathdict[key1]:
if self.params[key1][key2].lower() != 'none':
self.params[key1][key2] = os.path.abspath(self.params[key1][key2])
def getAttitudeMatrix(self):
#create an attitude_matrix from the distortion reference file model and other info
#calculate a local roll angle for the aperture
self.local_roll = set_telescope_pointing.compute_local_roll(self.params['Telescope']['rotation'],self.ra,self.dec,self.v2_ref,self.v3_ref)
#create attitude_matrix
attitude_matrix = rotations.attitude(self.refpix_pos['v2'],self.refpix_pos['v3'],self.ra,self.dec,self.local_roll)
return attitude_matrix
def darkints(self):
#check the number of integrations in the dark current file and compare with the
#requested number of integrations. Add/remove integrations if necessary
ndarkints,ngroups,ny,nx = self.dark.data.shape
reqints = self.params['Readout']['nint']
if reqints <= ndarkints:
#Fewer requested integrations than are in the input dark.
print("{} output integrations requested.".format(reqints))
self.dark.data = self.dark.data[0:reqints,:,:,:]
if self.dark.sbAndRefpix is not None:
self.dark.sbAndRefpix = self.dark.sbAndRefpix[0:reqints,:,:,:]
if self.dark.zeroframe is not None:
self.dark.zeroframe = self.dark.zeroframe[0:reqints,:,:]
elif reqints > ndarkints:
#More requested integrations than are in input dark.
print('Requested output has {} integrations, while input dark has only {}.'.format(reqints,ndarkints))
print('Adding integrations to input dark by making copies of the input.')
self.integration_copy(reqints,ndarkints)
def integration_copy(self,req,darkint):
#use copies of integrations in the dark current input to
#make up integrations in the case where the output has
#more integrations than the input
ncopies = np.int((req - darkint) / darkint)
extras = (req - darkint) % darkint
#full copies of the dark exposure (all integrations)
if ncopies > 0:
copydata = self.dark.data
if self.dark.sbAndRefpix is not None:
copysb = self.dark.sbAndRefpix
if self.dark.zeroframe is not None:
copyzero = self.dark.zero
for i in range(ncopies):
self.dark.data = np.vstack((self.dark.data,copydata))
if self.dark.sbAndRefpix is not None:
self.dark.sbAndRefpix = np.vstack((self.dark.sbAndRefpix,copysb))
if self.dark.zeroframe is not None:
self.dark.zeroframe = np.vstack((self.dark.zeroframe,copyzero))
#partial copy of dark (some integrations)
if extras > 0:
self.dark.data = np.vstack((self.dark.data,self.dark.data[0:extras,:,:,:]))
if self.dark.sbAndRefpix is not None:
self.dark.sbAndRefpix = np.vstack((self.dark.sbAndRefpix,self.dark.sbAndRefpix[0:extras,:,:,:]))
if self.dark.zeroframe is not None:
self.dark.zeroframe = np.vstack((self.dark.zeroframe,self.dark.zeroframe[0:extras,:,:]))
self.dark.header['NINTS'] = req
def nonsidereal_CRImage(self,file):
#create countrate image of non-sidereal sources that are being tracked.
#get countrate for mag 15 source, for scaling later
try:
if self.params['Readout']['pupil'][0].upper() == 'F':
usephot = 'pupil'
else:
usephot = 'filter'
mag15rate = self.countvalues[self.params['Readout'][usephot]]
except:
print("Unable to find mag 15 countrate for {} filter in {}.".format(self.params['Readout'][usephot],self.params['Reffiles']['phot']))
sys.exit()
totalCRList = []
#read in file containing targets
targs,pixFlag,velFlag = self.readMTFile(self.params['simSignals']['movingTargetToTrack'])
#We need to keep track of the proper motion of the target being tracked, because
#all other objects in the field of view will be moving at the same rate in the
#opposite direction
track_ra_vel = targs[0]['x_or_RA_velocity']
track_dec_vel = targs[0]['y_or_Dec_velocity']
#sort the targets by whether they are point sources, galaxies, extended
ptsrc_rows = []
galaxy_rows = []
extended_rows = []
for i,line in enumerate(targs):
if 'point' in line['object'].lower():
ptsrc_rows.append(i)
elif 'sersic' in line['object'].lower():
galaxy_rows.append(i)
else:
extended_rows.append(i)
#then re-use functions for the sidereal tracking situation
if len(ptsrc_rows) > 0:
ptsrc = targs[ptsrc_rows]
if pixFlag:
meta0 = 'position_pixel'
else:
meta0 = ''
if velFlag:
meta1 = 'velocity_pixel'
else:
meta1 = ''
meta2 = 'Point sources with non-sidereal tracking. File produced by ramp_simulator.py'
meta3 = 'from run using non-sidereal moving target list {}.'.format(self.params['simSignals']['movingTargetToTrack'])
ptsrc.meta['comments'] = [meta0,meta1,meta2,meta3]
ptsrc.write(os.path.join(self.params['Output']['directory'],'temp_non_sidereal_point_sources.list'),format='ascii',overwrite=True)
ptsrc = self.getPointSourceList(os.path.join(self.params['Output']['directory'],'temp_non_sidereal_point_sources.list'))
ptsrcCRImage = self.makePointSourceImage(ptsrc)
totalCRList.append(ptsrcCRImage)
if len(galaxy_rows) > 0:
galaxies = targs[galaxy_rows]
if pixFlag:
meta0 = 'position_pixel'
else:
meta0 = ''
if velFlag:
meta1 = 'velocity_pixel'
else:
meta1 = ''
meta2 = 'Galaxies (2d sersic profiles) with non-sidereal tracking. File produced by ramp_simulator.py'
meta3 = 'from run using non-sidereal moving target list {}.'.format(self.params['simSignals']['movingTargetToTrack'])
galaxies.meta['comments'] = [meta0,meta1,meta2,meta3]
galaxies.write(os.path.join(self.params['Output']['directory'],'temp_non_sidereal_sersic_sources.list'),format='ascii',overwrite=True)
#read in the PSF file with centered source
wfe = self.params['simSignals']['psfwfe']
wfegroup = self.params['simSignals']['psfwfegroup']
basename = self.params['simSignals']['psfbasename'] + '_'
if wfe == 0:
psfname=basename+self.params['simSignals'][usefilt].lower()+'_zero'
else:
psfname=basename+self.params['Readout'][usefilt].lower()+"_"+str(wfe)+"_"+str(wfegroup)
centerpsffile = self.params['simSignals']['psfpath'] + psfname + '_0p0_0p0.fits'
try:
centerpsf = fits.getdata(centerpsffile)
except:
print('Unable to read in PSF file for non sidereal moving target sersic work. {}'.format(centerpsffile))
sys.exit()
#crop the psf such that it is centered in its array
centerpsf = self.cropPSF(centerpsf)
#normalize the PSF to a total signal of 1.0
totalsignal = np.sum(centerpsf)
centerpsf = centerpsf / totalsignal
galaxyCRImage = self.makeGalaxyImage('temp_non_sidereal_sersic_sources.list',centerpsf)
totalCRList.append(galaxyCRImage)
if len(extended_rows) > 0:
extended = targs[extended_rows]
if pixFlag:
meta0 = 'position_pixel'
else:
meta0 = ''
if velFlag:
meta1 = 'velocity_pixel'
else:
meta1 = ''
meta2 = 'Extended sources with non-sidereal tracking. File produced by ramp_simulator.py'
meta3 = 'from run using non-sidereal moving target list {}.'.format(self.params['simSignals']['movingTargetToTrack'])
extended.meta['comments'] = [meta0,meta1,meta2,meta3]
extended.write(os.path.join(self.params['Output']['directory'],'temp_non_sidereal_extended_sources.list'),format='ascii',overwrite=True)
#read in the PSF file with centered source
wfe = self.params['simSignals']['psfwfe']
wfegroup = self.params['simSignals']['psfwfegroup']
basename = self.params['simSignals']['psfbasename'] + '_'
if wfe == 0:
psfname=basename+self.params['simSignals'][usefilt].lower()+'_zero'
else:
psfname=basename+self.params['Readout'][usefilt].lower()+"_"+str(wfe)+"_"+str(wfegroup)
centerpsffile = self.params['simSignals']['psfpath'] + psfname + '_0p0_0p0.fits'
try:
centerpsf = fits.getdata(centerpsffile)
except:
print('Unable to read in PSF file for non sidereal moving target sersic work. {}'.format(centerpsffile))
sys.exit()
#crop the psf such that it is centered in its array
centerpsf = self.cropPSF(centerpsf)
#normalize the PSF to a total signal of 1.0
totalsignal = np.sum(centerpsf)
centerpsf = centerpsf / totalsignal
#ext_src, extpixflag = self.readPointSourceFile(self.params['simSignals']['extended'])
extlist,extstamps = self.getExtendedSourceList('temp_non_sidereal_extended_sources.list')
#translate the extended source list into an image
extCRImage = self.makeExtendedSourceImage(extlist,extstamps)
#if requested, convolve the stamp images with the NIRCam PSF
if self.params['simSignals']['PSFConvolveExtended']:
extCRImage = s1.fftconvolve(extCRImage,centerpsf,mode='same')
totalCRList.append(extCRImage)
#Now combine into a final countrate image of non-sidereal sources (that are being tracked)
if len(totalCRList) > 0:
totalCRImage = totalCRList[0]
if len(totalCRList) > 1:
for i in xrange(1,len(totalCRList)):
totalCRImage += totalCRList[i]
else:
print("No non-sidereal countrate targets produced.")
print("You shouldn't be here.")
sys.exit()
return totalCRImage,track_ra_vel,track_dec_vel,velFlag
def combineSimulatedDataSources(self,inputtype,input1,input1_zero,mov_tar_ramp):
#inputtype can be 'countrate' in which case input needs to be made
#into a ramp before combining with mov_tar_ramp, or 'ramp' in which
#case you can combine directly. Use 'ramp' with
#moving_target MODE data, and 'countrate' with imaging MODE data
#Combine the countrate image with the moving target ramp.
if inputtype == 'countrate':
#First change the countrate image into a ramp
yd,xd = input1.shape
num_frames = self.params['Readout']['ngroup'] * (self.params['Readout']['nframe'] + self.params['Readout']['nskip'])
print("Countrate image of synthetic signals being converted to RAPID integration with {} frames.".format(num_frames))
input1_ramp = np.zeros((num_frames,yd,xd))
for i in xrange(num_frames):
input1_ramp[i,:,:] = input1 * self.frametime * (i+1)
else:
#if input1 is a ramp rather than a countrate image
input1_ramp = input1
#combine the input1 ramp and the moving target ramp, which are
#now in the same readout pattern
totalinput = input1_ramp + mov_tar_ramp
#totalzero = input1_zero + mov_tar_zero
return totalinput#,totalzero
def createGrismDirectData(self,input):
#Create the grism direct image or ramp to be saved
#Get the photflambda and photfnu values that go with
#the filter
try:
module = fits.getval(self.params['Reffiles']['linearized_darkfile'],'module')
except:
module = fits.getval(self.params['Reffiles']['dark'],'module')
zps = ascii.read(self.params['Reffiles']['flux_cal'])
if self.params['Readout']['pupil'][0] == 'F':
usephot = 'pupil'
else:
usephot = 'filter'
mtch = ((zps['Filter'] == self.params['Readout'][usephot]) & (zps['Module'] == module))
photflam = zps['PHOTFLAM'][mtch][0]
photfnu = zps['PHOTFNU'][mtch][0]
pivot = zps['Pivot_wave'][mtch][0]
arrayshape = input.shape
if len(arrayshape) == 2:
#make 3D for easier coding below
#input = np.expand_dims(input,axis=0)
units = 'e-/sec'
#indim = 2
yd,xd = arrayshape
tgroup = 0.
elif len(arrayshape) == 3:
units = 'e-'
#indim = 3
g,yd,xd = arrayshape
tgroup = self.frametime*(self.params['Readout']['nframe']+self.params['Readout']['nskip'])
grismDirectName = os.path.join(self.params['Output']['directory'],self.params['Output']['file'][0:-5] + '_' + self.params['Readout']['filter'] + '_GrismDirectData.fits')
xcent_fov = xd / 2
ycent_fov = yd / 2
kw = {}
kw['xcenter'] = xcent_fov
kw['ycenter'] = ycent_fov
kw['units'] = units
kw['TGROUP'] = tgroup
if self.params['Readout']['pupil'][0].upper() == 'F':
usefilt = 'pupil'
else:
usefilt = 'filter'
kw['filter'] = self.params['Readout'][usefilt]
kw['PHOTFLAM'] = photflam
kw['PHOTFNU'] = photfnu
kw['PHOTPLAM'] = pivot * 1.e4 #put into angstroms
self.saveSingleFits(input,grismDirectName,key_dict=kw)
print("Data to be used as input to make dispersed grism image(s) saved as {}".format(grismDirectName))
def changeReadPattern(self,array,outnframe,outnskip,outngroup):
#Assume an input integration with a RAPID readout pattern
#Average/skip frames to convert the integration to another pattern
ingroup,iny,inx = array.shape
newpatt = np.zeros((outngroup,iny,inx))
#total number of frames per group
deltaframe = outnskip + outnframe
#Loop over groups
for i in range(outngroup):
#average together the appropriate frames, skip the appropriate frames
frames = np.arange(i*deltaframe,i*deltaframe+outnframe)
#If averaging needs to be done
if outnframe > 1:
newpatt[i,:,:] = np.mean(array[frames,:,:],axis=0)
print("In changereadpattern, averaging frames {}".format(frames))
#If no averaging needs to be done
else:
newpatt[i,:,:] = array[frames[0],:,:]
return newpatt