-
Notifications
You must be signed in to change notification settings - Fork 2
/
procedures.py
9501 lines (8941 loc) · 626 KB
/
procedures.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/python
# -*- coding: utf-8 -*-
import os
os.environ['KMP_DUPLICATE_LIB_OK']='True' # to avoid "OMP: Error #15: Initializing libiomp5md.dll, but found libiomp5md.dll already initialized." when curve_fit is executed in windows
import sys
import time
import datetime
import operator
import copy
import random
import warnings
import numpy as np
import numpy.typing as npt
from typing import Union
from astropy.io import fits
from astropy.table import Table
import astropy.time as asttime
import astropy.coordinates as astcoords
import astropy.units as astunits
from astropy.utils.exceptions import AstropyUserWarning
#import astropy.constants as astconst
if sys.version_info[0] < 3: # Python 2
import matplotlib
matplotlib.use('agg') # Otherwise keeps hanging when using CERES; GUIs work, but plt.show() won't; matplotlib.interactive(True) doesn't show
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
from scipy.optimize import curve_fit # Curve-fit Gauss takes the longest: 20ms
# import scipy.interpolate as inter # Not used anymore
import scipy
import scipy.signal
try:
from tqdm import tqdm
except:
print('Error: tqdm could not be loaded. Did you activate the Anaconda hiflex environment?')
exit(1)
import json
# detect python version
if sys.version_info[0] < 3: # Python 2
from collections import OrderedDict as dict
import urllib2
else: # Python 3
import urllib
raw_input = input
if sys.version_info[0] == 3 and sys.version_info[1] >= 5:
try:
from deepCR import deepCR # Only available for python 3.5
except:
print('Warn: deepCR could not be loaded. Did you activate the Anaconda hiflex environment?')
deepCR = None
else:
deepCR = None
import plot_img_spec
import psutil
try:
import ephem
except:
print('Error: ephem could not be loaded. Did you activate the Anaconda hiflex environment?')
exit(1)
import math
"""try:
import multiprocess as multiprocessing # For Windows
print('Info: using multiprocess')
except:
import multiprocessing # Won't work under windows and Anaconda
print('Info: using multiprocessing')"""
import multiprocessing # Won't work under windows and Anaconda
import subprocess
# Necessary because of https://github.com/astropy/astropy/issues/9427 # commented out again on 20/5/2020 after "WARNING: IERSStaleWarning: IERS_Auto predictive values are older than 15.0 days but downloading the latest table did not find newer values [astropy.utils.iers.iers]"
try:
import astropy
except:
print('Error: astropy could not be loaded. Did you activate the Anaconda hiflex environment?')
exit(1)
""" # Stopped working on 24 May 2021
from astropy.utils.iers import conf as iers_conf
#iers_conf.iers_auto_url = 'https://astroconda.org/aux/astropy_mirror/iers_a_1/finals2000A.all' # is too old as of May 2020
iers_conf.iers_auto_url = 'https://datacenter.iers.org/data/9/finals2000A.all' # untested
iers_conf.iers_auto_url = 'ftp://cddis.gsfc.nasa.gov/pub/products/iers/finals2000A.all' # worked on 25 May 2020
iers_conf.auto_max_age = None """
success = False
for ii in range(5):
if sys.version_info[0] < 3: # Python 2
try:
import barycorrpy
success = True
break
except (urllib2.URLError, ValueError, astropy.utils.iers.iers.IERSRangeError) as e:
print('Warn: Cannot import barrycorrpy. Will try {0} more times. Error: {1}, Reason: {2}'.format(4-ii, e, e.reason))
try:
print('Warn: Cannot import barrycorrpy. Will try {0} more times. Error: {1}, Reason: {2}'.format(4-ii, e, e.reason))
except:
print('Warn: Cannot import barrycorrpy. Will try {0} more times. Error: {1}'.format(4-ii, e))
else:
try:
import barycorrpy
success = True
break
except (urllib.error.URLError, ValueError, astropy.utils.iers.iers.IERSRangeError) as e:
print('Warn: Cannot import barrycorrpy. Will try {0} more times. Error: {1}, Reason: {2}'.format(4-ii, e, e.reason))
if not success:
print('Error: barrycorrpy could not be loaded. It needs an active internet connection in order to download the IERS_B file. This failure will lead to a crash of the program later!'+os.linesep)
import glob
import pickle
import platform
import proc_gui
from proc_general import *
""" This was thought to be a way to improve calculation speed when calculating gaussians. However, even with declearing all variables it didn't speed up calculations
try:
import cy_oneD_gauss # Compiled code to calculate the Gaussian
except:
cy_oneD_gauss = None
"""
""" only needed for BJD calculation using jplephem and BVC calculation from CERES pipeline
import jplephem # jplehem from the pip install jplehem
import de423
# SSEphem package http://www.cv.nrao.edu/~rfisher/Python/py_solar_system.html coupled to SOFA http://www.iausofa.org/
# to clean up: mv SSEphem/ ssephem_update.py old/
baryc_dir = os.path.dirname(os.path.abspath(__file__))+'/SSEphem/'
sys.path.append(baryc_dir)
if not os.path.isfile(baryc_dir+'man_jplephem.py'):
os.system('cd {0} && ln jplephem.py man_jplephem.py'.format(baryc_dir))
time.sleep(1)
if not os.path.isfile(baryc_dir+'man_jplephem.py'):
print('Please run in a different termina the following line and then press Enter')
print('cd {0} && ln jplephem.py man_jplephem.py'.format(baryc_dir))
raw_input('')
import man_jplephem # jplephem in the SSEphem folder
ephemeris = 'DEc403' # To use the JPL DE403 ephemerides, https://en.wikipedia.org/wiki/Jet_Propulsion_Laboratory_Development_Ephemeris
"""
tqdm.monitor_interval = 0 #On the virtual machine at NARIT the code raises an exception otherwise
calimages = dict() # dictionary for all calibration images used by create_image_general and read_file_calibration
beg_ts = time.time()
GLOBALutils, correlation, lowess = None, None, None
class Constants: # from CERES
"Here I declare the constants I will use in the different functions"
"G,c: the gravity and speed of light constants in SI units; mearth and mmoon: the earth and moon masses in kg"
G,c = 6.673E-11,2.99792458E8
mearth, mmoon = 5.9736E24,7.349E22
"the mass of the planets,moon and sun in earth masses, the ennumeration is sun, moon, mercury,venus,mars,jupiter,saturn, uranus,neptune,pluto"
mplanets = np.array([332981.787,0.0123000371,0.05528,0.815,0.10745,317.83,95.19,14.536,17.147,0.0021])
"conversion from degrees to radians, and from hours to degrees, and from degrees to hours"
degtorad = np.pi/180.0
radtodeg = 180.0/np.pi
HtoDeg = 360.0/24.0
DegtoH = 1.0/HtoDeg
"Req: equatorial radius f: polar flatenning , w angular speedof earth in rad/s"
Req = 6378136.6 #in m
f = 1.0/298.256420
w = 7.2921158554E-5
DaystoYear = 1.0/365.256363004
MJD0 = 2400000.5
def get_timing(text=''):
global beg_ts
end_ts = time.time()
print(text+"elapsed time: %f" % (end_ts - beg_ts))
beg_ts = time.time()
def time_usage(func):
def wrapper(*args, **kwargs):
beg_ts = time.time()
retval = func(*args, **kwargs)
end_ts = time.time()
print("elapsed time: %f" % (end_ts - beg_ts))
return retval
return wrapper
def logger(message, show=True, printarrayformat=[], printarray=[], logfile='logfile', params=None):
"""
Saves the status information to a logfile
:param message: Text to log
:param show: if False than it will be logged only to the logfile but not printed on the screen
:param printarrayformat: Format of the columns in printarray, e.g. '%3.1f', '%2.2i'
:param printarray: Array with values to log
:param logfile: filename in which the information will be written
"""
if show:
print(message)
with open(logfile, 'a') as file:
if multiprocessing.current_process().name == 'MainProcess': pid = os.getpid()
else: pid = '{0}-{1}'.format(os.getppid(), os.getpid())
file.write('{1} - {2} - {3}{0}'.format( os.linesep, time.strftime("%Y%m%d%H%M%S", time.localtime()), pid, message ))
if len(printarrayformat) >0 and len(printarray) > 0:
for line in printarray:
text = ''
for i,printformat in enumerate(printarrayformat):
#print printformat,line[i]
text += printformat%line[i] + '\t'
file.write(text[:-1]+os.linesep)
if show:
print(text)
if message.find('Error') == 0:
print('\t-> exiting')
#if 'params' in locals() or 'params' in globals(): # this line can't compile, try except doesn't work either
if params is not None:
log_params(params)
exit(1)
def log_params(params, start=False):
"""
formats the dictionary to be saved in the logfile
"""
# Finding the python files
if start:
text = ''
list_of_files = glob.glob('{0}/*.py'.format(os.path.realpath(__file__).rsplit(os.sep,1)[0]))
list_of_files = sorted(list_of_files, key=os.path.getmtime, reverse=True)
for fname in list_of_files:
filedata1 = read_text_file(fname, no_empty_lines=False)
filedata2 = read_text_file(fname, no_empty_lines=True)
text += os.linesep+' {1} {2} {3} {4} {0}'.format( fname,
datetime.datetime.utcfromtimestamp(os.path.getmtime(fname)).strftime('%Y-%m-%dT%H:%M:%S'),
'%10.1i'%os.stat(fname).st_size, '%6.1i'%len(filedata1), '%6.1i'%len(filedata2) )
logger('python files with unix timestamp of modification, size in bytes, number of lines, and number of non-empty lines: '+text, show=False, logfile='logfile_params', params=params)
logger('Info: Using procedures file {0}'.format( os.path.realpath(__file__) ), show=False, params=params)
# Log the versions of the packages
text = os.linesep+' {: <20s} {}'.format( platform.system(), platform.release() )
text += os.linesep+' {: <20s} {}'.format('python', sys.version)
import types, pkg_resources
name, val = '', ''
for name, val in globals().items():
if isinstance(val, types.ModuleType):
try:
text += os.linesep+' {: <20s} {}'.format(val.__name__ , pkg_resources.get_distribution(val.__name__).version )
except:
text += os.linesep+' {}'.format(val.__name__)
logger('Info: imported packages, and, if available, the version'+text, show=False, logfile='logfile_params')
# Logging the parameters
paramstxt = dict()
for key in params.keys():
if type(params[key]) in [list, np.ndarray]:
paramstxt[key] = str(params[key])
else:
paramstxt[key] = params[key]
logger('params: '+json.dumps(paramstxt, sort_keys = False, indent = 4), show=False, logfile='logfile_params', params=params)
# Old way as real json, but not very compact
# logger('params: '+json.dumps(params, sort_keys = False, indent = 4), show=False, logfile='logfile_params')
def read_parameterfile(textfile):
# load text file (remove all white spaces)
if not os.path.exists(textfile):
logger('Error: The parameterfile {0} does not exist.'.format(textfile))
""" # Extra steps to check that the user didn't make mistakes in the file:
data = read_text_file(textfile, no_empty_lines=True)
data = convert_readfile(data, [str,str], delimiter='=', replaces=[' ', '\t',['\\',' ']], ignorelines=['#']) # this replaces too many spaces
for line in data:
if len(line) != 2:
logger(('Error: The line in file {0} containing the entries {1} has the wrong format. Expected was "parameter = value(s)" .'+\
'Please check if the "=" sign is there or if a comment "#" is missing.').format(textfile, line))"""
try:
keys, values = np.genfromtxt(textfile, dtype=str, comments='#', delimiter='=', filling_values='', autostrip=True, unpack=True)
except ValueError as error:
print(error)
logger(('Error: One line (see previous output (empty lines are missing in the line number counting)) in file {0} has the wrong format. Expected was "parameter = value(s)" .'+\
'Please check if the "=" sign is there or if a comment "#" is missing.').format(textfile))
# raise # just this to show the error
# raise ValueError # Don't do this, you'll lose the stack trace!
if len(keys) == 0 or len(values) == 0:
logger('Error: No values found when reading the parameterfile {0}.'.format(textfile))
"""# This bit was necesary before adding autostrip and unpack
keys, values = data[:, 0], data[:, 1]
for k in range(len(keys)):
keys[k] = keys[k].replace(' ', '')
values[k] = values[k].replace(' ', '')"""
textparams = dict(zip(keys, values))
return textparams
def read_cmdparams():
cmdargs = sys.argv[1:]
cmdparams = dict()
for arg in cmdargs:
if '=' in arg:
key, value = arg.replace(' ', '').split('=')
cmdparams[key] = value
# From Neil, helpful to make the below list of formating shorter
#VARIABLES = dict(filepath=str, plotpath=str, savepath=str, files=list,
# order_direction=str, xblur=float, yblur=float,
# fitmin=int, fitmax=int, width=int, minpixelsinorder=int)
#if key in VARIABLES:
# # try to convert to required type (using VARIABLES dict)
# try:
# if VARIABLES[key] in [list, np.ndarray]:
# ckwargs[key] = VARIABLES[key](value.split(','))
# else:
# ckwargs[key] = VARIABLES[key](value)
# except ValueError:
# emsg = [key, str(VARIABLES[key])]
# print('Command line input not understood for' +
# 'argument {0} must be {1}'.format(*emsg))
elif arg in ['prepare'] or 0 <= arg.lower().find('nogui') <= 2 or 0 <= arg.lower().find('norv') <= 2:
continue # Do nothing, just prevent the warning below
else:
logger('Warn: I dont know how to handle command line argument: {0}'.format(arg))
return cmdparams
def textfileargs(params, textfile=None):
"""
Imports parameters from text file located at textfile, parameters given in the command line and the parameters given with the programm code
Parameters in the programm code will be overwritten by parameters in the text file and parameters in the text file will be overwritten by parameters from the command line
:param params: dictionary, default parameter dictionary (the default values)
:param textfile: string, location of config file
:return params: dictionary with added parameters and updated default parameters
Please note, the parameters are only lists or single entries. No numpy arrays allowed
"""
# Set up standard parameters for
params['in_shift'] = -0
params['started_from_p3'] = 'False'
emsg = 'Error in config file: '
textparams = read_parameterfile(textfile)
params.update(textparams)
cmdparams = read_cmdparams()
params.update(cmdparams)
if 'configfile_fitsfiles' in params.keys():
if os.path.exists(params['configfile_fitsfiles']):
textparams = read_parameterfile(params['configfile_fitsfiles'])
params.update(textparams)
params.update(cmdparams) # Repeat, in case textparams overwrites cmdparams
# Add standard parameters, if they are missing, more at the end
params['logging_arc_line_identification_residuals_hist'] = params.get('logging_arc_line_identification_residuals_hist', 'arc_line_identification_residuals_hist.png')
params['logging_em_lines_bisector'] = params.get('logging_em_lines_bisector', 'wavelength_solution_emmission_lines_bisector.png')
params['logging_blaze_spec'] = params.get('logging_blaze_spec', 'blaze_spectrum.pdf')
# Added 20201022:
params['logging_crossdispersion_shift'] = params.get('logging_crossdispersion_shift', 'crossdispersions_shift.txt')
# Added 20201210:
params['logging_resolution_form'] = params.get('logging_resolution_form', 'wavelength_solution_resolution_on_detector.png')
# Added 20210129:
params['log_wavesolshift_details'] = params.get('log_wavesolshift_details', 'True')
# Added 20210323:
params['convert_rv_package_templates'] = params.get('convert_rv_package_templates', 'True')
# Added 20210323:
params['split_wavelength_solutions_at'] = params.get('split_wavelength_solutions_at', '[]')
# Added 20210423:
params['logging_em_lines_shift'] = params.get('logging_em_lines_shift', 'wavelength_solution_emmission_lines_shift_on_detector.png')
# Added 20210626:
params['traces_min_separation'] = params.get('traces_min_separation', '5')
list_txt = ['reference_catalog', 'use_catalog_lines', 'raw_data_file_endings', 'raw_data_mid_exposure_keys', 'raw_data_paths', 'raw_data_object_name_keys', 'cosmic_ray_settings']
list_int = ['arcshift_range', 'order_offset', 'px_offset', 'polynom_order_traces', 'polynom_order_intertraces',
'bin_search_apertures', 'bin_adjust_apertures', 'polynom_bck', 'dataset_rv_analysis', 'split_wavelength_solutions_at']
list_float = ['opt_px_range', 'px_offset_order', 'background_width_multiplier', 'sigmaclip_spectrum', 'traces_searchlimit_brightness']
list_abs = ['arcshift_range']
ints = ['polynom_order_apertures', 'rotate_frame']
floats = ['max_good_value', 'catalog_file_wavelength_muliplier', 'extraction_width_multiplier', 'arcextraction_width_multiplier',
'resolution_offset_pct', 'diff_pxs', 'maxshift_orders', 'wavelength_scale_resolution', 'width_percentile', 'raw_data_timezone_cor',
'altitude', 'latitude', 'longitude', 'in_shift', 'extraction_shift', 'extraction_min_ident_part_of_trace_percent',
'max_cores_used_pct', 'pxshift_between_wavesolutions', 'minimum_SNR', 'traces_min_separation']
bools = ['flip_frame', 'update_width_orders', 'GUI', 'started_from_p3', 'log_wavesolshift_details', 'convert_rv_package_templates']
text_selection = ['arcshift_side', 'extraction_precision']
results = ['path_extraction', 'path_extraction_single', 'logging_path', 'path_reduced', 'path_rv_ceres', 'path_rv_terra', 'path_rv_serval',
'path_harpsformat', 'master_blaze_spec_norm_filename'] #, 'configfile_fitsfiles' (excluded, as should be handled as conf.txt
full_filenames = ['badpx_mask_filename', 'original_master_traces_filename', 'original_master_wavelensolution_filename',
'configfile_fitsfiles', 'raw_data_file_list', 'terra_jar_file'] # deal with full filenames -> nothing to do
texts = ['editor', 'extracted_bitpix', 'site', 'object_file', 'raw_data_imtyp_keyword',
'raw_data_imtyp_bias', 'raw_data_imtyp_dark', 'raw_data_imtyp_flat', 'raw_data_imtyp_trace1', 'raw_data_imtyp_blaze', 'raw_data_imtyp_trace2',
'raw_data_exptim_keyword', 'raw_data_dateobs_keyword', 'raw_data_timezone_cor', 'blazercor_function', 'px_to_wavelength_file'] # -> nothing to do
paths, loggings = [], []
#list_raw = [] # not necessary anymore
for entry in params.keys(): # make the live easier by adding entries automatically to the lists above
if entry not in list_int and (entry.find('subframe') >= 0): # add to the list
list_int.append(entry)
if entry not in list_txt and (entry.find('_rawfiles') > 0 or entry.find('calibs_') > -1): # add to the list
list_txt.append(entry)
#if entry not in list_raw and (entry.find('_rawfiles') > 0):
# list_raw.append(entry)
if entry not in paths and (entry.find('path_') != -1 or entry.find('_path') != -1):
paths.append(entry)
if entry not in results and ((entry.find('master_') == 0 or entry.find('background_') == 0) and entry.find('_filename') > 0): # Put these files also into the result path
results.append(entry)
if entry not in loggings and (entry.find('logging_') == 0 and entry.find('logging_path') == -1):
loggings.append(entry)
if entry not in texts and (entry.find('raw_data_imtyp_') == 0 or (entry.find('raw_data_') == 0 and entry.find('_keyword') > 0)):
texts.append(entry)
if entry not in list_txt and (type(params[entry]) == str):
if params[entry].find(',') != -1:
list_txt.append(entry)
all_parameters = list(np.unique( list_txt + list_int + list_float + list_abs + ints + floats + bools + text_selection + results + full_filenames + paths + loggings + texts )) # + list_raw
trues = ['yes', 'true', '1']
falses = ['no', 'false', '0']
lefts = ['left','l']
rights = ['right','r']
centers = ['center','centre','c']
standards = ['standard']
linfits = ['linfit']
undeclared_params = ''
# Important settings first, as other things depend on them:
for entry in paths:
if entry in params.keys() and entry not in list_txt:
params[entry] = (params[entry]+os.sep).replace(os.sep+os.sep, os.sep) # Add a / at the end in case the user didn't
for entry in results:
if entry in params.keys(): # deal with result filenames/folders -> add result_path
params[entry] = params['result_path'] + params[entry]
# All parameters
for entry in params.keys():
if entry in list_txt+list_int+list_float: # Split the lists
temp = params[entry]
for i in ['[', ']']: # remove the brakets
temp = temp.replace(i,'')
temp = [x.strip() for x in temp.split(',')] # strip the entries from leading or ending spaces
for i in range(len(temp))[::-1]: # remove empty entries
if temp[i] == '':
del temp[i]
params[entry] = temp
for [list_no, funct, functxt] in [ [list_int, int, 'intergers'], [list_float, float, 'floats'], [list_abs, abs, 'numbers'] ]:
if entry in list_no: # deal with list of integers, floats, abs, numpy array (int and float could be also done with dtype, but then needs to be converted back to list
for i in range(len(params[entry])):
try:
params[entry][i] = funct(params[entry][i])
except:
logger(emsg + 'Parameter "{0}" (value of "{1}") must be a list of {2}. Error occured with index {3}.'.format(entry, params[entry], functxt, i), params=params)
for [ints_floats, funct, functxt] in [ [ints, int, 'intergers'], [floats, float, 'floats'] ]:
if entry in ints_floats: # deal with integers or floats
try:
params[entry] = funct(params[entry])
except:
logger(emsg + 'Parameter "{0}" (value of "{1}") must be a {2}'.format(entry, params[entry], functxt), params=params)
if entry in bools:
if params[entry].lower() in trues: params[entry] = True
elif params[entry].lower() in falses: params[entry] = False
else:
logger(emsg + 'Parameter "{0}" (value of "{1}") must be any of the following: {2}'.format(entry, params[entry], trues+falses), params=params)
if entry in text_selection:
if entry == 'arcshift_side':
if params[entry].lower() in lefts: params[entry] = -1
elif params[entry].lower() in rights: params[entry] = 1
elif params[entry].lower() in centers: params[entry] = 0
else:
logger(emsg + 'Parameter "{0}" (value of "{1}") must be one of the following: {2}'.format(entry, params[entry], lefts+rights+centers), params=params)
if entry == 'extraction_precision':
if params[entry].lower() in standards: params[entry] = standards[0]
elif params[entry].lower() in linfits: params[entry] = linfits[0]
else:
logger(emsg + 'Parameter "{0}" (value of "{1}") must be one of the following: {2}'.format(entry, params[entry], standards+linfits), params=params)
if entry in paths: # Create folders
if type(params[entry]).__name__ in ['list']:
for ii in range(len(params[entry])):
params[entry][ii] = (params[entry][ii]+os.sep).replace(os.sep+os.sep, os.sep) # Add a / at the end in case the user didn't
if entry in ['raw_data_paths', 'path_ceres', 'path_serval']: # raw_data_paths is list, hence has to be checked before
continue
if params[entry].lower() not in ['na'+os.sep, params['result_path']+'na'+os.sep, params['result_path']+os.sep+'na'+os.sep]:
make_directory(params[entry]) # Create folders, if necessary
#if entry in list_raw: # deal with lists of raw data filenames -> add path
# for i in range(len(params[entry])):
# params[entry][i] = params['raw_data_path'] + params[entry][i] # Not necessary anymore
if entry in loggings: # deal with logging filenames/folders -> add logging_path
params[entry] = params['logging_path'] + params[entry]
if entry not in all_parameters:
undeclared_params += '{0}, '.format(entry)
overdeclared_params = ''
for entry in all_parameters:
if entry not in params.keys():
overdeclared_params += '{0}, '.format(entry)
if len(overdeclared_params) > 0:
logger('Warn: The following parameters are expected, but do not appear in the configuration files. : {1}{0}\t!!! The missing parameters might cause crashes later !!!'.format(os.linesep, overdeclared_params[:-2]), params=params)
if len(undeclared_params) > 0:
logger('Warn: The following parameters appear in the configuration files, but the programm did not expect them: {0}'.format(undeclared_params[:-2]), params=params)
# Use standard parameters, if not given iin configuration; more at the beginning, e.g. when modification of paths is necessary
params['use_cores'] = int(multiprocessing.cpu_count()*params.get('max_cores_used_pct',80)/100.) # Use 80% of the CPU cores
params['dataset_rv_analysis'] = params.get('dataset_rv_analysis', [5, 6])
params['pxshift_between_wavesolutions'] = params.get('pxshift_between_wavesolutions', 0)
params['cosmic_ray_settings'] = params.get('cosmic_ray_settings', ['deepCR', 'ACS-WFC-F606W-2-32', '0.999'])
params['traces_searchlimit_brightness'] = params.get('traces_searchlimit_brightness', 50)
params['minimum_SNR'] = params.get('minimum_SNR', 1.0)
return params
def make_directory(directory, errormsg='Warn: Cannot create directory {0}'):
if not os.path.exists(directory):
try: # Create folders, if necessary
os.makedirs(directory)
except:
logger(errormsg.format(directory))
def log_returncode(code, explain=''):
"""
Logs a message if a subprocess fails
"""
if code != 0 and code != None:
logger('Warn: Subprocess returned with error code {0}. {1}'.format(code, explain))
return code
#def update_calibration_memory(key, value): # Not needed anymore with global calimages in hiflex.py and in procedures here
# """
# Add new information to the global variable calimages, which is not accessable from other python files
# :param key: string, key for the dictionary
# :param value: string, number, array, or anything: value for the dictionary
# """
# global calimages
# calimages[key] = copy.deepcopy(value)
def sort_for_multiproc_map(inputlist, number_cores):
"""
multiprocessing.Pool.map will give the first part of the list to the first subprocess, second part to second subprocess and so on.
So in order to get data reduced in a more linear way, resorting the list will be helpful
:param inputlist: list of anthing
:param number_cores: integer
:return outlist: same length as inputlist, only resorted (e.g. number_cores=3: inputlist=[1,'a',3,4,5,6,7,8] -> outlist=[1,4,7,'a',5,8,3,6]
"""
return inputlist # Do nothing
"""# !!! Doesn't work as expected
outlist = []
for ii in range(number_cores):
outlist.append([]) # [[]]*number will create empty arrays that are all the same, so if one is filled, also the others are filled
# Resort
for ii in range(len(inputlist)):
outlist[ii%number_cores].append(inputlist[ii])
# Combine into one list again
for entry in outlist[1:]:
outlist[0] += entry
outlist = outlist[0]
return outlist"""
def read_text_file(filename, no_empty_lines=False, warn_missing_file=True):
"""
Read a textfile and put it into a list, one entry per line
Empty lines means they consist only of tabs and spaces
"""
text = []
if os.path.isfile(filename) == True:
text1 = open(filename,'r').readlines()
for line in text1:
line = line.replace(os.linesep, '')
linetemp = line.replace('\t', '')
linetemp = linetemp.replace(' ', '')
if ( line == '' or linetemp == '') and no_empty_lines:
continue
text.append(line)
elif warn_missing_file:
logger('Warn: File {0} does not exist, assuming empty file'.format(filename))
return text
def add_text_file(text, filename):
"""
Adds a line or lines of text to a file without checking if the line already exists
"""
with open(filename, 'a') as file:
file.write(text+os.linesep)
def add_text_to_file(text, filename, warn_missing_file=True):
"""
If the text is not yet in the file, the text is added to the file
"""
oldtext = read_text_file(filename, warn_missing_file=warn_missing_file)
exists = False
for line in oldtext:
if line.find(text) != -1:
exists = True
break
if exists == False:
add_text_file(text, filename)
def convert_readfile(input_list, textformats, delimiter='\t', replaces=[], ignorelines=[], expand_input=False, shorten_input=False, ignore_badlines=False, replacewithnan=False):
"""
Can be used convert a read table into entries with the correct format. E.g integers, floats
Ignories the lines which have less entries than entries in textformats
replacements will be done before
:param input_list: 1D list or array of strings from a read file
:param textformats: 1D list or array of formats, e.g. [str, str, int, float, float].
Individual entries can consist of a list of entries, then the conversion will be done in the order given, e.g. [str, [float, int], [float, '%Y-%m-%dT%H:%M:%S.%f'], float]
If a string is given this will be used to convert it into a datetime object
If a float, int, or str is run on a datetime object, then the datetime object will be converted into a number
:param delimiter: string, used to split each line into the eleements
:param replaces: 1D list or array of strings and/or lists of 2 strings, contains the strings which should be replaced by '' (if string) or by the second entry (if list)
:param ignorelines: List of strings and/or lists. A list within ignorelines needs to consist of a string and the maximum position this string can be found.
If the string is found before the position, the entry of the input list is ignored
:param expand_input: bool, if the line in the input_list contains less elements than textformats the line will be filled up with ''s
:param shorten_input: bool, if the line in the input_list contains more elements than textformats the line will be shortened to len(textformats)
:param ignore_badlines: bool, if True will raise Warnings, if False will raise Errors (which will terminate code)
:param replacewithnan: bool, if True and conversion with textformats is not possible will replace with NaN, otherwise will raise Warning/Error
:retrun result_list: 2D list with numbers or strings, formated acording to textformats
"""
# Make the textformats, replaces, and ignorelines consistent
for index in range(len(textformats)):
if type(textformats[index]) != list:
textformats[index] = [textformats[index]] # convert the single entries into list, e.g. make a list of lists
for index in range(len(ignorelines)):
error = False
if type(ignorelines[index]) == str: # only one entry
ignorelines[index] = [ignorelines[index], 1E10]
elif type(ignorelines[index]) == list: # list with two entries: string, and position after which the string will be ignored
if len(ignorelines[index]) != 2:
error = True
else:
error = True
if error:
logger('Error: Programming error: ignorelines which where hand over to convert_readfile are wrong. '+\
'It has to be a list consisting of strings and/or 2-entry lists of string and integer. Please check ignorelines: {0}'.format(ignorelines))
for index in range(len(replaces)):
error = False
if type(replaces[index]) == str: # only one entry
replaces[index] = [replaces[index], '']
elif type(replaces[index]) == list: # list with two entries: search-string and replace-string
if len(replaces[index]) != 2:
error = True
else:
error = True
if error:
logger('Error: Programming error: replaces which where hand over to convert_readfile are wrong. '+\
'It has to be a list consisting of strings and/or 2-entry lists of strings. Please check replaces: {0}'.format(replaces))
# Convert the text
error = {False:'Error:',True:'Warn:'}[ignore_badlines]
result_list = []
for entry in input_list:
notuse = False
for ignore in ignorelines:
if entry.find(ignore[0]) <= ignore[1] and entry.find(ignore[0]) != -1:
notuse = True
break
if notuse:
continue
for replce in replaces:
entry = entry.replace(replce[0],replce[1])
entry = entry.split(delimiter)
if len(entry) < len(textformats): # add extra fields, if not enough
if expand_input:
entry += [''] * ( len(textformats) - len(entry) )
else:
continue
for index in range(len(textformats)):
for textformat in textformats[index]:
if type(textformat) == type:
if type(entry[index]) == datetime.datetime:
epoch = datetime.datetime.utcfromtimestamp(0)
entry[index] = (entry[index] - epoch).total_seconds() # (obsdate - epoch) is a timedelta
try:
entry[index] = textformat(entry[index])
except:
if replacewithnan:
entry[index] = np.nan
else:
logger(error+' Cannot convert {0} into a {1}. The problem happens in line {2} at index {3}.'.format(entry[index], textformat, entry, index))
elif type(textformat) == str:
try:
entry[index] = datetime.datetime.strptime(entry[index], textformat)
except:
logger(error+' Cannot convert {0} into a datetime object using the format {1}. The problem happens in list line {2} at index {3}.'.format(entry[index], textformat, entry, index))
else:
logger('Error: Programming error: textformats which where hand over to convert_readfile are wrong. It has to be a type or a string')
if shorten_input and len(entry) > len(textformats):
del entry[len(textformats):]
result_list.append(entry)
return result_list
def read_badpx_mask(params, imageshape):
"""
Reads the bad-pixel-mask and applies the corrections to it
:param params: Dictionary with all the parameters:
badpx_mask_filename is used in order to get the bad px mask. If it doesn't exist, than all pixels are fine
calibs is used in order to if a subframe of the image is used. In this case subframe is used to find the area used
:param imageshape: tuple of two integers, gives the shape of the image to which the badpixel mask should be applied
:return badpx_mask: 2D numpy array of the bad pixel mask
"""
filename = params['badpx_mask_filename']
if os.path.isfile(filename) == True:
badpx_mask = np.array(fits.getdata(filename), dtype=np.float64)
badpx_mask = rotate_flip_frame(badpx_mask, params )
if badpx_mask.shape != imageshape:
if 'subframe' in params.keys():
subframe = params['subframe']
if len(subframe) >= 4:
badpx_mask = badpx_mask[subframe[2]: subframe[0]+subframe[2], subframe[3]: subframe[1]+subframe[3]]
if badpx_mask.shape != imageshape:
logger('Warn: The loaded badpixel mask {0} has a different format than the image file: {1} vs {2} (image file)'.format(filename, badpx_mask.shape, imageshape))
else:
logger('Info: badpixel mask loaded: {0}'.format(filename))
else:
logger('Warn: No badpixel mask found ({0}). Assuming no bad pixels'.format(filename))
badpx_mask = np.ones(imageshape)
return badpx_mask
"""def read_background(params, filename): # not used anymore as not really useful
"#""
Reads the background map and applies the corrections to it
:param params: Dictionary with all the parameters
:param filename: path and name of the background file
:return background: 2D numpy array of the background, normalised to 1s exposure time
"#""
imtype = 'background'
params['calibs'] = params['calibs_read']
if '{0}_calibs_read'.format(imtype) in params.keys():
params['calibs'] = params['{0}_calibs_read'.format(imtype)]
im, im_head = read_file_calibration(params, filename)
exptime = im_head[params['raw_data_exptim_keyword']]
background = im/exptime
return background"""
def warn_images_not_same(ims, names):
"""
Check that two images have the same size
:param ims: list or array of arrays
:param names: list of strings, same length as ims
"""
if len(ims) <= 1:
return
if len(ims) != len(names):
logger('Warn: len(ims) != len(names), that seems like a coding error.')
problems = ''
for i in range(len(ims)-1):
for j in range(i,len(ims)):
if ims[i].shape != ims[j].shape:
problems += '\tImage: {0} ({2}) and Image {1} ({3}){4}'.format(names[i], names[j], ims[i].shape, ims[j].shape, os.linesep)
if problems != '':
logger('Error: The following images have not the same size, but should have. '+\
'This is most likely caused by a missing "subframe" in one of the parameters "calib*". Please check.{1} {0}'.format(problems[:-1]), os.linesep )
return
def read_file_calibration(params, filename, level=0, realrun=True):
"""
Reads the filename and applies the calibrations as given in params['calibs']
This works also if the files are stored with the following header: BITPIX = 16, BZERO = 32768
:param params: Dictionary with all the parameters
:param filename: path and name of the file
:return: 2D numpy array of the file, and the header of the file read
"""
global calimages
im, im_head = read_fits_file(filename, realrun=realrun)
if realrun:
ims = im.shape
while len(ims) > 2:
if ims[0] == 1: # e.g. MRES from NARIT
im = im[0,:,:]
elif ims[-1] == 1:
im = im[:,:,0]
else:
logger(('Error: The file is stored in a multi-demensional array, which I do not know how to handle. The size of the image is {0}. '+\
'\tThis requires a small adjustment to the code in procedure read_file_calibration.').format(ims), params=params)
ims = im.shape
if params['raw_data_exptim_keyword'] in im_head.keys():
exptime = im_head[params['raw_data_exptim_keyword']] #saved as float
else:
exptime = 0.0
logger('Warn: Cannot find the raw_data_exptim_keyword = {0} in the header. Assuming 0 seconds.'.format(params['raw_data_exptim_keyword'] ))
if realrun:
logger('Info: image loaded: {0}'.format(filename))
im = rotate_flip_frame(im, params )
filename_s = filename[max(0,len(filename)-(80-25-1)):] # Shorten the filename so it fits into the header, 21 failed, 25 worked, between not tested
im_head['HIERARCH HiFLEx orig'] = filename_s
filename_nopath = filename_s.rsplit(os.sep,1)[-1]
im_head['HIERARCH HiFLEx orid'] = filename_nopath
calimages['{0}_saturated'.format(filename_s)] = np.where( im > params['max_good_value'] ) # Find the saturated pixel in the original image, consists of (x,y) with x and y are arrays
for entry in params['calibs']:
if entry == '':
continue
logtxt, headtxt = [], []
if entry.lower().find('subframe') > -1 and realrun:
subframe = copy.deepcopy(params[entry])
if len(subframe) >= 4:
# First: extend the image, if necessary
if subframe[2] < 0:
im = np.vstack(( np.zeros((-subframe[2], im.shape[1])), im ))
subframe[2] = 0
if subframe[3] < 0:
im = np.hstack(( np.zeros((im.shape[0], -subframe[3])), im ))
subframe[3] = 0
if im.shape[0] < subframe[0]+subframe[2]:
im = np.vstack(( im, np.zeros((subframe[0]+subframe[2]-im.shape[0], im.shape[1])) ))
if im.shape[1] < subframe[1]+subframe[3]:
im = np.hstack(( im, np.zeros((im.shape[0], subframe[1]+subframe[3]-im.shape[1])) ))
# Last: cut the image
if im.shape != (subframe[0],subframe[1]): # only apply subframe if the file doesn't have the size already
im = im[subframe[2]: subframe[0]+subframe[2], subframe[3]: subframe[1]+subframe[3]]
ims = im.shape
logger('Info: {1}: subframe applied: {0} ({2})'.format(entry, level, params[entry]))
calimages['{0}_saturated'.format(filename_s)] = np.where( im > params['max_good_value'] ) # Find the saturated pixel in the original image, do it again if the image was cut
im_head['HIERARCH HiFLEx redu{0}a'.format(level)] = 'Subframe: {0}'.format(entry)
elif entry.lower().find('bias') > -1:
if entry not in calimages:
create_image_general(params, entry, level=level+1)
if realrun:
warn_images_not_same([im, calimages[entry]], [filename,entry])
if np.percentile(calimages[entry], 90) > 2000 or np.percentile(calimages[entry], 10) < -100:
logger('Warn: The bias ({0}) has unphysical values: 10%-percentile = {1}, 90%-percentile = {2}'.format(entry, np.percentile(calimages[entry], 10), np.percentile(calimages[entry], 90)))
if np.isinf(calimages[entry]).any() or np.isnan(calimages[entry]).any():
logger('Warn: The bias ({0}) has {1} INF and {2} NAN values. Try to recreate the master bias.'.format(entry, np.sum(np.isnan(calimages[entry])), np.sum(np.isnan(calimages[entry])) ))
im = im - calimages[entry]
logtxt, headtxt = ['bias correction applied'], ['redu{0}b'.format(level), 'Bias']
elif entry.lower().find('dark') > -1:
if entry == 'dark': #only add exposure time, if just dark is given
entry = entry+str(exptime)
if entry not in calimages:
create_image_general(params, entry, level=level+1)
if realrun:
warn_images_not_same([im, calimages[entry]], [filename,entry])
if np.percentile(calimages[entry], 90) > 2000 or np.percentile(calimages[entry], 10) < -100:
logger('Warn: The dark ({0}) has unphysical values: 10%-percentile = {1}, 90%-percentile = {2}'.format(entry, np.percentile(calimages[entry], 10), np.percentile(calimages[entry], 90)))
if np.isinf(calimages[entry]).any() or np.isnan(calimages[entry]).any():
logger('Warn: The dark ({0}) has {1} INF and {2} NAN values. Try to recreate the master dark.'.format(entry, np.sum(np.isnan(calimages[entry])), np.sum(np.isnan(calimages[entry])) ))
im = im - calimages[entry]
logtxt, headtxt = ['dark correction applied'], ['redu{0}c'.format(level), 'Dark']
elif entry.lower().find('rflat') > -1:
if entry not in calimages:
create_image_general(params, entry, level=level+1)
if realrun:
warn_images_not_same([im, calimages[entry]], [filename,entry])
im = im / (calimages[entry] / np.median(calimages[entry]) )
logtxt, headtxt = ['flat correction applied with normalised flat (rflat)'], ['redu{0}d'.format(level), 'Flat']
elif entry.lower().find('badpx_mask') > -1 and realrun:
if entry not in calimages:
calimages[entry] = read_badpx_mask(params, ims)
warn_images_not_same([im, calimages[entry]], [filename,entry])
im = im * calimages[entry]
nonzeroind = np.nonzero(1-calimages[entry]) #find the indexes where the badpx mask is zero
for i in range(len(nonzeroind[0])): #find a suitable amound of surrounding pixels
for j in range(1,5): # try to find it in the surrounding 8, 24, 48, 80, 120 pixel
section = im[max(0,nonzeroind[0][i]-j) : min(ims[0],nonzeroind[0][i]+j+1) , max(0,nonzeroind[1][i]-j) : min(ims[1],nonzeroind[1][i]+j+1)]
section = section[section!=0] #only use areas <> 0 (exclude bad pixel
if len(section) >= 5:
break
if len(section) == 0:
logger('Warn: cannot replace bad pixel ({0}, {1}) with surrounding area in {2}'.format(nonzeroind[0][i],nonzeroind[1][i],filename))
else:
im[nonzeroind[0][i],nonzeroind[1][i]] = np.median(section) #replace bad px with the median of each surrounding area
logger('Info: {1}: badpx correction applied: {0}'.format(entry, level))
im_head['HIERARCH HiFLEx redu{0}f'.format(level)] = 'Bad-pixel-mask: {0}'.format(entry)
elif entry.lower().find('cosmic_ray') > -1 and realrun:
cr_setting = params['cosmic_ray_settings']
if cr_setting[0].lower() == 'deepcr' and deepCR is not None:
if multiprocessing.current_process().name == 'MainProcess': n_jobs = params['use_cores']
else: n_jobs = 1
mdl = deepCR(mask=cr_setting[1], inpaint=cr_setting[1], device="CPU") # mdl could be initialised earlier
#logger('Step: Cosmic ray removal')
mask, im_clean = mdl.clean(im, threshold=float(cr_setting[2]), segment=True, n_jobs=n_jobs) # best threshold is highest value that generate mask covering full extent of CR; choose threshold by visualizing outputs.
pos = np.where(mask)
#for ii in range(pos[0].shape[0]):
# print(filename, ii, pos[0][ii],pos[1][ii])
fname = filename.rsplit(os.sep,1)
if mask.any() and os.path.exists(params['path_reduced']) and params['path_reduced'].lower() != 'na'+os.sep: # at least one entry is True/1
im_head['Comment'] = 'In this file {0} cosmic rays were marked.'.format(np.sum(mask)) +\
' The first image shows the cleaned image, the second the mask and the third the uncleaned image'
save_im_fits(params, [im_clean, mask, im], im_head, params['path_reduced']+fname[-1].replace('.fit','_cr.fit'))
im = im_clean
logger('Info: {1}: cosmic ray correction applied: {0}'.format(entry, level))
im_head['HIERARCH HiFLEx redu{0}g'.format(level)] = 'Cosmic ray correction: {0}'.format(entry)
im_head['HIERARCH HiFLEx deepCR'] = '{0} pixel from cosmic rays'.format(np.sum(mask))
elif entry.lower().find('background') > -1 and realrun: # was localbackground before 20200108, this search covers *background
if 'sci_trace' in calimages.keys() and 'cal_trace' in calimages.keys():
logger('Step: Performing the background fit')
sci_tr_poly, xlows, xhighs, widths = calimages['sci_trace']
cal_tr_poly, axlows, axhighs, awidths = calimages['cal_trace']
bck_px_sci = find_bck_px(im, sci_tr_poly, xlows, xhighs, widths, params['background_width_multiplier'][0])
#save_im_fits(params, bck_px_sci, im_head, 'bck_px_sci')
bck_px_cal = find_bck_px(im, cal_tr_poly, axlows, axhighs, awidths, params['background_width_multiplier'][1])
#save_im_fits(params, bck_px_cal, im_head, 'bck_px_cal')
bck_px = bck_px_sci * bck_px_cal # mask with 0 for all the traces and 1 for background data
bad_values = ( im*bck_px > np.percentile(im[bck_px==1],95) ) # exclude brightest data in the background data (emission lines or borders of the traces
#save_im_fits(params, bck_px, im_head, 'bck_px1')
bck_px[bad_values] = 0
#save_im_fits(params, bck_px, im_head, 'bck_px2')
# Some deviation at the red side with not many lines, computational heavy, binning in dispersion direction made it faster
bck_im = fit_2d_image(im, params['polynom_bck'][1], params['polynom_bck'][0], w=bck_px, binning=[int(im.shape[0]/200.),max(1,int(im.shape[1]/2000.))]) # binning of 3 in cross-disp for image > 6000px, otherwise RAM is not enough (needs 20GB for 6000x9000 frame)
#plot_img_spec.plot_image(im, ['savepaths'], 1, True, [0.05,0.95,0.95,0.05], 'orig')
#plot_img_spec.plot_image(bck_px, ['savepaths'], 1, True, [0.05,0.95,0.95,0.05], 'bck_px')
#plot_img_spec.plot_image(bck_im, ['savepaths'], 1, True, [0.05,0.95,0.95,0.05], 'bck_im')
#plot_img_spec.plot_image(im-bck_im, ['savepaths'], 1, True, [0.05,0.95,0.95,0.05], 'diff')
plot_img_spec.plot_image((im - bck_im)*bck_px, \
[params['logging_path']+'background_subtracted-'+filename.rsplit(os.sep,1)[-1]+'.png'],\
1, False, [0.05,0.95,0.95,0.05], 'difference between image and background fit')
im = im - bck_im
bck_px[ bck_px==0 ] = np.nan
bck_noise_std, bck_noise_var = measure_background_noise(im * bck_px)
if np.isnan(bck_noise_var):
bck_noise_var = -1
logger('Info: {1}: background correction applied: {0}'.format(entry, level))
im_head['HIERARCH HiFLEx redu{0}f'.format(level)] = 'Background: {0}'.format(entry)
im_head['HIERARCH HiFLEx BCKNOISE'] = round(bck_noise_std,8)
im_head['HIERARCH HiFLEx BNOISVAR'] = (round(bck_noise_var,8), 'Variation of noise through image') # Background noise variation can be very high, because some light of the traces remains
else:
logger('Warn: Could not apply the calibration step {0} because the science and/or calibration traces are not yet known.'.format(entry))
elif entry.lower().find('combine_sum') > -1 or entry.lower().find('combine_mean') > -1 or entry.lower().find('normalise') > -1:
'nothing to do, as for a different step'
elif realrun:
logger('Warn: do not know what to do with this correction: {0}'.format(entry))
if len(logtxt) > 0 and len(headtxt) > 0 and realrun:
#print "np.where(np.isnan(calimages[entry])), calimages[entry].shape", np.where(np.isnan(calimages[entry])), calimages[entry].shape, np.median(calimages[entry]), np.nanmedian(calimages[entry]), np.mean(calimages[entry]), np.nanstd(calimages[entry], ddof=1, axis=None), np.where( np.isinf(calimages[entry]) )
#print(entry, np.sum(np.isnan(im)), np.sum(np.isnan(calimages[entry])), np.sum(np.isinf(im)), np.sum(np.isinf(calimages[entry])) )
im_median = np.nanmedian(calimages[entry], axis=None)
im_std = np.nanstd(calimages[entry], ddof=1, axis=None)
im_median = int(round(im_median))
if np.isnan(im_std): im_std = -1000
im_std = int(round(im_std))
logger('Info: {4}: {3}: {0} (median={1}, std={2})'.format(entry, im_median, im_std, logtxt[0], level))
im_head['HIERARCH HiFLEx '+headtxt[0]] = '{3}: {0}, median={1}, std={2}'.format(entry, im_median, im_std, headtxt[1])
#logger('Info: image loaded and processed: {0}'.format(filename))
if os.path.exists(params['path_reduced']) and params['path_reduced'].lower() != 'na'+os.sep and realrun: # Save the reduced image
fname = filename.rsplit(os.sep,1)
save_im_fits(params, im, im_head, params['path_reduced']+fname[-1])
return im, im_head
def create_image_general(params, imtype, level=0, realrun=True):
"""
Reads or creates the imtype file. If the key and file for master_<imtype>_filename exists the file is read, otherwise the file is created by combining the <imtype>_rawfiles
:param params: Dictionary with all the parameters
:param imtype: type of images, e.g. flat, dark5.0, bias
:param level: Sometimes the calibration images has to be created before it can be applied, to destinguish between this and the current image, the level can be increased by one
:return im: 2D array of the combined file
:return im_head: header of the last read file
"""
global calimages
mem = psutil.virtual_memory() # svmem(total=33221091328, available=28485840896, percent=14.3, used=4202041344, free=25513508864, active=..., inactive=..., buffers=..., cached=.., shared=...)
loaded = False
if 'master_{0}_filename'.format(imtype) in params.keys():
if params['master_{0}_filename'.format(imtype)] != '' and os.path.isfile(params['master_{0}_filename'.format(imtype)]) == True:
if realrun:
logger('Info: Using existing {0}: {1}'.format(imtype,params['master_{0}_filename'.format(imtype)]))
params['calibs'] = params['calibs_read']
if '{0}_calibs_read'.format(imtype) in params.keys():
params['calibs'] = params['{0}_calibs_read'.format(imtype)]
im, im_head = read_file_calibration(params, params['master_{0}_filename'.format(imtype)], level=level, realrun=realrun)
loaded = True
if loaded == False:
if '{0}_calibs_create'.format(imtype) not in params.keys():
if 'standard_calibs_create' not in params.keys() and realrun:
logger('Error: Missing entry in the configuration file. Neigther "{0}_calibs_create" nor "standard_calibs_create" is given. Please update the configuration file(s).'.format(imtype), params=params)
params['{0}_calibs_create'.format(imtype)] = params['standard_calibs_create']
for i in range(len(params['{0}_calibs_create'.format(imtype)])): # make it safe from different user input
params['{0}_calibs_create'.format(imtype)][i] = params['{0}_calibs_create'.format(imtype)][i].lower()
params['{0}_calibs_create'.format(imtype)][i] = params['{0}_calibs_create'.format(imtype)][i].replace('normaliz', 'normalis')
params['{0}_calibs_create'.format(imtype)][i] = params['{0}_calibs_create'.format(imtype)][i].replace('normalisation', 'normalise')
im, med_fluxes, std_fluxes = None, [], []
if '{0}_rawfiles'.format(imtype) not in params.keys() and realrun:
logger('Error: The list of raw files for image type {0} is not defined in the configuration. Please check the configuration files.'.format(imtype), params=params)
if len(params['{0}_rawfiles'.format(imtype)]) == 0 and realrun:
logger('Error: The list of raw files for image type {0} is empty. Please check the configuration files.'.format(imtype), params=params)
num_imgs = len(params['{0}_rawfiles'.format(imtype)]) # how many images are expected
header_updates = np.zeros((num_imgs,2))
for im_index, imf in enumerate(params['{0}_rawfiles'.format(imtype)]): # Only works for maximum 40 images on neils machine
params['calibs'] = params['{0}_calibs_create'.format(imtype)] # get's overwritten when other files are being read
img, im_head = read_file_calibration(params, imf, level=level, realrun=realrun)
if realrun:
im_head, obsdate_midexp, obsdate_mid_float, jd_midexp = get_obsdate(params, im_head) # unix_timestamp of mid exposure time
header_updates[im_index,:] = [im_head['HIERARCH HiFLEx EXPOSURE'], obsdate_mid_float] # !!! Improve this calculation and write in the header so it can be used later by get_obsdate
med_flux = np.median(img, axis=None)
med_fluxes.append(med_flux)
mask_inf = ~np.isinf(img)
std_fluxes.append(np.nanstd(img[mask_inf], axis=None, ddof=1))
if 'normalise' in params['{0}_calibs_create'.format(imtype)]:
img = img/(med_flux+0.0)
if im is None: # Initiate the array with correct precission to avoid swapping
if 2. * num_imgs * np.prod(img.shape) > mem[1] * 0.49: # had to convert to float as otherwise long_scalar overflow might happen
prec = np.float16
logger(('Warn: The ammount of pictures will most likely cause swapping, which might cause the system to become unresponsive.'+\
' Available memory: {0} GB, need to process {1} Giga-Pixel ({2}x{3}x{4}),'+\
' each pixel needs more than 1 byte to store the value.').format(round(mem[1]/1024.**3,2),
round(1.*num_imgs*np.prod(img.shape)/1024.**3,2), num_imgs, img.shape[0], img.shape[1]) )
elif 4. * num_imgs * np.prod(img.shape) > mem[1] * 0.49:
prec = np.float16
elif 8. * num_imgs * np.prod(img.shape) > mem[1] * 0.49:
prec = np.float32
else:
prec = np.float64
#logger('Test: Precision {0} is used'.format(prec))
im = np.zeros( (num_imgs, img.shape[0], img.shape[1]) ).astype(prec)
#print 'Created:', im.dtype, im.itemsize, im.nbytes, sys.getsizeof(im), im.nbytes/7979408000., mem
im[im_index,:,:] = img
#print im.dtype, im.itemsize, im.nbytes, sys.getsizeof(im), im.nbytes/7979408000.
if realrun:
for i in range(len(med_fluxes)):
im_head['HIERARCH HiFLEx NORM_{0}'.format(i)] = (med_fluxes[i], 'Median flux in image {0}'.format(i))
for i in range(len(std_fluxes)):
im_head['HIERARCH HiFLEx STDV_{0}'.format(i)] = (round(std_fluxes[i],5), 'Stdev of flux')
if 'combine_mean' in params['{0}_calibs_create'.format(imtype)]:
im = combine_sum(im)/(len(im)+0.0)
im_head['HIERARCH HiFLEx redu07'] = 'Average of {0} images'.format(len(params['{0}_rawfiles'.format(imtype)]))
exposure_time = np.mean(header_updates[:,0]) # Average of the exposure times
elif 'combine_sum' in params['{0}_calibs_create'.format(imtype)]:
im = combine_sum(im)
im_head['HIERARCH HiFLEx redu07'] = 'Sum of {0} images'.format(len(params['{0}_rawfiles'.format(imtype)]))
exposure_time = np.sum(header_updates[:,0]) # Sum of the exposure times
else: # Median combine
im = combine_median(im)
#imt = combine_median(im)
#aa = np.where(np.isinf(imt))
#for ii in range(aa[0].shape):
# print aa[0],aa[1], im[:,aa[0],aa[1]]
im_head['HIERARCH HiFLEx redu07'] = 'Median of {0} images'.format(len(params['{0}_rawfiles'.format(imtype)]))
exposure_time = np.median(header_updates[:,0]) # Median of the exposure times
if 'normalise' in params['{0}_calibs_create'.format(imtype)]:
norm_factor = np.median(med_fluxes)