-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathtdose_utilities.py
4692 lines (3968 loc) · 244 KB
/
tdose_utilities.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
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
import numpy as np
import os
import datetime
import sys
import astropy
from astropy import wcs
from astropy import units
from astropy import convolution
import astropy.convolution as ac # convolve, convolve_fft, Moffat2DKernel, Gaussian2DKernel
import astropy.io.fits as afits
from astropy.coordinates import SkyCoord
from astropy import units as u
from astropy.modeling.models import Sersic1D
from astropy.modeling.models import Sersic2D
from astropy.nddata import Cutout2D
import subprocess
import glob
import shutil
import scipy.ndimage
import scipy.special
import scipy.integrate as integrate
import tdose_utilities as tu
import tdose_model_FoV as tmf
from scipy.stats import multivariate_normal
import matplotlib as mpl
from matplotlib.colors import LogNorm
mpl.use('Agg') # prevent pyplot from opening window; enables closing ssh session with detached screen running TDOSE
import matplotlib.pylab as plt
import pdb
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def load_setup(setupfile='./tdose_setup_template.txt',verbose=True):
"""
Return dictionary with the setups found in 'setupfile'
(both TDOSE run and modification setup files can be loaded)
--- INPUT ---
setupfile The name of the txt file containing the TDOSE setup to load
Template for relevant setup files can be generated with
tdose_load_setup.generate_setup_template() or
tdose_load_setup.generate_setup_template_modify()
verbose Toggle verbosity
--- EXAMPLE OF USE ---
import tdose_utilities as tu
setup = tu.load_setup(setupfile='./tdose_setup_template.txt')
setup_modify = tu.load_setup(setupfile='./tdose_setup_template_modify.txt')
"""
if verbose: print(' --- tdose_utilities.load_setup() --- ')
#------------------------------------------------------------------------------------------------------
if verbose: print((' - Loading setup for TDOSE in '+setupfile))
setup_arr = np.genfromtxt(setupfile,dtype=None,names=None)
setup_dic = {}
for ii in np.arange(int(setup_arr.shape[0])):
paramname = setup_arr[ii,0].astype(str)
if paramname in list(setup_dic.keys()):
sys.exit(' Setup parameter "'+paramname+'" appears multiple times in the setup file\n '+
setupfile)
try:
val = float(setup_arr[ii,1].astype(str))
except:
val = setup_arr[ii,1].astype(str)
# - - - treatment of individual paramters - - -
if ('extension' in paramname) & (type(val) == float): val = int(val)
if (type(val) == str) or (type(val) == np.str_):
if val.lower() == 'none':
val = None
elif val.lower() == 'true':
val = True
elif val.lower() == 'false':
val = False
if (type(val) == str) or (type(val) == np.str_):
dirs = ['sources_to_extract','model_cube_layers','cutout_sizes']
if (paramname in dirs) & ('/' in str(val)):
val = val
setup_dic[paramname] = val
continue
lists = ['modify_sources_list','nondetections','model_cube_layers','sources_to_extract','plot_1Dspec_xrange','plot_1Dspec_yrange',
'plot_S2Nspec_xrange','plot_S2Nspec_yrange','cutout_sizes','aperture_size']
if (paramname in lists) & (val != 'all') & (val.lower() != 'none') & (val[0] == '['):
val = [float(vv) for vv in val.split('[')[-1].split(']')[0].split(',')]
setup_dic[paramname] = val
continue
if ('psf_sigma' in paramname):
if '/' in val:
sigmasplit = val.split('/')
if len(sigmasplit) != 2:
pass
else:
val = float(sigmasplit[0]) / float(sigmasplit[1])
setup_dic[paramname] = val
continue
setup_dic[paramname] = val
if verbose: print(' - Checking main keys are available; if not, adding them with None values')
checkkeys = ['nondetections','gauss_guess']
for ck in checkkeys:
if ck not in list(setup_dic.keys()):
setup_dic[ck] = None
if verbose: print(' - Returning dictionary containing setup parameters')
return setup_dic
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def generate_setup_template(outputfile='./tdose_setup_template.txt',clobber=False,verbose=True):
"""
Generate setup text file template
--- INPUT ---
outputfile The name of the output which will contain the TDOSE setup template
clobber Overwrite files if they exist
verbose Toggle verbosity
--- EXAMPLE OF USE ---
import tdose_utilities as tu
filename = './tdose_setup_template_new.txt'
tu.generate_setup_template(outputfile=filename,clobber=False)
setup = tu.load_setup(setupfile=filename)
"""
if verbose: print(' --- tdose_utilities.generate_setup_template() --- ')
#------------------------------------------------------------------------------------------------------
if os.path.isfile(outputfile) & (clobber == False):
sys.exit(' ---> Outputfile already exists and clobber=False ')
else:
if verbose: print((' - Will store setup template in '+outputfile))
if os.path.isfile(outputfile) & (clobber == True):
if verbose: print(' - Output already exists but clobber=True so overwriting it ')
setuptemplate = """
#-------------------------------------------------START OF TDOSE SETUP-------------------------------------------------
#
# Template for Three Dimensional Optimal Spectral Extracion (TDOSE, http://github.com/kasperschmidt/TDOSE) setup file
# Template was generated with tdose_utilities.generate_setup_template() on %s
# Setup file can be run with tdose.perform_extraction() or tdose.perform_extractions_in_parallel()
#
# - - - - - - - - - - - - - - - - - - - - - - - - - - DATA INPUT - - - - - - - - - - - - - - - - - - - - - - - - - - -
data_cube /path/datacube.fits # Path and name of fits file containing data cube to extract spectra from
cube_extension DATA_DCBGC # Name or number of fits extension containing data cube
variance_cube /path/variancecube.fits # Path and name of fits file containing variance cube to use for extraction
variance_extension VARCUBE # Name or number of fits extension containing noise cube
ref_image /path/referenceimage.fits # Path and name of fits file containing image to use as reference when creating source model
img_extension 0 # Name or number of fits extension containing reference image
wht_image /path/refimage_wht.fits # Path and name of fits file containing weight map of reference image (only cut out; useful for galfit modeling)
wht_extension 0 # Name or number of fits extension containing weight map
source_catalog /path/tdose_sourcecat.fits # Path and name of source catalog containing sources to extract spectra for
sourcecat_IDcol id # Column containing source IDs in source_catalog
sourcecat_xposcol x_image # Column containing x pixel position in source_catalog
sourcecat_yposcol y_image # Column containing y pixel position in source_catalog
sourcecat_racol ra # Column containing ra position in source_catalog (used to position cutouts if model_cutouts = True)
sourcecat_deccol dec # Column containing dec position in source_catalog (used to position cutouts if model_cutouts = True)
sourcecat_fluxcol fluxscale # Column containing a flux scale used for the modeling if no gauss_guess is provided
sourcecat_parentIDcol None # Column containing parent source IDs grouping source IDs into objects. Set to None to used id column
# corresponding to assigning each source to a single object
# if not None the parentid is used to group source models when storing 1D spectra. All models keep sources separate.
# - - - - - - - - - - - - - - - - - - - - - - - - OUTPUT DIRECTORIES - - - - - - - - - - - - - - - - - - - - - - - - -
models_directory /path/tdose_models/ # Directory to store the modeling output from TDOSE in
cutout_directory /path/tdose_cutouts/ # Directory to store image and cube cutouts in if model_cutouts=True
spec1D_directory /path/tdose_spectra/ # Output directory to store spectra in.
# - - - - - - - - - - - - - - - - - - - - - - - - - - CUTOUT SETUP - - - - - - - - - - - - - - - - - - - - - - - - - -
model_cutouts True # Perform modeling and spectral extraction on small cutouts of the cube and images to reduce run-time
cutout_sizes /path/tdose_setup_cutoutsizes.txt # Size of cutouts [ra,dec] in arcsec around each source to model.
# To use source-specific cutouts provide ascii file containing ID xsize[arcsec] and ysize[arcsec].
# - - - - - - - - - - - - - - - - - - - - - - - - SOURCE MODEL SETUP - - - - - - - - - - - - - - - - - - - - - - - - -
model_image_ext tdose_modelimage # Name extension of fits file containing reference image model. To ignored use None
model_param_reg tdose_modelimage_ds9 # Name extension of DS9 region file for reference image model. To ignored use None
model_image_cube_ext tdose_modelimage_cubeWCS # Name extension of fits file containing model image after conversion to cube WCS. To ignored use None.
source_model gauss # The source model to use for sources. Choices are:
# gauss Each source is modeled as a multivariate gaussian using the source_catalog input as starting point
# galfit The sources in the field-of-view are defined based on GALFIT header parameters; if all components are # Not enabled yet
# Gaussians an analytical convolution is performed. Otherwise numerical convolution is used. # Not enabled yet
# modelimg A model image exists, e.g., obtained with Galfit, in modelimg_directory. To disentangle/de-blend individual
# components, a model cube and parent_ids should be provided (see comments to modelimg_directory). If a model
# image is provded, TDOSE assumes it to represent the 1 object in the field-of-view.
# If the model image is not found a gaussian model of the FoV (source_model=gauss) is performed instead.
# aperture A simple aperture extraction on the datacubes is performed, i.e., no modeling of sources.
# - - - - - - - - - - - - - - - - - - - - - - - - GAUSS MODEL SETUP - - - - - - - - - - - - - - - - - - - - - - - - - -
gauss_guess /path/sextractoroutput.fits # To base initial guess of gaussian parameters on a SExtractor output provide SExtractor output fits file here
# If gauss_initguess=None the positions and flux scale provided in source_catalog will be used.
gauss_guess_idcol ID # Column of IDs in gauss_guess SExtractor file
gauss_guess_racol RA # Column of RAs in gauss_guess SExtractor file
gauss_guess_deccol DEC # Column of Decs in gauss_guess SExtractor file
gauss_guess_aimg A_IMAGE # Column of major axis in gauss_guess SExtractor file
gauss_guess_bimg B_IMAGE # Column of minor axis in gauss_guess SExtractor file
gauss_guess_angle THETA_IMAGE # Column of angle in gauss_guess SExtractor file
gauss_guess_fluxscale ACS_F814W_FLUX # Column of flux in gauss_guess SExtractor file to us for scaling
gauss_guess_fluxfactor 3 # Factor to apply to flux scale in initial Gauss parameter guess
gauss_guess_Nsigma 1 # Number of sigmas to include in initial Gauss parameter guess
max_centroid_shift 10 # The maximum shift of the centroid of each source allowed in the gaussian modeling. Given in pixels to
# set bounds ypix_centroid +/- max_centroid_shift and xpix_centroid +/- max_centroid_shift
# If none, no bounds are put on the centroid position of the sources.
# - - - - - - - - - - - - - - - - - - - - - - - - GALFIT MODEL SETUP - - - - - - - - - - - - - - - - - - - - - - - - -
galfit_directory /path/models_galfit/ # If source_model = galfit provide path to directory containing galfit models.
# TDOSE will look for galfit_*ref_image*_output.fits (incl. the cutout string if model_cutouts=True)
# If no model is found a source_model=gauss run on the object will be performed instead.
galfit_model_extension 2 # Fits extension containing galfit model with model parameters of each source in header.
# - - - - - - - - - - - - - - - - - - - - - - - - MODEL IMAGE SETUP - - - - - - - - - - - - - - - - - - - - - - - - -
modelimg_directory /path/models_cutouts/ # If source_model = modelimg provide the path to directory containing the individual source models
# TDOSE will look for model_*ref_image*.fits (incl. the cutout string if model_cutouts=True). If no model is found the object is skipped
# If a model image named model_*ref_image*_cube.fits is found, TDOSE assumes this file contains a cube with the individual model
# components isolated in individual layers of the cube. TDOSE will use this model instead of one generated within TDOSE.
# Parent IDs in the source catalog can be used to define what components belong to the object of interest (i.e., to extract a spectrum for)
# GALFIT models can be converted to TDOSE-suited model-cubes with tdose_utilities.galfit_convertmodel2cube()
# A TDOSE-suited model-cube can be build from individual 2D models with tdose_utilities.build_modelcube_from_modelimages()
modelimg_extension 0 # Fits extension containing model
# - - - - - - - - - - - - - - - - - - - - - - - - APERTURE MODEL SETUP - - - - - - - - - - - - - - - - - - - - - - - -
aperture_size 1.5 # Radius of apertures (float or list) to use given in arc seconds. For longer list of
# object-specific apertures provide ascii file containing ID and aperturesize[arcsec].
# - - - - - - - - - - - - - - - - - - - - - - - - - PSF MODEL SETUP - - - - - - - - - - - - - - - - - - - - - - - - - -
psf_type gauss # Select PSF model to build. Choices are:
# gauss Model the PSF as a symmetric Gaussian with sigma = FWHM/2.35482
# kernel_gauss An astropy.convolution.Gaussian2DKernel() used for numerical convolution # Not enabled yet
# kernel_moffat An astropy.convolution.Moffat2DKernel() used for numerical convolution # Not enabled yet
psf_FWHM_evolve linear # Evolution of the FWHM from blue to red end of data cube. Choices are:
# linear FWHM wavelength dependence described as FWHM(lambda) = p0[''] + p1[''/A] * (lambda - p2[A])
psf_FWHMp0 0.940 # p0 parameter to use when determining wavelength dependence of PSF
psf_FWHMp1 -3.182e-5 # p1 parameter to use when determining wavelength dependence of PSF
psf_FWHMp2 7050 # p2 parameter to use when determining wavelength dependence of PSF
psf_savecube True # To save fits file containing the PSF cube set psf_savecube = True
# This cube is used for the "source_model = modelimg" numerical PSF convolution
# - - - - - - - - - - - - - - - - - - - - - - - - - - - NON_DETECTIONS - - - - - - - - - - - - - - - - - - - - - - - -
nondetections None # List of IDs of sources in source_catalog that are not detected in the reference image or which
# have low flux levels in which cases the Gaussian modeling is likely to be inaccurate.
# For long list of objects provide ascii file containing ids.
# If source_model = gauss then sources will be extracted by replacing models within ignore_radius
# with a single point source in the reference image model, which will then
# be convolved with the PSF specified when extracting, as usual.
# If source_model = modelimg TDOSE assumes that the input model already represents the desired extraction model
# of the non-detection. I.e., if the object should be extracted as a (PSF
# convolved) point source, the model image should include a point source.
# Hence, for source_model = modelimg the keyword nondetections is ignored.
ignore_radius 0.3 # Models within a radius of ignore_radius [arcsec] of the non-detection location will be replaced with a
# point source for extractions with source_model = gauss before convolving with the PSF and adjusting the flux
# leves in each model cube layer.
# - - - - - - - - - - - - - - - - - - - - - - - - - CUBE MODEL SETUP - - - - - - - - - - - - - - - - - - - - - - - - -
model_cube_layers all # Layers of data cube to model [both end layers included]. If 'all' the full cube will be modeled.
# To model source-specific layers provide ascii file containing ID layerlow and layerhigh.
# If layerlow=all and layerhigh=all all layers will be modeled for particular source
model_cube_optimizer matrix # The optimizer to use when matching flux levels in cube layers:
# matrix Optimize fluxes analytically using matrix algebra to minimize chi squared of
# the equation set comparing model and data in each layer.
# nnls Optimize fluxes using Scipy's non-negative least squares solver restricting
# flux scales to >= 0 (assuming source models are non-negative too).
# curvefit Optimize fluxes numerically using least square fitting from scipy.optimize.curve_fit().
# Only enabled for analytic convolution of Gaussian source models.
# lstsq Optimize fluxes analytically using scipy.linalg.lstsq().
model_cube_ext tdose_modelcube # Name extension of fits file containing model data cube.
residual_cube_ext tdose_modelcube_residual # Name extension of fits file containing residual between model data cube and data. To ignored use None.
source_model_cube_ext tdose_source_modelcube # Name extension of fits file containing source model cube (used to modify data cube).
# - - - - - - - - - - - - - - - - - - - - - - - - SPECTRAL EXTRACTION - - - - - - - - - - - - - - - - - - - - - - - - -
sources_to_extract [8685,9262,10195,29743] # Sources in source_catalog to extract 1D spectra for.
# If sourcecat_parentIDcol is not None all associated spectra are included in stored object spectra
# If set to 'all', 1D spectra for all sources in source_catalog is produced (without grouping according to parents).
# For long list of objects provide ascii file containing containing ids (here parent grouping will be performed)
spec1D_name tdose_spectrum # Name extension to use for extracted 1D spectra
# - - - - - - - - - - - - - - - - - - - - - - - - - - - PLOTTING - - - - - - - - - - - - - - - - - - - - - - - - - - -
plot_generate True # Indicate whether to generate plots or not
plot_1Dspec_ext fluxplot # Name extension of pdf file containing plot of 1D spectrum
plot_1Dspec_xrange [4800,9300] # Range of x-axes (wavelength) for plot of 1D spectra
plot_1Dspec_yrange [-100,1500] # Range of y-axes (flux) for plot of 1D spectra
plot_1Dspec_shownoise True # Indicate whether to show the noise envelope in plot or not
plot_S2Nspec_ext S2Nplot # Name extension of pdf file containing plot of S/N spectrum
plot_S2Nspec_xrange [4800,9300] # Range of x-axes (wavelength) for plot of S2N spectra
plot_S2Nspec_yrange [-1,15] # Range of y-axes (S2N) for plot of S2N spectra
#--------------------------------------------------END OF TDOSE SETUP--------------------------------------------------
""" % (tu.get_now_string())
fout = open(outputfile,'w')
fout.write(setuptemplate)
fout.close()
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def generate_setup_template_modify(outputfile='./tdose_setup_template_modify.txt',clobber=False,verbose=True):
"""
Generate setup text file template for modifying data cubes
--- INPUT ---
outputfile The name of the output which will contain the TDOSE setup template
clobber Overwrite files if they exist
verbose Toggle verbosity
--- EXAMPLE OF USE ---
import tdose_utilities as tu
filename = './tdose_setup_template_modify_new.txt'
tu.generate_setup_template_modify(outputfile=filename,clobber=True)
setup = tu.load_setup(setupfile=filename)
"""
if verbose: print(' --- tdose_utilities.generate_setup_template_modify() --- ')
#------------------------------------------------------------------------------------------------------
if os.path.isfile(outputfile) & (clobber == False):
sys.exit(' ---> Outputfile already exists and clobber=False ')
else:
if verbose: print((' - Will store setup template in '+outputfile))
if os.path.isfile(outputfile) & (clobber == True):
if verbose: print(' - Output already exists but clobber=True so overwriting it ')
setuptemplate = """
#---------------------------------------------START OF TDOSE MODIFY SETUP---------------------------------------------
#
# Template for TDOSE (http://github.com/kasperschmidt/TDOSE) setup file for modifying data cubes.
# Generated with tdose_utilities.generate_setup_template_modify() on %s
# Cube modifications are performed with tdose_modify_cube.perform_modification(setupfile=setup_file_modify)
#
# - - - - - - - - - - - - - - - - - - - - - - - - - MODIFYING CUBE - - - - - - - - - - - - - - - - - - - - - - - - - -
data_cube /path/datacube.fits # Path and name of fits file containing data cube to modify
cube_extension DATA_DCBGC # Name or number of fits extension containing data cube
source_model_cube /path/tdose_source_modelcube.fits # Path and name of fits file containing source model cube
source_extension DATA_DCBGC # Name or number of fits extension containing source model cube
modified_cube_dir /path/to/output/ # Path of output directory to store modified cube in
modified_cube tdose_modified_datacube # Name extension of file containing modified data cube.
modify_sources_list [1,2,5] # List of IDs of sources to remove from data cube using source model cube.
# Corresponds to indices of source model cube so expects [0,Nmodelcomp-1]
# For long list of IDs provide path and name of file containing IDs (only)
sources_action remove # Indicate how to modify the data cube. Chose between:
# 'remove' Sources in modify_sources_list are removed from data cube
# 'keep' All sources except the sources in modify_sources_list are removed from data cube
#----------------------------------------------END OF TDOSE MODIFY SETUP----------------------------------------------
""" % (tu.get_now_string())
fout = open(outputfile,'w')
fout.write(setuptemplate)
fout.close()
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def duplicate_setup_template(outputdirectory,infofile,infohdr=2,infofmt="S250",
loopcols=['data_cube','cube_extension'],
namebase='MUSEWide_tdose_setup',clobber=False,verbose=True):
"""
Take a setup template generated with generate_setup_template() and duplicate it filling
it with information from a provided infofile, e.g., fill update PSF info, field names,
image names, source lists, etc.
--- INPUT ---
outputdirectory Directory to store setup templates in
infofile File containing info to replace values in template setup with
infohdr Number of header (comment) lines in infofile before the expected list of column names
infofmt Format of columns in infofile (format for all columns are needed; not just loopcols)
If just a single format string is provided, this will be used for all columns.
loopcols The name of the columns in the loopcols to perform replacements for. The columns should
correspond to keywords in the TDOSE setup file. The first column of the file should be
named 'setupname' and will be used to name the duplicated setup file (appending it to namebase).
if 'all', all columns in infofile will be attempted replaced.
namebase Name base to use for the setup templates
clobber Overwrite files if they exist
verbose Toggle verbosity
--- EXAMPLE OF USE ---
import tdose_utilities as tu
outputdir = '/Users/kschmidt/work/TDOSE/muse_tdose_setups/'
infofile = outputdir+'musewide_infofile.txt'
tu.duplicate_setup_template(outputdir,infofile,namebase='MUSEWide_tdose_setup',clobber=False,loopcols=['setupname','data_cube','cube_extension'])
"""
if verbose: print(' --- tdose_utilities.duplicate_setup_template_MUSEWide() --- ')
filename = outputdirectory+namebase+'.txt'
tu.generate_setup_template(outputfile=filename,clobber=clobber)
if ',' not in infofmt: #if a single common format is given count columns in infofile
copen = np.genfromtxt(infofile,skip_header=infohdr,names=True)
Ncol = len(copen[0])
infofmt = ','.join([infofmt]*Ncol)
copen = np.genfromtxt(infofile,skip_header=infohdr,names=True,dtype=infofmt)
if loopcols == 'all':
if verbose: print(' - loopcals="all" so will attempt replacement of all columns in infofile')
loopcols = np.asarray(copen.dtype.names).tolist()
Nfiles = len(copen[loopcols[0]])
if verbose: print((' - Performing replacements and generating the '+str(Nfiles)+' TDOSE setup templates ' \
'described in \n '+infofile))
for setupnumber in np.arange(int(Nfiles)):
replacements = copen[setupnumber]
newsetup = outputdirectory+namebase+'_'+replacements['setupname'].astype(str)+'.txt'
if os.path.isfile(newsetup) & (clobber == False):
if verbose: print(' - File '+newsetup+' already exists and clobber = False so moving on to next duplication ')
continue
else:
fout = open(newsetup,'w')
with open(filename,'r') as fsetup:
for setupline in fsetup:
if setupline.startswith('#'):
if "Generated with tdose_utilities.generate_setup_template()" in setupline:
nowstring = tu.get_now_string()
fout.write("# Generated with tdose_utilities.duplicate_setup_template() on "+nowstring+' \n')
else:
fout.write(setupline)
elif setupline == '\n':
fout.write(setupline)
else:
vals = setupline.split()
if vals[0] in loopcols:
replaceline = setupline.replace(' '+vals[1]+' ',' '+copen[vals[0]][setupnumber].astype(str)+' ')
else:
replaceline = setupline.replace(' '+vals[1]+' ',' NO_REPLACEMENT ')
newline = replaceline.split('#')[0]+'#'+\
'#'.join(setupline.split('#')[1:]) # don't include comment replacements
fout.write(newline)
fout.close()
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def build_2D_cov_matrix(sigmax,sigmay,angle,verbose=True):
"""
Build a covariance matrix for a 2D multivariate Gaussian
--- INPUT ---
sigmax Standard deviation of the x-compoent of the multivariate Gaussian
sigmay Standard deviation of the y-compoent of the multivariate Gaussian
angle Angle to rotate matrix by in degrees (clockwise) to populate covariance cross terms
verbose Toggle verbosity
--- EXAMPLE OF USE ---
import tdose_utilities as tu
covmatrix = tu.build_2D_cov_matrix(3,1,35)
"""
if verbose: print((' - Build 2D covariance matrix with varinaces (x,y)=('+str(sigmax)+','+str(sigmay)+\
') and then rotated '+str(angle)+' degrees'))
cov_orig = np.zeros([2,2])
cov_orig[0,0] = sigmay**2.0
cov_orig[1,1] = sigmax**2.0
angle_rad = (180.0-angle) * np.pi/180.0 # The (90-angle) makes sure the same convention as DS9 is used
c, s = np.cos(angle_rad), np.sin(angle_rad)
rotmatrix = np.matrix([[c, -s], [s, c]])
cov_rot = np.dot(np.dot(rotmatrix,cov_orig),np.transpose(rotmatrix)) # performing rot * cov * rot^T
return cov_rot
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def normalize_2D_cov_matrix(covmatrix,verbose=True):
"""
Calculate the normalization foctor for a multivariate gaussian from it's covariance matrix
However, not that gaussian returned by tu.gen_2Dgauss() is normalized for scale=1
--- INPUT ---
covmatrix covariance matrix to normaliz
verbose Toggle verbosity
"""
detcov = np.linalg.det(covmatrix)
normfac = 1.0 / (2.0 * np.pi * np.sqrt(detcov) )
return normfac
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def gen_noisy_cube(cube,type='poisson',gauss_std=0.5,verbose=True):
"""
Generate noisy cube based on input cube.
--- INPUT ---
cube Data cube to be smoothed
type Type of noise to generate
poisson Generates poissonian (integer) noise
gauss Generates gaussian noise for a gaussian with standard deviation gauss_std=0.5
gauss_std Standard deviation of noise if type='gauss'
verbose Toggle verbosity
--- EXAMPLE OF USE ---
import tdose_utilities as tu
datacube = np.ones(([3,3,3])); datacube[0,1,1]=5; datacube[1,1,1]=6; datacube[2,1,1]=8
cube_with_noise = tu.gen_noisy_cube(datacube,type='gauss',gauss_std='0.5')
"""
if verbose: print((' - Generating "'+str(type)+'" noise on data cube'))
if type == 'poisson':
cube_with_noise = np.random.poisson(lam=cube, size=None)
elif type == 'gauss':
cube_with_noise = cube + np.random.normal(loc=np.zeros(cube.shape),scale=gauss_std, size=None)
else:
sys.exit(' ---> type="'+type+'" is not valid in call to mock_cube_sources.generate_cube_noise() ')
return cube_with_noise
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def gen_psfed_cube(cube,type='gauss',type_param=[0.5,1.0],use_fftconvolution=False,verbose=True):
"""
Smooth cube with a 2D kernel provided by 'type', i.e., applying a model PSF smoothing to cube
--- INPUT ---
cube Data cube to be smoothed
type Type of smoothing kernel to apply
gauss Use 2D gaussian smoothing kernel
type_param expected: [stdev,(stdev_wave_scale)]
moffat Use a 2D moffat profile to represent the PSF
type_param expected: [gamma,alpha,(gamma_wave_scale,alpha_wave_scale)]
NB: If *wave_scale inputs are provided a list of scales to apply at each wavelength layer
(z-direction) of data cube is expected, hence, adding a wavelength dependence to the kernels.
type_param List of parameters for the smoothing kernel.
For expected paramters see notes in description of "type" keyword above.
use_fftconvolution Perform convolution in Foruire space with FFT
verbose Toggle verbosity
--- EXAMPLE OF USE ---
import tdose_utilities as tu
datacube = np.ones(([3,3,3])); datacube[0,1,1]=5; datacube[1,1,1]=6; datacube[2,1,1]=8
cube_smoothed = tu.gen_psfed_cube(datacube,type='gauss',type_param=[10.0,[1.1,1.3,1.5]])
--- EXAMPLE OF USE ---
"""
if verbose: print((' - Applying a '+type+' PSF to data cube'))
Nparam = len(type_param)
Nlayers = cube.shape[0]
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if type == 'gauss':
if Nparam == 1:
if verbose: print(' No wavelength dependence; duplicating kernel for all layers')
kernel = ac.Gaussian2DKernel(type_param[0])
kernels = [kernel]*Nlayers
elif Nparam == 2:
if verbose: print(' Wavelength dependence; looping over layers to generate kernels')
if Nlayers != len(type_param[1]):
sys.exit(' ---> The number of wavelength scalings provided ('+str(len(type_param[1]))+
') is different from the number of layers in cube ('+str(Nlayers)+')')
kernels = []
for ll in np.arange(int(Nlayers)):
kernel = ac.Gaussian2DKernel(type_param[0]*type_param[1][ll])
kernels.append(kernel)
else:
sys.exit(' ---> Invalid number of paramters provided ('+str(Nparam)+')')
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elif type == 'moffat':
if Nparam == 2:
if verbose: print(' No wavelength dependence; duplicating kernel for all layers')
kernel = ac.Moffat2DKernel(type_param[0],type_param[1])
kernels = [kernel]*Nlayers
elif Nparam == 4:
if verbose: print(' Wavelength dependence; looping over layers to generate kernels')
if (Nlayers != len(type_param[2])) or (Nlayers != len(type_param[3])):
sys.exit(' ---> The number of wavelength scalings provided ('+str(len(type_param[2]))+
' and '+str(len(type_param[3]))+
') are different from the number of layers in cube ('+str(Nlayers)+')')
kernels = []
for ll in np.arange(int(Nlayers)):
kernel = ac.Moffat2DKernel(type_param[0]*type_param[2][ll],type_param[1]*type_param[3][ll])
kernels.append(kernel)
else:
sys.exit(' ---> Invalid number of paramters provided ('+str(Nparam)+')')
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
else:
sys.exit(' ---> type="'+type+'" is not valid in call to mock_cube_sources.gen_smoothed_cube() ')
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if verbose: print((' - Applying convolution kernel ('+type+') to each wavelength layer '))
cube_psfed = tu.perform_2Dconvolution(cube,kernels,use_fftconvolution=use_fftconvolution,verbose=True)
return cube_psfed
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def perform_2Dconvolution(cube,kernels,use_fftconvolution=False,verbose=True):
"""
Perform 2D convolution in data cube layer by layer
--- INPUT ---
cube Data cube to convolve
kernels List of (astropy) kernels to apply on each (z/wavelengt)layer of the cube
use_fftconvolution To convolve in FFT space set this keyword to True
verbose Toggle verbosity
--- EXAMPLE OF USE ---
# see tdose_utilities.gen_psfed_cube()
"""
csh = cube.shape
cube_convolved = np.zeros(csh)
for zz in np.arange(int(csh[0])): # looping over wavelength layers of cube
layer = cube[zz,:,:]
if use_fftconvolution:
layer_convolved = ac.convolve_fft(layer, kernels[zz], boundary='fill')
else:
layer_convolved = ac.convolve(layer, kernels[zz], boundary='fill')
cube_convolved[zz,:,:] = layer_convolved
return cube_convolved
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def gen_aperture(imgsize,ypos,xpos,radius,pixval=1,showaperture=False,verbose=True):
"""
Generating an aperture image
--- INPUT ---
imgsize The dimensions of the array to return. Expects [y-size,x-size].
The aperture will be positioned in the center of a (+/-x-size/2., +/-y-size/2) sized array
ypos Pixel position in the y direction
xpos Pixel position in the x direction
radius Radius of aperture in pixels
showaperture Display image of generated aperture
verbose Toggle verbosity
--- EXAMPLE OF USE ---
import tdose_utilities as tu
apertureimg = tu.gen_aperture([20,40],10,5,10,showaperture=True)
apertureimg = tu.gen_aperture([2000,4000],900,1700,150,showaperture=True)
"""
if verbose: print(' - Generating aperture in image (2D array)')
y , x = np.ogrid[-ypos+1.:imgsize[0]-ypos+1., -xpos+1.:imgsize[1]-xpos+1.] # +1s make sure pixel indication starts at pixel 1,1
mask = x*x + y*y <= radius**2.
aperture = np.zeros(imgsize)
if verbose: print((' - Assigning pixel value '+str(pixval)+' to aperture'))
aperture[mask] = pixval
if showaperture:
if verbose: print(' - Displaying resulting image of aperture (added background noise)')
noisimg = np.random.normal(0,pixval/5.,imgsize)
noisimg[mask] = pixval
plt.imshow(noisimg,interpolation='none')
plt.grid()
plt.title('Generated aperture')
plt.show()
plt.ion()
return aperture
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def gen_2Dgauss(size,cov,scale,method='scipy',show2Dgauss=False,savefits=False,verbose=True):
"""
Generating a 2D gaussian with specified parameters
--- INPUT ---
size The dimensions of the array to return. Expects [ysize,xsize].
The 2D gauss will be positioned in the center of the array
cov Covariance matrix of gaussian, i.e., variances and rotation
Can be build with cov = build_2D_cov_matrix(stdx,stdy,angle)
scale Scaling the 2D gaussian. By default scale = 1 returns normalized 2D Gaussian.
I.e., np.trapz(np.trapz(gauss2D,axis=0),axis=0) = 1
method Method to use for generating 2D gaussian:
'scipy' Using the class multivariate_normal from the scipy.stats library
'matrix' Use direct matrix expression for PDF of 2D gaussian (slow!)
show2Dgauss Save plot of generated 2D gaussian
savefits Save generated profile to fits file
verbose Toggler verbosity
--- EXAMPLE OF USE ---
import tdose_utilities as tu
covmatrix = tu.build_2D_cov_matrix(4,1,5)
gauss2Dimg = tu.gen_2Dgauss([20,40],covmatrix,5,show2Dgauss=True)
gauss2Dimg = tu.gen_2Dgauss([9,9],covmatrix,1,show2Dgauss=True)
sigmax = 3.2
sigmay = 1.5
covmatrix = tu.build_2D_cov_matrix(sigmax,sigmay,0)
scale = 1 # returns normalized gaussian
Nsigwidth = 15
gauss2DimgNorm = tu.gen_2Dgauss([sigmay*Nsigwidth,sigmax*Nsigwidth],covmatrix,scale,show2Dgauss=True,savefits=True)
covmatrix = tu.build_2D_cov_matrix(4,2,45)
scale = 1 # returns normalized gaussian
gauss2DimgNorm = tu.gen_2Dgauss([33,33],covmatrix,scale,show2Dgauss=True,savefits=True)
"""
if verbose: print(' - Generating multivariate_normal object for generating 2D gauss using ')
if method == 'scipy':
if verbose: print(' scipy.stats.multivariate_normal.pdf() ')
mvn = multivariate_normal([0, 0], cov)
if verbose: print(' - Setting up grid to populate with 2D gauss PDF')
#x, y = np.mgrid[-np.ceil(size[0]/2.):np.floor(size[0]/2.):1.0, -np.ceil(size[1]/2.):np.floor(size[1]/2.):1.0] #LT170707
x, y = np.mgrid[-np.floor(size[0]/2.):np.ceil(size[0]/2.):1.0, -np.floor(size[1]/2.):np.ceil(size[1]/2.):1.0]
pos = np.zeros(x.shape + (2,))
pos[:, :, 0] = x; pos[:, :, 1] = y
gauss2D = mvn.pdf(pos)
elif method == 'matrix':
if verbose: print(' loop over matrix expression ')
gauss2D = np.zeros([np.int(np.ceil(size[0])),np.int(np.ceil(size[1]))])
mean = np.array([np.floor(size[0]/2.),np.floor(size[1]/2.)])
norm = 1/np.linalg.det(np.sqrt(cov))/2.0/np.pi
for xpix in np.arange(size[1]):
for ypix in np.arange(size[0]):
coordMmean = np.array([int(ypix),int(xpix)]) - mean
MTXexpr = np.dot(np.dot(np.transpose(coordMmean),np.linalg.inv(cov)),coordMmean)
gauss2D[int(ypix),int(xpix)] = norm * np.exp(-0.5 * MTXexpr)
if float(size[0]/2.) - float(int(size[0]/2.)) == 0.0:
ypos = np.asarray(size[0])/2.0-1.0
else:
ypos = np.floor(np.asarray(size[0])/2.0)
if float(size[1]/2.) - float(int(size[1]/2.)) == 0.0:
xpos = np.asarray(size[1])/2.0-1.0
else:
xpos = np.floor(np.asarray(size[1])/2.0)
gauss2D = tu.shift_2Dprofile(gauss2D,[ypos,xpos],showprofiles=False,origin=0)
if verbose: print((' - Scaling 2D gaussian by a factor '+str(scale)))
gauss2D = gauss2D*scale
if show2Dgauss:
savename = './Generated2Dgauss.pdf'
if verbose: print((' - Saving resulting image of 2D gaussian to '+savename))
plt.clf()
centerdot = gauss2D*0.0
center = [int(gauss2D.shape[0]/2.),int(gauss2D.shape[1]/2.)]
centerdot[center[1],center[0]] = 2.0*np.max(gauss2D)
print((' - Center of gaussian (pixelized - marked in plot):'+str(center)))
print((' - Center of gaussian (subpixel) :'+str([ypos,xpos])))
plt.imshow(gauss2D-centerdot,interpolation=None,origin='lower')
plt.colorbar()
plt.title('Generated 2D Gauss')
plt.savefig(savename)
plt.clf()
if savefits:
fitsname = './Generated2Dgauss.fits'
hduimg = afits.PrimaryHDU(gauss2D)
hdus = [hduimg]
hdulist = afits.HDUList(hdus) # turn header into to hdulist
hdulist.writeto(fitsname,overwrite=True) # write fits file
if verbose: print((' - Saved image of shifted profile to '+fitsname))
return gauss2D
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def gen_2Dsersic(size,parameters,normalize=False,show2Dsersic=False,savefits=False,verbose=True):
"""
Generating a 2D sersic with specified parameters using astropy's generator
--- INPUT ---
size The dimensions of the array to return. Expects [ysize,xsize].
The 2D gauss will be positioned in the center of the array
parameters List of the sersic parameters.
Expects [amplitude,effective radius, Sersic index,ellipticity,rotation angle]
The amplitude is the central surface brightness within the effective radius (Ftot/2 is within r_eff)
The rotation angle should be in degrees, counterclockwise from the positive x-axis.
normalize Normalize the profile so sum(profile img) = 1.
show2Dsersic Save plot of generated 2D Sersic
savefits Save generated profile to fits file
verbose Toggler verbosity
--- EXAMPLE OF USE ---
import tdose_utilities as tu
size = [30,40]
size = [31,41]
parameters = [1,6.7,1.7,1.0-0.67,17.76-90]
sersic2D = tu.gen_2Dsersic(size,parameters,show2Dsersic=True,savefits=True)
size = [30,30]
size = [31,31]
parameters = [1,5,1.7,0.5,45]
sersic2D = tu.gen_2Dsersic(size,parameters,show2Dsersic=True,savefits=True)
"""
x, y = np.meshgrid(np.arange(size[1]), np.arange(size[0]))
if float(size[0]/2.) - float(int(size[0]/2.)) == 0.0:
ypos = np.asarray(size[0])/2.0-0.5
else:
ypos = np.floor(np.asarray(size[0])/2.0)
if float(size[1]/2.) - float(int(size[1]/2.)) == 0.0:
xpos = np.asarray(size[1])/2.0-0.5
else:
xpos = np.floor(np.asarray(size[1])/2.0)
model = Sersic2D(amplitude=parameters[0], r_eff=parameters[1], n=parameters[2], ellip=parameters[3],
theta=parameters[4]*np.pi/180., x_0=xpos, y_0=ypos)
sersic2D = model(x, y)
if normalize:
sersic2D = sersic2D / np.sum(sersic2D)
if show2Dsersic:
plt.clf()
savename = './Generated2Dsersic.pdf'
if verbose: print((' - Displaying resulting image of 2D sersic in '+savename))
centerdot = sersic2D*0.0
center = [int(sersic2D.shape[0]/2.),int(sersic2D.shape[1]/2.)]
# centerdot[center[1],center[0]] = 2.0*np.max(sersic2D)
print((' - Center of Sersic (pixelized - marked in plot): '+str(center)))
plt.imshow(sersic2D,interpolation=None,origin='lower')
plt.colorbar()
plt.title('Generated 2D Sersic')
plt.savefig(savename)
plt.clf()
if savefits:
fitsname = './Generated2Dsersic.fits'
hduimg = afits.PrimaryHDU(sersic2D)
hdus = [hduimg]
hdulist = afits.HDUList(hdus) # turn header into to hdulist
hdulist.writeto(fitsname,overwrite=True) # write fits file
if verbose: print((' - Saved image of shifted profile to '+fitsname))
return sersic2D
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def get_2DsersicIeff(value,reff,sersicindex,axisratio,boxiness=0.0,returnFtot=False):
"""
Get the surface brightness value at the effective radius of a 2D sersic profile (given GALFIT Sersic parameters).
Ieff is calculated using ewuations (4) and (5) in Peng et al. (2010), AJ 139:2097.
This Ieff is what is referred to as 'amplitude' in astropy.modeling.models.Sersic2D
used in tdose_utilities.gen_2Dsersic()
--- INPUT ---
value If returnFtot=False "value" corresponds to Ftot of the profile (total flux for profile integrated
til r=infty) and Ieff will be returned.
If instead returnFtot=True "value" should provide Ieff so Ftot can be returned
reff Effective radius
sersicindex Sersic index of profile
axisratio Ratio between the minor and major axis (0<axisratio<1)
boxiness The boxiness of the profile
returnFtot If Ftot is not known, but Ieff is, set returnFtot=True to return Ftot instead (providing Ieff to "value")
--- EXAMPLE OF USE ---
Ieff = 1.0
reff = 25.0
sersicindex = 4.0
axisratio = 1.0
Ftot_calc = tu.get_2DsersicIeff(Ieff,reff,sersicindex,axisratio,returnFtot=True)
Ieff_calc = tu.get_2DsersicIeff(Ftot_calc,reff,sersicindex,axisratio)
size = 1000
x,y = np.meshgrid(np.arange(size), np.arange(size))
mod = Sersic2D(amplitude = Ieff, r_eff = reff, n=sersicindex, x_0=size/2.0, y_0=size/2.0, ellip=1-axisratio, theta=-1)
img = mod(x, y)
hducube = afits.PrimaryHDU(img)
hdus = [hducube]
hdulist = afits.HDUList(hdus)
hdulist.writeto('/Volumes/DATABCKUP2/TDOSEextractions/models_cutouts/model_sersic_spherical.fits',clobber=True)
"""
gam2n = scipy.special.gamma(2.0*sersicindex)
kappa = scipy.special.gammaincinv(2.0*sersicindex,0.5)
Rfct = np.pi * (boxiness + 2.) / (4. * scipy.special.beta(1./(boxiness+2.),1.+1./(boxiness+2.)) )
factor = 2.0 * np.pi * reff**2.0 * np.exp(kappa) * sersicindex * kappa**(-2*sersicindex) * gam2n * axisratio / Rfct
if returnFtot:
Ftot = value * factor
return Ftot
else:
Ieff = value / factor
return Ieff
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def shift_2Dprofile(profile,position,padvalue=0.0,showprofiles=False,origin=1,splineorder=3,savefits=False,verbose=True):
"""
Shift 2D profile to given position in array by rolling it in x and y.
Can move by sub-pixel amount using interpolation
--- INPUT ---
profile profile to shift
position position to move center of image (profile) to: [ypos,xpos]
padvalue the values to padd the images with when shifting profile
origin The orging of the position values. If 0-based pixels postions the
center calculation is updated to refelect this.
showprofiles Save plot of profile when shifted?
splineorder Order of spline interpolation to use when shifting
savefits Save a fitsfile of the shifted profile
verbose Toggle verbosity
--- EXAMPLE OF USE ---
profile = np.ones([35,35])
profile[17,17] = 5.0
fitsname = './Shifted2Dprofile_initial.fits'
hduimg = afits.PrimaryHDU(profile)
hdus = [hduimg]
hdulist = afits.HDUList(hdus)
hdulist.writeto(fitsname,clobber=True)
profile_shifted = tu.shift_2Dprofile(profile,[20.5,20.5],padvalue=0.0,showprofiles=False,origin=1,splineorder=3,savefits=True)
"""
profile_dim = profile.shape
yposition = np.asarray(position[0])
xposition = np.asarray(position[1])
if origin == 1:
yposition = yposition - 1.0
xposition = xposition - 1.0
ycenter_img = profile_dim[0]/2.-0.5 # sub-pixel center to use as reference when estimating shift
xcenter_img = profile_dim[1]/2.-0.5 # sub-pixel center to use as reference when estimating shift
yshift = np.float(yposition)-ycenter_img
xshift = np.float(xposition)-xcenter_img
profile_shifted = scipy.ndimage.interpolation.shift(profile, [yshift,xshift], output=None, order=splineorder,
mode='nearest', cval=0.0, prefilter=True)
if showprofiles:
plt.clf()
savename = './Shifted2Dprofile.pdf'
vmaxval = np.max(profile_shifted)
plt.imshow(profile_shifted,interpolation=None,origin='lower') # ,vmin=-vmaxval, vmax=vmaxval
plt.colorbar()
plt.title('Positioned Source')
plt.savefig(savename)
plt.clf()
if verbose: print((' - Saved image of shifted profile to '+savename))
if savefits:
fitsname = './Shifted2Dprofile.fits'
hduimg = afits.PrimaryHDU(profile_shifted)
hdus = [hduimg]
hdulist = afits.HDUList(hdus) # turn header into to hdulist
hdulist.writeto(fitsname,overwrite=True) # write fits file
if verbose: print((' - Saved image of shifted profile to '+fitsname))
return profile_shifted
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def roll_2Dprofile(profile,position,padvalue=0.0,showprofiles=False):
"""
Move 2D profile to given position in array by rolling it in x and y.
Note that the roll does not handle sub-pixel moves.
tu.shift_2Dprofile() does this using interpolation
--- INPUT ---
profile profile to shift
position position to move center of image (profile) to: [ypos,xpos]
padvalue the values to padd the images with when shifting profile
showprofiles Show profile when shifted?
--- EXAMPLE OF USE ---
tu.roll_2Dprofile(gauss2D,)
"""
profile_dim = profile.shape
yroll = np.int(np.round(position[0]-profile_dim[0]/2.))
xroll = np.int(np.round(position[1]-profile_dim[1]/2.))
profile_shifted = np.roll(np.roll(profile,yroll,axis=0),xroll,axis=1)
if showprofiles:
vmaxval = np.max(profile_shifted)
plt.imshow(profile_shifted,interpolation='none',vmin=-vmaxval, vmax=vmaxval)
plt.title('Positioned Source')
plt.show()
if yroll < 0:
profile_shifted[yroll:,:] = padvalue
else:
profile_shifted[:yroll,:] = padvalue
if xroll < 0:
profile_shifted[:,xroll:] = padvalue
else:
profile_shifted[:,:xroll] = padvalue
if showprofiles:
plt.imshow(profile_shifted,interpolation='none',vmin=-vmaxval, vmax=vmaxval)
plt.title('Positioned Source with 0s inserted')
plt.show()
return profile_shifted
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def get_now_string(withseconds=False):
"""
Retruning a string containing a formated version of the current data and time
--- INPUNT ---
withseconds To include seconds in the outputted string set this keyword to True
"""
if withseconds:
nowstr = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S.%f")
else:
nowstr = datetime.datetime.now().strftime("%Y-%m-%d %H:%M")
return nowstr
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def gen_gridcomponents(imgsize):
"""
Generate grid compoents, i.e. x and y indecese for a given image size
--- INPUT ---
imgsize size of image to generate grid points for (y,x)
"""
x = np.linspace(0, imgsize[1]-1, imgsize[1])