GrisFitsFile.py 13.5 KB
Newer Older
Carl Schaffer's avatar
Carl Schaffer committed
1
2
3
4
5
6
"""
Created on 2018-04-09
@author Carl Schaffer
@mail   carl.schaffer@leibniz-kis.de
"""

Carl Schaffer's avatar
Carl Schaffer committed
7
import datetime  # date handling
8
9
import os  # path and file manipulations
import sys  # error handling
10
from glob import glob
11
from os.path import basename, exists, join, dirname
12
from warnings import warn
13

Carl Schaffer's avatar
linting    
Carl Schaffer committed
14
import astropy.io.fits as fitsio  # fits reading and writing
15
import astropy.units as u
16
import matplotlib.pyplot as plt
Carl Schaffer's avatar
Carl Schaffer committed
17
import numpy as np
18
from astropy.coordinates import SkyCoord
19
from astropy.wcs import WCS
Carl Schaffer's avatar
Carl Schaffer committed
20

21
22
23
from ..generic.fits import FitsFile
from ..util import calc_mu, calc_theta
from ..util.util import date_from_fn, gris_obs_mode
Carl Schaffer's avatar
Carl Schaffer committed
24

Carl Schaffer's avatar
Carl Schaffer committed
25

26
class GrisFitsFile(FitsFile):
Carl Schaffer's avatar
Carl Schaffer committed
27
28
29
30
31
32
33
34
    """Class for handling fits files for GRIS db
    When the constructor is called, the filename is stored and basic
    metadata parameters such as date and runnumber are extracted from
    the filename. Once the class function GrisFitsFile.parse() is
    called, metadata from the fits file is available as a dictionary
    stored at GrisFitsFile.properties. Any errors ocurring during
    opening and parsing of the file are stored as strings in
    GrisFitsFile.errors
Carl Schaffer's avatar
Carl Schaffer committed
35

36
37
38
39
40
    Args:
      filename: filename of target fits file

    Returns:

Carl Schaffer's avatar
Carl Schaffer committed
41
    """
Carl Schaffer's avatar
Carl Schaffer committed
42

Carl Schaffer's avatar
Carl Schaffer committed
43
44
    def __init__(self, filename):
        """Constructor"""
45
        super().__init__(filename)
Carl Schaffer's avatar
Carl Schaffer committed
46
        assert exists(filename)
Carl Schaffer's avatar
Carl Schaffer committed
47
48
49
50
        self.properties = {}
        self.errors = []

        self.filename = filename
51
        self.name = basename(filename)
52
        self.path = self.filename
53

54
        self.dirname = dirname(self.path)
Carl Schaffer's avatar
Carl Schaffer committed
55
56
57
58
        self.parse_filename()
        self.is_parsed = False

    def parse_filename(self):
59
        """Extract metadata from filename"""
Carl Schaffer's avatar
Carl Schaffer committed
60
        # Extract time from Filename
61
        name = os.path.splitext(basename(self.filename))[0].split("_")
Carl Schaffer's avatar
linting    
Carl Schaffer committed
62
63
64
        year = int(name[1][:4])
        month = int(name[1][4:6])
        day = int(name[1][6:])
Carl Schaffer's avatar
Carl Schaffer committed
65
66

        # Get time
Carl Schaffer's avatar
linting    
Carl Schaffer committed
67
68
69
        # hh = int(name[2][:2])
        # mm = int(name[2][2:4])
        # outstring = int(name[2][4:6])
Carl Schaffer's avatar
Carl Schaffer committed
70

Carl Schaffer's avatar
linting    
Carl Schaffer committed
71
        # self.date = datetime.datetime(year,month,day,hh,mm,outstring)
Carl Schaffer's avatar
Carl Schaffer committed
72
        self.date = datetime.datetime(year, month, day, 0, 0, 0)
Carl Schaffer's avatar
Carl Schaffer committed
73

74
        # Extract Run-, Map-, and Stepnumber, as well as mode string
