-
-
Notifications
You must be signed in to change notification settings - Fork 85
/
core.py
1535 lines (1207 loc) · 49.1 KB
/
core.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
# Licensed under a 3-clause BSD style license - see LICENSE.rst
# This module implements the base CCDPROC functions
from __future__ import (absolute_import, division, print_function,
unicode_literals)
import numbers
import numpy as np
from astropy.extern import six
from astropy.units.quantity import Quantity
from astropy import units as u
from astropy.modeling import fitting
from astropy import stats
from astropy.nddata import StdDevUncertainty
from astropy.wcs.utils import proj_plane_pixel_area
from scipy import ndimage
from .ccddata import CCDData
from .utils.slices import slice_from_string
from .log_meta import log_to_metadata
__all__ = ['background_deviation_box', 'background_deviation_filter',
'ccd_process', 'cosmicray_median', 'cosmicray_lacosmic',
'create_deviation', 'flat_correct', 'gain_correct', 'rebin',
'sigma_func', 'subtract_bias', 'subtract_dark', 'subtract_overscan',
'transform_image', 'trim_image', 'wcs_project', 'Keyword']
# The dictionary below is used to translate actual function names to names
# that are FITS compliant, i.e. 8 characters or less.
_short_names = {
'background_deviation_box': 'bakdevbx',
'background_deviation_filter': 'bakdfilt',
'ccd_process': 'ccdproc',
'cosmicray_median': 'crmedian',
'create_deviation': 'creatvar',
'flat_correct': 'flatcor',
'gain_correct': 'gaincor',
'subtract_bias': 'subbias',
'subtract_dark': 'subdark',
'subtract_overscan': 'suboscan',
'trim_image': 'trimim',
'transform_image': 'tranim',
'wcs_project': 'wcsproj'
}
@log_to_metadata
def ccd_process(ccd, oscan=None, trim=None, error=False, master_bias=None,
dark_frame=None, master_flat=None, bad_pixel_mask=None,
gain=None, readnoise=None, oscan_median=True, oscan_model=None,
min_value=None, dark_exposure=None, data_exposure=None,
exposure_key=None, exposure_unit=None,
dark_scale=False):
"""Perform basic processing on ccd data.
The following steps can be included:
* overscan correction (:func:`subtract_overscan`)
* trimming of the image (:func:`trim_image`)
* create deviation frame (:func:`create_deviation`)
* gain correction (:func:`gain_correct`)
* add a mask to the data
* subtraction of master bias (:func:`subtract_bias`)
* subtraction of a dark frame (:func:`subtract_dark`)
* correction of flat field (:func:`flat_correct`)
The task returns a processed `~ccdproc.CCDData` object.
Parameters
----------
ccd : `~ccdproc.CCDData`
Frame to be reduced.
oscan : `~ccdproc.ccddata.CCDData`, str or None, optional
For no overscan correction, set to None. Otherwise proivde a region
of ccd from which the overscan is extracted, using the FITS
conventions for index order and index start, or a
slice from ccd that contains the overscan.
Default is ``None``.
trim : str or None, optional
For no trim correction, set to None. Otherwise proivde a region
of ccd from which the image should be trimmed, using the FITS
conventions for index order and index start.
Default is ``None``.
error : bool, optional
If True, create an uncertainty array for ccd.
Default is ``False``.
master_bias : `~ccdproc.CCDData` or None, optional
A master bias frame to be subtracted from ccd.
Default is ``None``.
dark_frame : `~ccdproc.CCDData` or None, optional
A dark frame to be subtracted from the ccd.
Default is ``None``.
master_flat : `~ccdproc.CCDData` or None, optional
A master flat frame to be divided into ccd.
Default is ``None``.
bad_pixel_mask : `numpy.ndarray` or None, optional
A bad pixel mask for the data. The bad pixel mask should be in given
such that bad pixels havea value of 1 and good pixels a value of 0.
Default is ``None``.
gain : `~astropy.units.Quantity` or None, optional
Gain value to multiple the image by to convert to electrons.
Default is ``None``.
readnoise : `~astropy.units.Quantity` or None, optional
Read noise for the observations. The read noise should be in
electrons.
Default is ``None``.
oscan_median : bool, optional
If true, takes the median of each line. Otherwise, uses the mean.
Default is ``True``.
oscan_model : `~astropy.modeling.Model` or None, optional
Model to fit to the data. If None, returns the values calculated
by the median or the mean.
Default is ``None``.
min_value : float or None, optional
Minimum value for flat field. The value can either be None and no
minimum value is applied to the flat or specified by a float which
will replace all values in the flat by the min_value.
Default is ``None``.
dark_exposure : `~astropy.units.Quantity` or None, optional
Exposure time of the dark image; if specified, must also provided
``data_exposure``.
Default is ``None``.
data_exposure : `~astropy.units.Quantity` or None, optional
Exposure time of the science image; if specified, must also provided
``dark_exposure``.
Default is ``None``.
exposure_key : `~ccdproc.Keyword`, str or None, optional
Name of key in image metadata that contains exposure time.
Default is ``None``.
exposure_unit : `~astropy.units.Unit` or None, optional
Unit of the exposure time if the value in the meta data does not
include a unit.
Default is ``None``.
dark_scale : bool, optional
If True, scale the dark frame by the exposure times.
Default is ``False``.
Returns
-------
occd : `~ccdproc.CCDData`
Reduded ccd.
Examples
--------
1. To overscan, trim and gain correct a data set::
>>> import numpy as np
>>> from astropy import units as u
>>> from ccdproc import CCDData
>>> from ccdproc import ccd_process
>>> ccd = CCDData(np.ones([100, 100]), unit=u.adu)
>>> nccd = ccd_process(ccd, oscan='[1:10,1:100]',
... trim='[10:100, 1:100]', error=False,
... gain=2.0*u.electron/u.adu)
"""
# make a copy of the object
nccd = ccd.copy()
# apply the overscan correction
if isinstance(oscan, CCDData):
nccd = subtract_overscan(nccd, overscan=oscan,
median=oscan_median,
model=oscan_model)
elif isinstance(oscan, six.string_types):
nccd = subtract_overscan(nccd, fits_section=oscan,
median=oscan_median,
model=oscan_model)
elif oscan is None:
pass
else:
raise TypeError('oscan is not None, a string, or CCDData object.')
# apply the trim correction
if isinstance(trim, six.string_types):
nccd = trim_image(nccd, fits_section=trim)
elif trim is None:
pass
else:
raise TypeError('trim is not None or a string.')
# create the error frame
if error and gain is not None and readnoise is not None:
nccd = create_deviation(nccd, gain=gain, readnoise=readnoise)
elif error and (gain is None or readnoise is None):
raise ValueError(
'gain and readnoise must be specified to create error frame.')
# apply the bad pixel mask
if isinstance(bad_pixel_mask, np.ndarray):
nccd.mask = bad_pixel_mask
elif bad_pixel_mask is None:
pass
else:
raise TypeError('bad_pixel_mask is not None or numpy.ndarray.')
# apply the gain correction
if isinstance(gain, Quantity):
nccd = gain_correct(nccd, gain)
elif gain is None:
pass
else:
raise TypeError('gain is not None or astropy.units.Quantity.')
# subtracting the master bias
if isinstance(master_bias, CCDData):
nccd = nccd.subtract(master_bias)
elif master_bias is None:
pass
else:
raise TypeError(
'master_bias is not None, numpy.ndarray, or a CCDData object.')
# subtract the dark frame
if isinstance(dark_frame, CCDData):
nccd = subtract_dark(nccd, dark_frame, dark_exposure=dark_exposure,
data_exposure=data_exposure,
exposure_time=exposure_key,
exposure_unit=exposure_unit,
scale=dark_scale)
elif dark_frame is None:
pass
else:
raise TypeError(
'dark_frame is not None or a CCDData object.')
# test dividing the master flat
if isinstance(master_flat, CCDData):
nccd = flat_correct(nccd, master_flat, min_value=min_value)
elif master_flat is None:
pass
else:
raise TypeError(
'master_flat is not None, numpy.ndarray, or a CCDData object.')
return nccd
@log_to_metadata
def create_deviation(ccd_data, gain=None, readnoise=None):
"""
Create a uncertainty frame. The function will update the uncertainty
plane which gives the standard deviation for the data. Gain is used in
this function only to scale the data in constructing the deviation; the
data is not scaled.
Parameters
----------
ccd_data : `~ccdproc.CCDData`
Data whose deviation will be calculated.
gain : `~astropy.units.Quantity` or None, optional
Gain of the CCD; necessary only if ``ccd_data`` and ``readnoise``
are not in the same units. In that case, the units of ``gain``
should be those that convert ``ccd_data.data`` to the same units as
``readnoise``.
Default is ``None``.
readnoise : `~astropy.units.Quantity` or None, optional
Read noise per pixel.
Default is ``None``.
{log}
Raises
------
UnitsError
Raised if ``readnoise`` units are not equal to product of ``gain`` and
``ccd_data`` units.
Returns
-------
ccd : `~ccdproc.CCDData`
CCDData object with uncertainty created; uncertainty is in the same
units as the data in the parameter ``ccd_data``.
"""
if gain is not None and not isinstance(gain, Quantity):
raise TypeError('gain must be a astropy.units.Quantity.')
if readnoise is None:
raise ValueError('must provide a readnoise.')
if not isinstance(readnoise, Quantity):
raise TypeError('readnoise must be a astropy.units.Quantity.')
if gain is None:
gain = 1.0 * u.dimensionless_unscaled
if gain.unit * ccd_data.unit != readnoise.unit:
raise u.UnitsError("units of data, gain and readnoise do not match.")
# Need to convert Quantity to plain number because NDData data is not
# a Quantity. All unit checking should happen prior to this point.
gain_value = float(gain / gain.unit)
readnoise_value = float(readnoise / readnoise.unit)
var = (gain_value * ccd_data.data + readnoise_value ** 2) ** 0.5
ccd = ccd_data.copy()
# ensure uncertainty and image data have same unit
var /= gain_value
ccd.uncertainty = StdDevUncertainty(var)
return ccd
@log_to_metadata
def subtract_overscan(ccd, overscan=None, overscan_axis=1, fits_section=None,
median=False, model=None):
"""
Subtract the overscan region from an image.
Parameters
----------
ccd : `~ccdproc.CCDData`
Data to have overscan frame corrected.
overscan : `~ccdproc.CCDData` or None, optional
Slice from ``ccd`` that contains the overscan. Must provide either
this argument or ``fits_section``, but not both.
Default is ``None``.
overscan_axis : 0, 1 or None, optional
Axis along which overscan should combined with mean or median. Axis
numbering follows the *python* convention for ordering, so 0 is the
first axis and 1 is the second axis.
If overscan_axis is explicitly set to None, the axis is set to
the shortest dimension of the overscan section (or 1 in case
of a square overscan).
Default is ``1``.
fits_section : str or None, optional
Region of ``ccd`` from which the overscan is extracted, using the FITS
conventions for index order and index start. See Notes and Examples
below. Must provide either this argument or ``overscan``, but not both.
Default is ``None``.
median : bool, optional
If true, takes the median of each line. Otherwise, uses the mean.
Default is ``False``.
model : `~astropy.modeling.Model` or None, optional
Model to fit to the data. If None, returns the values calculated
by the median or the mean.
Default is ``None``.
{log}
Raises
------
TypeError
A TypeError is raised if either ``ccd`` or ``overscan`` are not the
correct objects.
Returns
-------
ccd : `~ccdproc.CCDData`
CCDData object with overscan subtracted.
Notes
-----
The format of the ``fits_section`` string follow the rules for slices that
are consistent with the FITS standard (v3) and IRAF usage of keywords like
TRIMSEC and BIASSEC. Its indexes are one-based, instead of the
python-standard zero-based, and the first index is the one that increases
most rapidly as you move through the array in memory order, opposite the
python ordering.
The 'fits_section' argument is provided as a convenience for those who are
processing files that contain TRIMSEC and BIASSEC. The preferred, more
pythonic, way of specifying the overscan is to do it by indexing the data
array directly with the ``overscan`` argument.
Examples
--------
Creating a 100x100 array containing ones just for demonstration purposes::
>>> import numpy as np
>>> from astropy import units as u
>>> arr1 = CCDData(np.ones([100, 100]), unit=u.adu)
The statement below uses all rows of columns 90 through 99 as the
overscan::
>>> no_scan = subtract_overscan(arr1, overscan=arr1[:, 90:100])
>>> assert (no_scan.data == 0).all()
This statement does the same as the above, but with a FITS-style section::
>>> no_scan = subtract_overscan(arr1, fits_section='[91:100, :]')
>>> assert (no_scan.data == 0).all()
Spaces are stripped out of the ``fits_section`` string.
"""
if not (isinstance(ccd, CCDData) or isinstance(ccd, np.ndarray)):
raise TypeError('ccddata is not a CCDData or ndarray object.')
if ((overscan is not None and fits_section is not None) or
(overscan is None and fits_section is None)):
raise TypeError('specify either overscan or fits_section, but not '
'both.')
if (overscan is not None) and (not isinstance(overscan, CCDData)):
raise TypeError('overscan is not a CCDData object.')
if (fits_section is not None and
not isinstance(fits_section, six.string_types)):
raise TypeError('overscan is not a string.')
if fits_section is not None:
overscan = ccd[slice_from_string(fits_section, fits_convention=True)]
if overscan_axis is None:
overscan_axis = 0 if overscan.shape[1] > overscan.shape[0] else 1
if median:
oscan = np.median(overscan.data, axis=overscan_axis)
else:
oscan = np.mean(overscan.data, axis=overscan_axis)
if model is not None:
of = fitting.LinearLSQFitter()
yarr = np.arange(len(oscan))
oscan = of(model, yarr, oscan)
oscan = oscan(yarr)
if overscan_axis == 1:
oscan = np.reshape(oscan, (oscan.size, 1))
else:
oscan = np.reshape(oscan, (1, oscan.size))
else:
if overscan_axis == 1:
oscan = np.reshape(oscan, oscan.shape + (1,))
else:
oscan = np.reshape(oscan, (1,) + oscan.shape)
subtracted = ccd.copy()
# subtract the overscan
subtracted.data = ccd.data - oscan
return subtracted
@log_to_metadata
def trim_image(ccd, fits_section=None):
"""
Trim the image to the dimensions indicated.
Parameters
----------
ccd : `~ccdproc.CCDData`
CCD image to be trimmed, sliced if desired.
fits_section : str or None, optional
Region of ``ccd`` from which the overscan is extracted; see
`~ccdproc.subtract_overscan` for details.
Default is ``None``.
{log}
Returns
-------
trimmed_ccd : `~ccdproc.CCDData`
Trimmed image.
Examples
--------
Given an array that is 100x100,
>>> import numpy as np
>>> from astropy import units as u
>>> arr1 = CCDData(np.ones([100, 100]), unit=u.adu)
the syntax for trimming this to keep all of the first index but only the
first 90 rows of the second index is
>>> trimmed = trim_image(arr1[:, :90])
>>> trimmed.shape
(100, 90)
>>> trimmed.data[0, 0] = 2
>>> arr1.data[0, 0]
1.0
This both trims *and makes a copy* of the image.
Indexing the image directly does *not* do the same thing, quite:
>>> not_really_trimmed = arr1[:, :90]
>>> not_really_trimmed.data[0, 0] = 2
>>> arr1.data[0, 0]
2.0
In this case, ``not_really_trimmed`` is a view of the underlying array
``arr1``, not a copy.
"""
if (fits_section is not None and
not isinstance(fits_section, six.string_types)):
raise TypeError("fits_section must be a string.")
trimmed = ccd.copy()
if fits_section:
python_slice = slice_from_string(fits_section, fits_convention=True)
trimmed = trimmed[python_slice]
return trimmed
@log_to_metadata
def subtract_bias(ccd, master):
"""
Subtract master bias from image.
Parameters
----------
ccd : `~ccdproc.CCDData`
Image from which bias will be subtracted.
master : `~ccdproc.CCDData`
Master image to be subtracted from ``ccd``.
{log}
Returns
-------
result : `~ccdproc.CCDData`
CCDData object with bias subtracted.
"""
result = ccd.subtract(master)
result.meta = ccd.meta.copy()
return result
@log_to_metadata
def subtract_dark(ccd, master, dark_exposure=None, data_exposure=None,
exposure_time=None, exposure_unit=None,
scale=False):
"""
Subtract dark current from an image.
Parameters
----------
ccd : `~ccdproc.CCDData`
Image from which dark will be subtracted.
master : `~ccdproc.CCDData`
Dark image.
dark_exposure : `~astropy.units.Quantity` or None, optional
Exposure time of the dark image; if specified, must also provided
``data_exposure``.
Default is ``None``.
data_exposure : `~astropy.units.Quantity` or None, optional
Exposure time of the science image; if specified, must also provided
``dark_exposure``.
Default is ``None``.
exposure_time : str or `~ccdproc.Keyword` or None, optional
Name of key in image metadata that contains exposure time.
Default is ``None``.
exposure_unit : `~astropy.units.Unit` or None, optional
Unit of the exposure time if the value in the meta data does not
include a unit.
Default is ``None``.
scale: bool, optional
If True, scale the dark frame by the exposure times.
Default is ``False``.
{log}
Returns
-------
result : `~ccdproc.CCDData`
Dark-subtracted image.
"""
if not (isinstance(ccd, CCDData) and isinstance(master, CCDData)):
raise TypeError("ccd and master must both be CCDData objects.")
if (data_exposure is not None and
dark_exposure is not None and
exposure_time is not None):
raise TypeError("specify either exposure_time or "
"(dark_exposure and data_exposure), not both.")
if data_exposure is None and dark_exposure is None:
if exposure_time is None:
raise TypeError("must specify either exposure_time or both "
"dark_exposure and data_exposure.")
if isinstance(exposure_time, Keyword):
data_exposure = exposure_time.value_from(ccd.header)
dark_exposure = exposure_time.value_from(master.header)
else:
data_exposure = ccd.header[exposure_time]
dark_exposure = master.header[exposure_time]
if not (isinstance(dark_exposure, Quantity) and
isinstance(data_exposure, Quantity)):
if exposure_time:
try:
data_exposure *= exposure_unit
dark_exposure *= exposure_unit
except TypeError:
raise TypeError("must provide unit for exposure time.")
else:
raise TypeError("exposure times must be astropy.units.Quantity "
"objects.")
if scale:
master_scaled = master.copy()
# data_exposure and dark_exposure are both quantities,
# so we can just have subtract do the scaling
master_scaled = master_scaled.multiply(data_exposure / dark_exposure)
result = ccd.subtract(master_scaled)
else:
result = ccd.subtract(master)
result.meta = ccd.meta.copy()
return result
@log_to_metadata
def gain_correct(ccd, gain, gain_unit=None):
"""Correct the gain in the image.
Parameters
----------
ccd : `~ccdproc.CCDData`
Data to have gain corrected.
gain : `~astropy.units.Quantity` or `~ccdproc.Keyword`
gain value for the image expressed in electrons per adu.
gain_unit : `~astropy.units.Unit` or None, optional
Unit for the ``gain``; used only if ``gain`` itself does not provide
units.
Default is ``None``.
{log}
Returns
-------
result : `~ccdproc.CCDData`
CCDData object with gain corrected.
"""
if isinstance(gain, Keyword):
gain_value = gain.value_from(ccd.header)
elif isinstance(gain, numbers.Number) and gain_unit is not None:
gain_value = gain * u.Unit(gain_unit)
else:
gain_value = gain
result = ccd.multiply(gain_value)
return result
@log_to_metadata
def flat_correct(ccd, flat, min_value=None):
"""Correct the image for flat fielding.
The flat field image is normalized by its mean before flat correcting.
Parameters
----------
ccd : `~ccdproc.CCDData`
Data to be flatfield corrected.
flat : `~ccdproc.CCDData`
Flatfield to apply to the data.
min_value : float or None, optional
Minimum value for flat field. The value can either be None and no
minimum value is applied to the flat or specified by a float which
will replace all values in the flat by the min_value.
Default is ``None``.
{log}
Returns
-------
ccd : `~ccdproc.CCDData`
CCDData object with flat corrected.
"""
# Use the min_value to replace any values in the flat
use_flat = flat
if min_value is not None:
flat_min = flat.copy()
flat_min.data[flat_min.data < min_value] = min_value
use_flat = flat_min
# divide through the flat
flat_corrected = ccd.divide(use_flat)
# multiply by the mean of the flat
flat_corrected = flat_corrected.multiply(use_flat.data.mean() *
use_flat.unit)
flat_corrected.meta = ccd.meta.copy()
return flat_corrected
@log_to_metadata
def transform_image(ccd, transform_func, **kwargs):
"""Transform the image.
Using the function specified by transform_func, the transform will
be applied to data, uncertainty, and mask in ccd.
Parameters
----------
ccd : `~ccdproc.CCDData`
Data to be flatfield corrected.
transform_func : function
Function to be used to transform the data.
kwargs :
Additional keyword arguments to be used by the transform_func.
{log}
Returns
-------
ccd : `~ccdproc.CCDData`
A transformed CCDData object.
Notes
-----
At this time, transform will be applied to the uncertainy data but it
will only transform the data. This will not properly handle uncertainties
that arise due to correlation between the pixels.
These should only be geometric transformations of the images. Other
methods should be used if the units of ccd need to be changed.
Examples
--------
Given an array that is 100x100::
>>> import numpy as np
>>> from astropy import units as u
>>> arr1 = CCDData(np.ones([100, 100]), unit=u.adu)
The syntax for transforming the array using
`scipy.ndimage.shift`::
>>> from scipy.ndimage.interpolation import shift
>>> from ccdproc import transform_image
>>> transformed = transform_image(arr1, shift, shift=(5.5, 8.1))
"""
# check that it is a ccddata object
if not (isinstance(ccd, CCDData)):
raise TypeError('ccd is not a CCDData.')
# check that transform is a callable function
if not hasattr(transform_func, '__call__'):
raise TypeError('transform is not a function.')
# make a copy of the object
nccd = ccd.copy()
# transform the image plane
nccd.data = transform_func(nccd.data, **kwargs)
# transform the uncertainty plane if it exists
if nccd.uncertainty is not None:
nccd.uncertainty.array = transform_func(nccd.uncertainty.array,
**kwargs)
# transform the mask plane
if nccd.mask is not None:
mask = transform_func(nccd.mask, **kwargs)
nccd.mask = (mask > 0)
return nccd
@log_to_metadata
def wcs_project(ccd, target_wcs, target_shape=None, order='bilinear'):
"""
Given a CCDData image with WCS, project it onto a target WCS and
return the reprojected data as a new CCDData image.
Any flags, weight, or uncertainty are ignored in doing the
reprojection.
Parameters
----------
ccd : `~ccdproc.CCDData`
Data to be projected.
target_wcs : `~astropy.wcs.WCS` object
WCS onto which all images should be projected.
target_shape : two element list-like or None, optional
Shape of the output image. If omitted, defaults to the shape of the
input image.
Default is ``None``.
order : str, optional
Interpolation order for re-projection. Must be one of:
+ 'nearest-neighbor'
+ 'bilinear'
+ 'biquadratic'
+ 'bicubic'
Default is ``'bilinear'``.
{log}
Returns
-------
ccd : `~ccdproc.CCDData`
A transformed CCDData object.
"""
from reproject import reproject_interp
if not (ccd.wcs.is_celestial and target_wcs.is_celestial):
raise ValueError('one or both WCS is not celestial.')
if target_shape is None:
target_shape = ccd.shape
projected_image_raw, _ = reproject_interp((ccd.data, ccd.wcs),
target_wcs,
shape_out=target_shape,
order=order)
reprojected_mask = None
if ccd.mask is not None:
reprojected_mask, _ = reproject_interp((ccd.mask, ccd.wcs),
target_wcs,
shape_out=target_shape,
order=order)
# Make the mask 1 if the reprojected mask pixel value is non-zero.
# A small threshold is included to allow for some rounding in
# reproject_interp.
reprojected_mask = reprojected_mask > 1e-8
# The reprojection will contain nan for any pixels for which the source
# was outside the original image. Those should be masked also.
output_mask = np.isnan(projected_image_raw)
if reprojected_mask is not None:
output_mask = output_mask | reprojected_mask
# Need to scale counts by ratio of pixel areas
area_ratio = (proj_plane_pixel_area(target_wcs) /
proj_plane_pixel_area(ccd.wcs))
# If nothing ended up masked, don't create a mask.
if not output_mask.any():
output_mask = None
nccd = CCDData(area_ratio * projected_image_raw, wcs=target_wcs,
mask=output_mask,
header=ccd.header, unit=ccd.unit)
return nccd
def sigma_func(arr, axis=None):
"""
Robust method for calculating the deviation of an array. ``sigma_func``
uses the median absolute deviation to determine the standard deviation.
Parameters
----------
arr : `~ccdproc.CCDData` or `numpy.ndarray`
Array whose deviation is to be calculated.
axis : int, tuple of ints or None, optional
Axis or axes along which the function is performed.
If ``None`` it is performed over all the dimensions of
the input array. The axis argument can also be negative, in this case
it counts from the last to the first axis.
Default is ``None``.
Returns
-------
uncertainty : float
uncertainty of array estimated from median absolute deviation.
"""
return stats.median_absolute_deviation(arr) * 1.482602218505602
def setbox(x, y, mbox, xmax, ymax):
"""
Create a box of length mbox around a position x,y. If the box will
be out of [0,len] then reset the edges of the box to be within the
boundaries.
Parameters
----------
x : int
Central x-position of box.
y : int
Central y-position of box.
mbox : int
Width of box.
xmax : int
Maximum x value.
ymax : int
Maximum y value.
Returns
-------
x1 : int
Lower x corner of box.
x2 : int
Upper x corner of box.
y1 : int
Lower y corner of box.
y2 : int
Upper y corner of box.
"""
mbox = max(int(0.5 * mbox), 1)
y1 = max(0, y - mbox)
y2 = min(y + mbox + 1, ymax - 1)
x1 = max(0, x - mbox)
x2 = min(x + mbox + 1, xmax - 1)
return x1, x2, y1, y2
def background_deviation_box(data, bbox):
"""
Determine the background deviation with a box size of bbox. The algorithm
steps through the image and calculates the deviation within each box.
It returns an array with the pixels in each box filled with the deviation
value.
Parameters
----------
data : `numpy.ndarray` or `numpy.ma.MaskedArray`
Data to measure background deviation.
bbox : int
Box size for calculating background deviation.
Raises
------
ValueError
A value error is raised if bbox is less than 1.
Returns
-------
background : `numpy.ndarray` or `numpy.ma.MaskedArray`
An array with the measured background deviation in each pixel.
"""
# Check to make sure the background box is an appropriate size
# If it is too small, then insufficient statistics are generated
if bbox < 1:
raise ValueError('bbox must be greater than 1.')
# make the background image
barr = data * 0.0 + data.std()
ylen, xlen = data.shape
for i in range(int(0.5 * bbox), xlen, bbox):
for j in range(int(0.5 * bbox), ylen, bbox):
x1, x2, y1, y2 = setbox(i, j, bbox, xlen, ylen)
barr[y1:y2, x1:x2] = sigma_func(data[y1:y2, x1:x2])
return barr
def background_deviation_filter(data, bbox):
"""
Determine the background deviation for each pixel from a box with size of
bbox.
Parameters
----------
data : `numpy.ndarray`
Data to measure background deviation.