Repository navigation
Expand file tree
/
Copy pathcore.py
More file actions
3622 lines (3031 loc) · 129 KB
/
Copy pathcore.py
File metadata and controls
3622 lines (3031 loc) · 129 KB
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"""
import inspect
import logging
import math
import numbers
import operator
import sys
import warnings
from functools import partial
import array_api_compat
import array_api_extra as xpx
import numpy as np
from astropy import nddata, stats
from astropy import units as u
from astropy.modeling import fitting
from astropy.nddata import CCDData, NDData, StdDevUncertainty
from astropy.nddata.bitmask import (
bitfield_to_boolean_mask as _bitfield_to_boolean_mask,
)
from astropy.units.quantity import Quantity
from astropy.utils import deprecated, deprecated_renamed_argument
from astropy.utils.exceptions import AstropyUserWarning
from astropy.wcs.utils import proj_plane_pixel_area
from numpy import mgrid as np_mgrid
from numpy.ma import nomask as np_ma_nomask
from packaging import version as pkgversion
from scipy import ndimage
from . import _blocks
from ._ccddata_wrapper_for_array_api import (
_copy_ccddata,
_set_mask,
_slice_ccddata,
_unwrap_ccddata_for_array_api,
_wrap_ccddata_for_array_api,
)
from ._nanfuncs import _promote_to_real
from ._nanfuncs import median as _nanfuncs_median
from ._nanfuncs import nanmad as _nanfuncs_nanmad
from ._nanfuncs import nanmedian as _nanfuncs_nanmedian
from ._windowfilters import _window_any, _window_rank, _window_reduce
from .log_meta import log_to_metadata
from .utils.slices import slice_from_string
logger = logging.getLogger(__name__)
__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",
"median_filter",
"ccdmask",
"bitfield_to_boolean_mask",
"HostCopyWarning",
]
# 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",
}
def _is_array(arr):
"""
Check whether an object is an array by tring to find a namespace
for it.
Parameters
----------
arr : object
Object to be tested.
Returns
-------
is_array : bool
``True`` if arr is an array, ``False`` otherwise.
"""
try:
array_api_compat.array_namespace(arr)
except TypeError:
return False
return True
def _native_numpy(arr):
"""
Return a NumPy array with the same data in native byte order.
Parameters
----------
arr : `numpy.ndarray`
A NumPy array, possibly in non-native byte order.
Returns
-------
`numpy.ndarray`
``arr`` converted to its dtype's native scalar type. No copy is made
if one is not required.
Notes
-----
Data read from a FITS file is big-endian (e.g. dtype ``>f8``), a byte
order that only NumPy understands; other array namespaces reject such
an array outright. Requesting ``dtype=arr.dtype.type`` from
``xp.asarray`` works around that in NumPy and JAX but hands a NumPy
type object to the namespace, which ``array_api_strict`` warns about.
Converting to native byte order in NumPy first and handing the
resulting plain array to ``xp.asarray`` with no ``dtype`` avoids both
problems.
"""
return arr.astype(arr.dtype.type, copy=False)
def _to_numpy(arr):
"""
Return ``arr`` as a NumPy array, copying only when that is required.
Parameters
----------
arr : array-like
An array from any array-API namespace, on any device.
Returns
-------
`numpy.ndarray`
A NumPy array with the same shape, dtype and values as ``arr``.
Notes
-----
This is the deliberate host-side conversion for code that requires
NumPy, such as writing a FITS file through `astropy.io.fits`. Callers
must not hand the result back to the user in place of the input: the
input stays in its own array namespace.
``np.asarray`` only works for arrays on the namespace's default device,
so the array is first moved there, asking for that device through the
standard ``__array_namespace_info__().default_device()``. The move is a
no-op when the array is already there, and is skipped when the
namespace reports `None` as its default device (as JAX does), since
there is then nothing to move to.
"""
xp = array_api_compat.array_namespace(arr)
default_device = xp.__array_namespace_info__().default_device()
if default_device is not None:
arr = array_api_compat.to_device(arr, default_device)
return np.asarray(arr)
class HostCopyWarning(AstropyUserWarning):
"""
A non-NumPy array was copied to the host CPU to run a NumPy-only
operation. Results are returned in the namespace and on the device of
the primary input.
Silence it with::
import warnings
import ccdproc
warnings.filterwarnings("ignore", category=ccdproc.HostCopyWarning)
"""
def _from_numpy(arr, like=None, *, xp=None):
"""
Return ``arr`` in the array namespace and on the device of ``like``, or
on the default device of ``xp`` if there is no ``like``.
Parameters
----------
arr : `numpy.ndarray` or None
The NumPy array to convert. `None` is passed through unchanged, so
that an absent mask needs no special case at the call site.
like : array or None, optional
The array whose namespace and device the result should have. If
`None`, the result is on the default device of ``xp``, which must
then be given.
xp : array namespace, optional
Array namespace to convert into. If not provided, the namespace is
determined from ``like``.
Returns
-------
array or None
``arr`` as an array of ``xp`` on the device of ``like`` (on the
default device of ``xp`` if ``like`` is `None`), or `None` if
``arr`` is `None`. If ``arr`` and ``like`` both have a real
floating dtype, the result has the dtype of ``like``; otherwise,
including when there is no ``like``, it has the dtype of ``arr``,
unless ``xp`` cannot represent it: JAX without 64-bit mode
silently converts float64 to float32.
Notes
-----
This is the inverse of `_to_numpy`: it returns the result of a
NumPy-only operation to the caller's array namespace and device, so
that the caller never receives a NumPy array in place of what it
passed in.
The NumPy-only libraries decide the dtype of their own output:
``reproject`` and ``astropy.modeling`` always return float64 and
astroscrappy always returns float32. Floating input keeps its dtype
instead (#1023), so a floating result is cast to the floating dtype of
``like``. A boolean result, such as a mask, and the result for integer
input keep their own dtype: casting a floating result back to an
integer dtype would wrap negative values and could not hold the NaN
that ``wcs_project`` needs.
"""
if arr is None:
return None
xp = xp or array_api_compat.array_namespace(like)
device = None
dtype = None
if like is not None:
device = array_api_compat.device(like)
if np.isdtype(arr.dtype, "real floating") and xp.isdtype(
like.dtype, "real floating"
):
dtype = like.dtype
return xp.asarray(arr, dtype=dtype, device=device)
def _ccddata_from_numpy(ccd, like=None, *, xp=None):
"""
Convert the data and uncertainty of a NumPy-backed CCDData into ``xp``.
Parameters
----------
ccd : `~astropy.nddata.CCDData`
An image whose data and uncertainty are NumPy arrays, such as one
just read from a FITS file. It is changed in place.
like : array or None, optional
As for `_from_numpy`.
xp : array namespace, optional
As for `_from_numpy`.
Returns
-------
`~astropy.nddata.CCDData`
``ccd``, its data and uncertainty now arrays of ``xp`` on the device
of ``like`` (on the default device of ``xp`` if ``like`` is `None`),
each with the dtype it had.
Notes
-----
Arrays read from a FITS file may be big-endian, which only NumPy
accepts, so each is put in native byte order with `_native_numpy`
before it is converted.
The conversion does not go through `_from_numpy`, whose floating dtype
rule is for results computed from ``like`` (#1023). A file is not such
a result, so its arrays keep the dtype they were read with: a
``float64`` file combined with a ``float32`` image in memory must not
lose precision before ``combine`` casts it to its ``dtype``.
The mask is left in NumPy: astropy's `~astropy.nddata.CCDData` mask
setter converts any mask to NumPy, so a caller that needs the mask in
``xp`` moves it there itself.
"""
xp = xp or array_api_compat.array_namespace(like)
device = None if like is None else array_api_compat.device(like)
ccd.data = xp.asarray(_native_numpy(ccd.data), device=device)
if ccd.uncertainty is not None:
ccd.uncertainty.array = xp.asarray(
_native_numpy(ccd.uncertainty.array), device=device
)
return ccd
def _is_internal_frame(frame):
"""
Return `True` if ``frame`` belongs to ccdproc's own internal machinery.
A frame counts as internal when it is not in ``ccdproc.tests`` and its
module either is, or is inside, plain ``ccdproc``, or is astropy's
``deprecated_renamed_argument``/``deprecated`` decorator wrapper module
(``cosmicray_lacosmic`` is wrapped by it).
"""
module = frame.f_globals.get("__name__", "")
if module.startswith("ccdproc.tests"):
return False
return (
module == "ccdproc" or module.startswith("ccdproc.")
) or module == "astropy.utils.decorators"
def _caller_stacklevel():
"""
Compute a `warnings.warn` ``stacklevel`` that reaches the user's code.
Walk the call stack, starting one frame above the caller of this
function, until the first frame outside ccdproc's internal machinery
(see `_is_internal_frame`) is found. The returned level is counted
relative to the `warnings.warn` call itself, so it can be passed to it
directly as ``stacklevel``.
Notes
-----
A fixed ``stacklevel`` breaks as soon as a public function is called
through another layer of ccdproc, such as ``ccd_process`` calling
``subtract_overscan``, or wrapped by a decorator, such as
``cosmicray_lacosmic``'s ``@deprecated_renamed_argument``: the constant
that is right for a direct call points at the wrong line, or even the
wrong file, once there is an extra frame in between. Walking the stack
instead finds the real caller regardless of how many such frames sit
between it and the ``warnings.warn`` call. ``ccdproc.tests`` is
deliberately *not* treated as internal, since the test suite calls these
functions directly and expects the warning attributed to the test's own
call site, not to pytest or unittest machinery further up the stack.
"""
# Frame 0 is this function; frame 1 is its caller, the function that
# calls `warnings.warn`. ``stacklevel=1`` in `warnings.warn` means "the
# call to `warn` itself", so start counting from there.
frame = sys._getframe(1)
level = 1
while frame is not None:
if not _is_internal_frame(frame):
return level
frame = frame.f_back
level += 1
return level
def _warn_host_copy(function_name, *arrays):
"""
Warn once if any of ``arrays`` is not NumPy.
Parameters
----------
function_name : str
Name of the public function doing the copy, as it appears in the
warning message.
*arrays : array or None
Every array the caller is about to copy to the host. Entries that
are `None`, or that are not arrays at all (such as a bare Python
scalar), are skipped. A NumPy scalar counts as a NumPy array.
Nothing is warned about when all of them are NumPy, since no copy
happens then.
Notes
-----
Call this before converting any of ``arrays``. Exactly one warning is
issued per call of the public function, naming the namespace of the
first non-NumPy array, however many arrays end up crossing to the host.
Python's own once-per-location default filter then collapses repeated
calls from the same place in the user's code.
The decision is made from the arrays that actually cross, not from the
namespace of the main input, because those can differ: a NumPy ``ccd``
can come with a non-NumPy ``inbkg`` or ``invar``.
The ``stacklevel`` is computed by `_caller_stacklevel` rather than
passed in as a constant, because the right value depends on how the
public function was reached: directly, through another ccdproc
function such as ``ccd_process``, or through astropy's
``deprecated_renamed_argument`` decorator on ``cosmicray_lacosmic``.
Walking the stack finds the actual caller in every case, including
calls made from ``ccdproc.tests``, which are treated as external so
that the test suite's own call sites get the warning, not pytest.
"""
for arr in arrays:
if not _is_array(arr):
continue
xp = array_api_compat.array_namespace(arr)
if array_api_compat.is_numpy_namespace(xp):
continue
warnings.warn(
f"{function_name} runs on the host CPU, so {xp.__name__} array "
"data was copied to numpy.",
HostCopyWarning,
stacklevel=_caller_stacklevel(),
)
return
def _namespace_dtype(dtype, xp):
"""
Return the dtype object of an array namespace that corresponds to ``dtype``.
Parameters
----------
dtype : dtype-like
Anything NumPy accepts as a dtype (a builtin such as ``int``, a
string such as ``"float32"``, a NumPy scalar type or `numpy.dtype`),
or a dtype object of ``xp`` itself.
xp : array namespace
The namespace whose dtype is wanted.
Returns
-------
dtype
The dtype of ``xp`` with the same name as NumPy's reading of
``dtype``, or ``dtype`` itself if NumPy cannot interpret it (it is
then presumably already a dtype of ``xp``) or ``xp`` has no dtype of
that name.
Notes
-----
``xp.asarray`` and ``xp.astype`` require the namespace's own dtype
objects. A builtin type such as ``int`` is a valid NumPy dtype, but
array-api-strict rejects it, so NumPy is used here only to resolve the
name and the dtype is then looked up on ``xp``.
"""
try:
name = np.dtype(dtype).name
except TypeError:
return dtype
return getattr(xp, name, dtype)
def _namespace_from_module(module):
"""
Normalise a user-supplied array namespace or module.
Parameters
----------
module : array namespace or module
A raw module such as `numpy` or ``dask.array``, or an
array-api-compat namespace.
Returns
-------
array namespace
The array-api-compat namespace of ``module``.
Notes
-----
A raw module may lack array-API spellings this package relies on
(``xp.bool``, the ``device`` keyword), so it is round-tripped through
one of its own arrays to get the array-api-compat namespace. A namespace
that is already a compat one comes back unchanged.
"""
return array_api_compat.array_namespace(module.asarray(0))
def _library_name(xp):
"""Return the name of the array library behind the array namespace ``xp``."""
return xp.__name__.removeprefix("array_api_compat.")
def _namespace_of(**arrays):
"""
Return the array namespace shared by ``arrays``, checking that they agree.
Parameters
----------
**arrays : array, scalar or None
The array inputs of one call, keyed by the argument names to use in
an error message. `None` and scalars are skipped; see Notes. At
least one must be an array.
Returns
-------
array namespace
The array namespace of the arrays.
Raises
------
TypeError
If the arrays belong to different array namespaces, or if none of
``arrays`` is an array.
ValueError
If the arrays are on different devices.
Notes
-----
Depending on the pair of libraries, an array operation mixing two
namespaces or devices silently converts one array to the other's
library, possibly through NumPy on the host, keeps whichever library
happens to be on the left, or fails with an error that names neither
argument. Checking up front replaces all of those with one error that
names the arguments.
A scalar broadcasts against an array of any namespace, so it is not
checked. Scalars are Python numbers, NumPy scalars and 0-d NumPy arrays,
including a scalar `~astropy.units.Quantity`: each holds a single value,
such as ``np.std`` of a NumPy array returns.
Only the data arrays are meant to be passed in. astropy's
`~astropy.nddata.CCDData` converts any mask it is given to NumPy
whatever the library of its data, so checking masks would reject
ordinary input; uncertainties are not checked either.
"""
present = {
name: arr
for name, arr in arrays.items()
if _is_array(arr)
and not (isinstance(arr, np.generic | np.ndarray) and np.ndim(arr) == 0)
}
if not present:
raise TypeError(
f"none of {', '.join(arrays)} is an array; at least one must be."
)
names = list(present)
first = names[0]
xp = array_api_compat.array_namespace(present[first])
device = array_api_compat.device(present[first])
for name in names[1:]:
other_xp = array_api_compat.array_namespace(present[name])
if other_xp is not xp:
raise TypeError(
f"{name} comes from {_library_name(other_xp)} but {first} comes "
f"from {_library_name(xp)}; all array inputs must come from the "
"same array library."
)
other_device = array_api_compat.device(present[name])
if other_device != device:
raise ValueError(
f"{name} is on device {other_device} but {first} is on device "
f"{device}; all array inputs must be on the same device."
)
return xp
def _percentile_fallback(array, percentiles, xp=None):
"""
Try calculating percentile using namespace, otherwise fall back to
an implmentation that uses sort. As of the 2023 version of the array API
there is no percentile function in the API but there is a sort function.
Parameters
----------
array : array_like
Array from which to calculate the percentile.
percentiles : float or list-like
Percentile to calculate.
xp : array namespace, optional
Array namespace to use for calculations. If not provided, the
namespace will be determined from the array.
Returns
-------
percentile : float or list-like
Calculated percentile.
"""
xp = xp or array_api_compat.array_namespace(array)
try:
return xp.percentile(array, percentiles)
except AttributeError:
pass
percentile_array = xp.asarray(
percentiles,
dtype=xp.float64,
device=array_api_compat.device(array),
)
if bool(xp.any((percentile_array < 0) | (percentile_array > 100))):
raise ValueError("Percentiles must be in the range [0, 100]")
# Fall back to using sort
sorted_array = xp.sort(array)
indexes = xp.astype(
sorted_array.shape[0] * percentile_array / 100,
xp.int64,
)
indexes = xp.minimum(
indexes,
xp.asarray(
sorted_array.shape[0] - 1,
dtype=xp.int64,
device=array_api_compat.device(sorted_array),
),
)
return sorted_array[indexes]
def _median_fallback(array, axis, xp=None):
"""
Try calculating the median using the namespace, otherwise fall back to
`ccdproc._nanfuncs.median`. As of the 2023 version of the array API
there is no median function in the API.
Parameters
----------
array : array_like
Array from which to calculate the median.
axis : int
Axis along which to calculate the median.
xp : array namespace, optional
Array namespace to use for calculations. If not provided, the
namespace will be determined from the array.
Returns
-------
median : array
Median of ``array`` along ``axis``.
"""
xp = xp or array_api_compat.array_namespace(array)
try:
return xp.median(array, axis=axis)
except AttributeError:
# median is not part of the array API standard; fall back to an
# implementation built on nanmedian, which also matches
# numpy.median's NaN-propagating semantics.
return _nanfuncs_median(array, axis=axis, xp=xp)
# ---------------------------------------------------------------------------
# Window-filter dispatch
#
# scipy.ndimage's window filters are numpy-only, so numpy input keeps going
# to them -- unchanged, and with ndimage's own selection algorithm -- while
# every other namespace is served by ccdproc._windowfilters, which is written
# purely in terms of the array API. The four ``_dispatch_*`` below are the only
# places that choice is made, so that the call sites read the same whichever
# backend they run on.
#
# The two implementations agree exactly on finite input. They part company on
# NaN: ndimage sorts NaNs in with the rest of a window, the native rank
# filters exclude them. That difference is what issue #984 will use to drop
# masked pixels out of the median, at which point these wrappers grow a mask
# argument and numpy input with a mask starts taking the native path too.
# ---------------------------------------------------------------------------
def _namespace_or_numpy(data):
"""
Resolve the array namespace of ``data``, treating array-likes as numpy.
Parameters
----------
data : array or array-like
An array from any array API namespace, or a plain array-like such
as a nested list or tuple that `numpy.asarray` accepts.
Returns
-------
array namespace
``array_api_compat.array_namespace(data)`` for an array, `numpy`
for anything ``array_namespace`` rejects.
Raises
------
TypeError
If ``data`` defines ``__array_namespace__`` but its namespace
lookup fails anyway, so that a broken array is reported rather
than coerced through numpy.
Notes
-----
A bare ``array_namespace`` raises ``TypeError: list is not a supported
array type`` for the nested lists and tuples that `scipy.ndimage`
accepts and that every earlier ccdproc therefore accepted too, so those
fall back to numpy. The ``__array_namespace__`` gate cannot come
*first*, because array_api_compat recognises dask arrays by module and
they have no such attribute.
"""
try:
return array_api_compat.array_namespace(data)
except TypeError:
if hasattr(data, "__array_namespace__"):
raise
return np
def _dispatch(ndimage_func, native_func, data, *args, xp, mode="reflect"):
"""
Run one window filter: `scipy.ndimage` on numpy input, native otherwise.
Parameters
----------
ndimage_func : callable
The `scipy.ndimage` filter, called as
``ndimage_func(data, *extra, size=size, mode=mode)``.
native_func : callable
Its `ccdproc._windowfilters` counterpart, called as
``native_func(data, size, *extra, mode=mode, xp=xp)``.
data : array
Array to filter, belonging to ``xp``.
*args
``(*extra, size)``: the filter's own positional arguments, if any
(a percentile, a reduction callable), followed by the window size,
in ndimage's argument order.
xp : array namespace
Namespace of ``data``; a numpy namespace selects ``ndimage_func``.
mode : str, optional
Boundary mode, ``'reflect'`` by default as in ndimage.
Returns
-------
array
The filtered array, in ``data``'s namespace.
Notes
-----
ndimage calls a ``generic_filter`` callable once per window with the
values flattened, where `ccdproc._windowfilters._window_reduce` calls
it as ``func(stack, axis=-1)`` on a band at a time; a reduction taking
an ``axis`` keyword, as `sigma_func` does, suits both.
"""
*extra, size = args
if array_api_compat.is_numpy_namespace(xp):
return ndimage_func(data, *extra, size=size, mode=mode)
return native_func(data, size, *extra, mode=mode, xp=xp)
# Module globals, one per filter, so that a test can monkeypatch one.
_dispatch_median_filter = partial(
_dispatch, ndimage.median_filter, partial(_window_rank, percentile=50.0)
)
_dispatch_percentile_filter = partial(
_dispatch, ndimage.percentile_filter, _window_rank
)
_dispatch_maximum_filter = partial(_dispatch, ndimage.maximum_filter, _window_any)
_dispatch_generic_filter = partial(_dispatch, ndimage.generic_filter, _window_reduce)
def _mad_fallback(data, axis, ignore_nan, xp=None):
"""
Median absolute deviation written purely in terms of the array API.
This is the non-numpy branch of `sigma_func`;
`astropy.stats.median_absolute_deviation` is numpy-only, so any other
namespace computes the deviation here, on the device of ``data``.
Parameters
----------
data : array
Array whose median absolute deviation is to be calculated. Integer
and boolean input is promoted to the namespace's default real
floating dtype.
axis : int, tuple of int, list of int or None
Axis or axes along which the deviation is computed. ``None``
flattens ``data`` first; a tuple or list reduces over all the
listed axes. Negative values count from the last axis.
ignore_nan : bool
If `True`, NaNs are ignored; otherwise a NaN anywhere in a reduced
slice makes that slice's result NaN, as `numpy.median` does.
xp : array namespace, optional
Namespace to use. If not provided, it is determined from ``data``.
Returns
-------
mad : array
Median absolute deviation of ``data`` along ``axis``, with the
reduced axes removed (0-d when ``axis`` is `None`), in the namespace
and on the device of ``data``.
Raises
------
ValueError
If a tuple or list ``axis`` is out of bounds for ``data`` or lists
an axis more than once (including via a negative alias).
TypeError
If a tuple or list ``axis`` contains a bool or a non-integer entry.
Notes
-----
The namespace's own ``nanmedian``/``median`` is used when it has one;
otherwise the sort-based fallbacks from `ccdproc._nanfuncs` are used
instead, which cost O(n log n) along the reduced axes rather than the
O(n) of a native median, and yield NaN for an entirely NaN slice
without the ``RuntimeWarning`` numpy emits.
"""
xp = xp or array_api_compat.array_namespace(data)
def med(d, axis):
try:
return (xp.nanmedian if ignore_nan else xp.median)(d, axis=axis)
except AttributeError:
fallback = _nanfuncs_nanmedian if ignore_nan else _nanfuncs_median
return fallback(d, axis=axis, xp=xp)
# The deviation itself is _nanfuncs.nanmad, the same computation
# Combiner._nanmadstd uses; only the median tier passed in differs.
# nanmad's _setup owns the axis handling (``None`` flattens, a tuple
# or list is merged into a single trailing axis) and the promotion to
# the namespace's default real floating dtype, so ``med`` only ever
# sees a single integer axis.
return _nanfuncs_nanmad(data, axis=axis, xp=xp, median=med)
@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,
gain_corrected=True,
):
"""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 `~astropy.nddata.CCDData` object.
Parameters
----------
ccd : `~astropy.nddata.CCDData`
Frame to be reduced.
oscan : `~astropy.nddata.CCDData`, str or None, optional
For no overscan correction, set to None. Otherwise provide 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 provide 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 : `~astropy.nddata.CCDData` or None, optional
A master bias frame to be subtracted from ccd. The unit of the
master bias frame should match the unit of the image **after
gain correction** if ``gain_corrected`` is True.
Default is ``None``.
dark_frame : `~astropy.nddata.CCDData` or None, optional
A dark frame to be subtracted from the ccd. The unit of the
master dark frame should match the unit of the image **after
gain correction** if ``gain_corrected`` is True.
Default is ``None``.
master_flat : `~astropy.nddata.CCDData` or None, optional
A master flat frame to be divided into ccd. The unit of the
master flat frame should match the unit of the image **after
gain correction** if ``gain_corrected`` is True.
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 have a 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``.
gain_corrected : bool, optional
If True, the ``master_bias``, ``master_flat``, and ``dark_frame``
have already been gain corrected. Default is ``True``.
Raises
------
TypeError
If ``ccd`` is not a `~astropy.nddata.CCDData`, or if the data of the
images given as ``ccd``, ``oscan``, ``master_bias``, ``dark_frame``
and ``master_flat`` come from different array libraries.
ValueError
If the data of those images are on different devices.
Returns
-------
occd : `~astropy.nddata.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 astropy.nddata 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)
"""
if not isinstance(ccd, CCDData):
raise TypeError("ccd is not a CCDData object.")
# make a copy of the object
nccd = _copy_ccddata(ccd)
# Check the images here, rather than leaving it to the step functions,
# so that an error names the arguments of ccd_process.
_namespace_of(
**{
name: image.data if isinstance(image, CCDData) else None
for name, image in [
("ccd", nccd),
("oscan", oscan),
("master_bias", master_bias),
("dark_frame", dark_frame),
("master_flat", master_flat),
]
}
)
# apply the overscan correction
if isinstance(oscan, CCDData):
nccd = subtract_overscan(
nccd, overscan=oscan, median=oscan_median, model=oscan_model
)
elif isinstance(oscan, str):
nccd = subtract_overscan(