75
        self.runnumber = int(name[4])
76
        self.mapnumber = int(name[5])
Carl Schaffer's avatar
linting    
Carl Schaffer committed
77
        self.slitposnumber = name[6]
78
79
        self.modestring = name[3]

80
    @property
Carl Schaffer's avatar
Carl Schaffer committed
81
    def verbose_runname(self):
82
83
84
85
86
87
        """
        Get verbose run name such as 24apr14.002

        Returns:
            outstring: verbose run name
        """
88
89
90
91
        outstring = self.date.strftime("%d%b%y").lower()
        outstring += f".{int(self.runnumber):03d}"
        return outstring

Carl Schaffer's avatar
Carl Schaffer committed
92
    @property
93
    def _coords_from_simple_header(self):
94
95
96
97
98
99
100
101
102
103
104
105
106
        """coordinates for slit, the Keywords 'SLITPOSX' and 'SLITPOSY' are assumed to be the center of the slit
         while the keyword 'SLITORIE' describes the rotation of the slit w.r.t. the helioprojective axes. The
         algorithm further assumes that each pixel represents a square area with side length of 'STEPSIZE'.
         These assumptions are used to calculate the coordinates for the center of each pixel within the data array.

        Returns:
            (X,Y):ndarray(float) x and y coordinates for each pixel of the slit


        WARNING! The code below is not checked for correctness! We do not know what direction SLITORIE references
        off and coordinates provided by GREGOR are very error-prone.
        """

Carl Schaffer's avatar
Carl Schaffer committed
107
        h = self.header
Carl Schaffer's avatar
Carl Schaffer committed
108
109
110
        slit_length_pixels = h["NAXIS2"]
        center_pixel = slit_length_pixels // 2

111
112
113
        x_slitcenter = h["SLITPOSX"] if h["SLITPOSX"] != "" else 0
        y_slitcenter = h["SLITPOSY"] if h["SLITPOSY"] != "" else 0
        angle = h["SLITORIE"] if h["SLITORIE"] != "" else 0
Carl Schaffer's avatar
Carl Schaffer committed
114
115
116

        delta_pix = [i - center_pixel for i in range(slit_length_pixels)]
        delta_pix = np.array(delta_pix)
117
118
        arcsec_per_pixel_x = h["STEPSIZE"]  # assuming square pixels
        arcsec_per_pixel_y = h["STEPSIZE"] if h["STEPSIZE"] else 0.135  # assuming square pixels
Carl Schaffer's avatar
Carl Schaffer committed
119

120
121
        X = x_slitcenter - np.sin(angle) * delta_pix * arcsec_per_pixel_x
        Y = y_slitcenter - np.cos(angle) * delta_pix * arcsec_per_pixel_y
Carl Schaffer's avatar
Carl Schaffer committed
122

123
124
125
126
127
128
129
130
131
        message = "Using FITS file coordinates, they are probably wrong by at least 100 arcsec!" \
                  " Also the slit angle is definitely off!"
        warn(message)

        coord_array = np.array([X, Y]) * u.arcsec
        # reorder axes for consistency with other coordinate retrieval methods
        coord_array = np.swapaxes(coord_array, 0, 1)

        return coord_array
132

133
134
135
136
137
138
139
140
    @property
    def coords(self):
        try:
            coords = self._coords_from_wcs
        except AssertionError:
            coords = self._coords_from_simple_header
        return coords

141
142
143
144
145
146
147
    @property
    def old_header(self):
        """Check if the header of this file has already been processed"""
        old_header = True
        if any([k in self.header for k in ['SOLARNET']]):
            old_header = False
        return old_header
Carl Schaffer's avatar
Carl Schaffer committed
148

149
150
151
152
153
154
155
    @property
    def mu(self):
        """
        Calculate mu from coordinates
        Returns: mu

        """
156
        coords = self.coords.value
157
        X, Y = coords[:, 0], coords[:, 1]
