CIMP/atbdngproc.py at develop · mmiesch/CIMP · GitHub

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

"""
The purpose of this script is to process a bunch of images for use with NoiseGate to model and remove noise.  By construction, it uses the same `process` function as `atbdnoise.py` but applies it to many images.
"""
import os
import matplotlib.pyplot as plt
import numpy as np
import sunpy.visualization.colormaps as cm

from astropy.io import fits
from CIMP import Snapshot as snap
from CIMP import Enhance
from sunpy.net import attrs as a

#------------------------------------------------------------------------------

def process(snap, background = 'ratio', point = 'none', \
            detail = 'none', noise = 'none', equalize = False, \
            downsample = False, clip = None, rmin = 0.0, rmax = None, \
            params = None, rescale_output = False):

    if background == 'subtract':
        snap.mask_background(rmin = rmin, rmax = rmax, nonzero = False)
        snap.subtract_background(rescale=False)
    elif background == 'ratio':
        snap.mask_background(rmin = rmin, rmax = rmax, nonzero = True)
        snap.background_ratio(rescale=False)
    else:
        snap.mask_background(rmin = rmin, rmax = rmax, nonzero = False)

    snap.mask_annulus(rmin = rmin, rmax = rmax)

    if downsample:
        snap.downsample()

    snap.enhance(clip = clip, point = point, detail = detail, noise_filter = noise, \
                 equalize = equalize, params = params)

    # hit it with another mask after processing
    snap.mask_annulus(rmin=rmin, rmax = rmax, missingval = np.nanmin(snap.data))

    if rescale_output:
        snap.rescale()

    return

#------------------------------------------------------------------------------
# parameters common to all images

# L0.5 LASCO data
instrument = a.Instrument.lasco
detector = a.Detector.c3
dir='/home/mark.miesch/data/lasco_monthly/c3/2012_04/'
bgfile = dir+'background.fts'
rmin = 0.16
rmax = 1.0
params = None
rescale = False

# pick the directories to do.  Note that the background is computed for
# the month so it is probably most accurate in the middle of the month
dlist = ['14','15','16']

outdir = '/home/mark.miesch/data/lasco_monthly/c3/2012_04_processed/'

#------------------------------------------------------------------------------

fig = 2

if fig == 1:

    background = 'ratio'
    downsample = True
    clip       = (1.0,1.4)
    point      = 'omr'
    detail     = 'mgn'
    params     = (0.8,1.5)
    noise      = 'None'
    equalize   = False

elif fig == 2:

    # this is just background subtraction and normalization by exposure time.
    # This is intended to produce an L3 test data set for the CCOR ATBD

    outdir = '/home/mark.miesch/data/lasco_monthly/c3/lasco_c3_L0.5/'

    background = 'ratio'
    downsample = False
    clip       = None
    point      = 'none'
    detail     = 'none'
    noise      = 'None'
    equalize   = False

else:
    print("pick a valid figure number")
    exit()

#------------------------------------------------------------------------------
# loop over all directories for the month

idx = np.uint64(0)
for d in dlist:
    day = dir+d
    for file in os.listdir(day):
        fpath = day+'/'+file
        try:
            assert(os.path.isfile(fpath))
            assert("median" not in file)
            assert(fpath != bgfile)
            x = snap.snapshot(file = fpath, bgfile = bgfile, \
                instrument = instrument, detector = detector, \
                normalize = True)
            process(x, background = background, point = point, detail = detail, \
                    noise = noise, equalize = equalize, downsample = downsample, \
                    clip = clip, rmin = rmin, rmax = rmax, params = params, \
                    rescale_output = rescale)
            outfile = outdir+'image'+str(idx).zfill(3)+'.fts'
            print(outfile)
            header0 = x.header
            del header0['HISTORY']
            hdu_out = fits.PrimaryHDU(x.data,header0)
            hdu_out.writeto(outfile, overwrite = True)
            idx += np.uint64(1)
        except Exception as e:
            print(f"{e}\nSkipping file {file}")
            pass

        hdu_out.writeto(outfile, overwrite = True)


#------------------------------------------------------------------------------