158
159
160
161
162
163
164
165
166
167
168
169
        x = np.mean(X)
        y = np.mean(Y)
        mu = calc_mu(x, y)
        return mu

    @property
    def theta(self):
        """
        Calculate theta from coordinates
        Returns: theta

        """
170
        coords = self.coords.value
171
        X, Y = coords[:, 0], coords[:, 1]
172
173
174
175
176
        x = np.mean(X)
        y = np.mean(Y)
        theta = calc_theta(x, y)
        return theta

Carl Schaffer's avatar
Carl Schaffer committed
177
178
    @property
    def slit_orientation(self):
179
180
181
182
183
        """Get slit_orientation, defaults to 0. Not sure what this angle means...

        Returns:
            slit_orientation: float
        """
184
185
186
187

        slit_orientation = self._get_value(["SLITORIE"])
        slit_orientation = slit_orientation if slit_orientation != "" else 0
        return slit_orientation
Carl Schaffer's avatar
Carl Schaffer committed
188
189

    @property
190
191
    def number_of_maps(self):
        """
192
193

        Returns:
194
            nmaps: int number of maps in associated observation
195
        """
Carl Schaffer's avatar
Carl Schaffer committed
196
        return self._get_value(["NMAPS", "SERIES"])
Carl Schaffer's avatar
Carl Schaffer committed
197
198
199

    @property
    def map_index(self):
200
201
202
203
204
        """

        Returns:
            i_map: int map index within observation
        """
Carl Schaffer's avatar
Carl Schaffer committed
205
        return self._get_value(["IMAP", "ISERIE"])
Carl Schaffer's avatar
Carl Schaffer committed
206
207
208

    @property
    def n_steps(self):
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
        """
        Determine the number of slit scanning steps in the associated map.

        Returns:
            nsteps: number of steps in the associated map

        """
        # get initial guess from header
        nsteps = self._get_value(["STEPS", "NSTEPS"])

        # check if other files matching the same run are present. If so, check the running step number
        # for the same map to determine the actual number of steps. It is assumed, that all files belonging to the
        # same map are stored in the same directory
        pattern = (
            f"gris_*_*_*_{self.runnumber:03d}_{self.mapnumber:03d}_*.fits"
        )
        matching_files = glob(join(self.dirname, pattern))
        max_steps = int(
Carl Schaffer's avatar
Carl Schaffer committed
227
228
229
230
231
232
            max(
                [
                    basename(m).split(".")[0].split("_")[-1]
                    for m in matching_files
                ]
            )
233
234
235
236
237
238
        )

        if max_steps < nsteps:
            nsteps = max_steps

        return nsteps
Carl Schaffer's avatar
Carl Schaffer committed
239

240
241
    @property
    def obs_time(self):
242
243
244
245
246
        """

        Returns:
            obs_time: observation time given in filename
        """
Carl Schaffer's avatar
Carl Schaffer committed
247
        return date_from_fn(self.filename)
248

Carl Schaffer's avatar
Carl Schaffer committed
249
250
    @property
    def step_index(self):
251
252
253
254
255
        """

        Returns:
            step_index: index of scanning step within current map
        """
Carl Schaffer's avatar
Carl Schaffer committed
256
        return self._get_value(["ISTEP"])
Carl Schaffer's avatar
Carl Schaffer committed
257

258
259
    @property
    def step_angle(self):
260
261
262
263
264
265
        """

        Returns:
            step_angle: scanning angle
        """

Carl Schaffer's avatar
Carl Schaffer committed
266
        return self._get_value(["STEPANGL"])
267
268
269
270

    @property
    def wavelength_step_size(self):
        """Get step length along wavelength axis. We use the values from the continuum correction fitting"""
Carl Schaffer's avatar
Carl Schaffer committed
271
        return self._kw_mean(["FF1WLDSP", "FF2WLDSP"])
272
273
274
275
276

    @property
    def wavelength_region(self):
        """Get wavelength region covered by file. Returns tuple of (upper, lower)"""
        step_size = self.wavelength_step_size
277
        n_steps_wl = self.data.shape[-1]
278
279

        # average for lower border:
Carl Schaffer's avatar
Carl Schaffer committed
280
281
        lower = self._kw_mean(["FF1WLOFF", "FF2WLOFF"])
        # add wl steps to calculate upper border
282
283
284
285
        upper = lower + step_size * n_steps_wl

        return lower, upper

286
287
288
289
290
291
292
293
294
    @property
    def derotator_in(self):
        """
        Determine whether the derotator is inserted

        Returns:
            bool: derotator in

        """
295
296
297
298
299
300
        try:
            rotval = self._get_value(["ROTCODE"])
        except KeyError:
            # if no derotator keywords are scontained, assume it doesn't exist yet
            return False

301
302
303
304
305
        if rotval == 0:
            return True
        else:
            return False

Carl Schaffer's avatar
Carl Schaffer committed
306
307
308
309
310
311
312
313
314
315
316
317
    @property
    def states(self):
        """
        Determine number of polarimetric states

        Returns:
            n_states int: number of states
        """

        # Number of polarimetric states is stored in the last axis. Determine
        # number of axes from NAXIS keyword and query value for last axis.
        naxes = self.query("NAXIS")
318
        keyword = f"NAXIS{naxes}"
Carl Schaffer's avatar
Carl Schaffer committed
319
320
321
322
        states = self.query(keyword)

        return states

323
324
    @property
    def obs_mode(self):
325
        """Get observation mode"""
326
        self.parse()
327
328
329
330
        # extract obs mode keywords
        states = self.states
        n_maps = self.number_of_maps
        n_steps = self.n_steps
331

332
        mode = gris_obs_mode(states, n_maps, n_steps)
333
        return mode
Carl Schaffer's avatar
Carl Schaffer committed
334

Carl Schaffer's avatar
Carl Schaffer committed
335
336
337
    def parse(self):
        """extract data and metadata from file"""
        if self.is_parsed:
Carl Schaffer's avatar
Carl Schaffer committed
338
            return
339
340
341
342
343
344
345
346
347

        def is_valid(value):
            valid = True
            if isinstance(value, fitsio.card.Undefined):
                valid = False
            elif value == "-NAN":
                valid = False

            return valid
Carl Schaffer's avatar
Carl Schaffer committed
348
349

        # Read file
Carl Schaffer's avatar
Carl Schaffer committed
350
        try:
Carl Schaffer's avatar
Carl Schaffer committed
351
            header = self.header
352

Carl Schaffer's avatar
Carl Schaffer committed
353
            # Add header to self.properties, cast strange data types where
354
355
            # necessary, fill invalid values with empty strings

Carl Schaffer's avatar
Carl Schaffer committed
356
            for key in header.keys():
357
358
359
                value = header[key]
                # check validity:
                value = value if is_valid(value) else ""
Carl Schaffer's avatar
Carl Schaffer committed
360
                if key in ["COMMENT", "HISTORY"]:
361
                    self.properties[key] = str(value)
Carl Schaffer's avatar
Carl Schaffer committed
362
                else:
363
364
                    self.properties[key] = value

Carl Schaffer's avatar
Carl Schaffer committed
365
            self.is_parsed = True
Carl Schaffer's avatar
Carl Schaffer committed
366
        except:
Carl Schaffer's avatar
linting    
Carl Schaffer committed
367
            exception = sys.exc_info()[0]
Carl Schaffer's avatar
Carl Schaffer committed
368
            self.errors.append("File %s" % self.filename + str(exception))
Carl Schaffer's avatar
Carl Schaffer committed
369
370

    def make_bson(self):
Carl Schaffer's avatar
Carl Schaffer committed
371
372
        """construct BSON document for mongodb storage"""
        return self.properties
Carl Schaffer's avatar
Carl Schaffer committed
373

374
    @property
375
    def _coords_from_wcs(self):
376
377
378
379
380
381
        """
        Retrieve spatial coordinates for each pixel in for slices projected into the spatial dimensions.
        Returns a data cube shaped s the spatial dimensions
        Returns: arcsec, cube of (NAXIS1, NAXIS2, 2) Values are Helioprojective X and Y in arcseconds
        """

382
383
        assert not self.old_header, f"Can't get coordinates from {self.path} ! Run header conversion first!"

384
385
386
        wcs = WCS(self.header)
        pixel_shape = wcs.pixel_shape
        n_y = pixel_shape[1]
387
        slit_list = np.zeros((n_y, self["NAXIS"]))
388
389
        slit_list[:, 1] = np.arange(n_y)
        slit_list[:, 0] = 0
390
        converted = np.array(wcs.pixel_to_world_values(slit_list))
Carl Schaffer's avatar
typo    
Carl Schaffer committed
391
        degrees = u.degree * converted[:, :2]
392
        arcsec = degrees.to(u.arcsec)
393

394
395
        return arcsec

396
    def plot_slit(self):
397
        m = self.full_disk_map
398

399
        # retrieve different data sources
400
        coordinates = self.coords
401
402
403
404

        X, Y = coordinates[:, 0], coordinates[:, 1]

        # determine fov limits for plotting
405

406
        plotting_padding = 300 * u.arcsec
407
408
409
410
        x_min = min(X) - plotting_padding
        x_max = max(X) + plotting_padding
        y_min = min(Y) - plotting_padding
        y_max = max(Y) + plotting_padding
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432

        # Create submap of FOV and plot
        submap = m.submap(SkyCoord(x_min, y_min, frame=m.coordinate_frame),
                          SkyCoord(x_max, y_max, frame=m.coordinate_frame))
        fig = plt.figure()  # noqa F841
        ax = plt.subplot(projection=submap)
        im = submap.plot()  # noqa:F841

        # create SkyCoord objects and overplot
        slit_wcs = SkyCoord(X, Y, frame=submap.coordinate_frame)
        ax.plot_coord(slit_wcs, 'g-', label="WCS Coordinates")

        # Create Legend and annotations
        title = "Gris Slit {0}_{1:03d}_{2:03d}_{3:03d}"
        title = title.format(self.obs_time.strftime('%Y%m%d'), self.runnumber, self.mapnumber, self.step_index)
        plt.title(title)

        annotation = f"Spectrum taken at \n{self.obs_time.strftime('%Y-%m-%d %H:%M:%S')}"
        annotation += f"\nHMI Reference taken at \n{submap.date.strftime('%Y-%m-%d %H:%M:%S')}"
        plt.text(0.05, 0.8, annotation, horizontalalignment='left', transform=ax.transAxes)
        plt.tight_layout(pad=3.5)
        _ = plt.legend()
433
434
435

        return fig, ax

Carl Schaffer's avatar
Carl Schaffer committed
436
    def __str__(self):
Carl Schaffer's avatar
linting    
Carl Schaffer committed
437
438
        outstring = "GrisFitsFile %s from %s:\n\tRun: %s\n\tMap: %s\n\tSlitposition: %s\n"
        outstring = outstring % (
Carl Schaffer's avatar
Carl Schaffer committed
439
            os.path.basename(self.filename),
Carl Schaffer's avatar
Carl Schaffer committed
440
441
            self.date,
            self.runnumber,
442
            self.map_index,
Carl Schaffer's avatar
Carl Schaffer committed
443
            self.slitposnumber,
Carl Schaffer's avatar
Carl Schaffer committed
444
        )
Carl Schaffer's avatar
linting    
Carl Schaffer committed
445
        return outstring
446
447
448
449
450
451


if __name__ == '__main__':
    gf = GrisFitsFile(sys.argv[1])
    gf.parse()
    print(gf)