GrisFitsFile.py 13.4 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 math import radians
12
from os.path import basename, exists, join, dirname
13
from warnings import warn
14

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

23
24
from ..generic.fits import FitsFile
from ..util.util import date_from_fn, gris_obs_mode
Carl Schaffer's avatar
Carl Schaffer committed
25

Carl Schaffer's avatar
Carl Schaffer committed
26

27
class GrisFitsFile(FitsFile):
Carl Schaffer's avatar
Carl Schaffer committed
28
29
30
31
32
    """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
Carl Schaffer's avatar
typo    
Carl Schaffer committed
33
    stored at GrisFitsFile.properties. Any errors occurring during
Carl Schaffer's avatar
Carl Schaffer committed
34
35
    opening and parsing of the file are stored as strings in
    GrisFitsFile.errors
Carl Schaffer's avatar
Carl Schaffer committed
36

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

    Returns:

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

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

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

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

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

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

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

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

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

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

Carl Schaffer's avatar
Carl Schaffer committed
93
    @property
94
    def _coords_from_simple_header(self):
95
96
97
98
99
100
101
102
103
104
105
106
107
        """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
108
        h = self.header
Carl Schaffer's avatar
Carl Schaffer committed
109
110
111
        slit_length_pixels = h["NAXIS2"]
        center_pixel = slit_length_pixels // 2

112
113
        x_slitcenter = h["SLITPOSX"] if h["SLITPOSX"] != "" else 0
        y_slitcenter = h["SLITPOSY"] if h["SLITPOSY"] != "" else 0
114
        p0_angle = radians(h['P0ANGLE'])
115

116
117
        x_slitcenter, y_slitcenter = rotate_around_first_coord(np.array([x_slitcenter]), np.array([y_slitcenter]), p0_angle)
        print(x_slitcenter, y_slitcenter, h['SLITPOSX'],h['SLITPOSY'])
118
        slit_orientation_angle = h["SLITORIE"] if h["SLITORIE"] != "" else 0
Carl Schaffer's avatar
Carl Schaffer committed
119
120
121

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

127
128
        X = x_slitcenter - np.sin(slit_orientation_angle) * delta_pix * arcsec_per_pixel_x
        Y = y_slitcenter - np.cos(slit_orientation_angle) * delta_pix * arcsec_per_pixel_y
Carl Schaffer's avatar
Carl Schaffer committed
129

Carl Schaffer's avatar
Carl Schaffer committed
130
131
        message = (
            "Using FITS file coordinates, they are probably wrong by at least 100 arcsec!"
132
            " Also the slit slit_orientation_angle is definitely off!"
Carl Schaffer's avatar
Carl Schaffer committed
133
        )
134
135
136
137
138
139
140
        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
141

142
143
144
145
146
147
148
149
    @property
    def coords(self):
        try:
            coords = self._coords_from_wcs
        except AssertionError:
            coords = self._coords_from_simple_header
        return coords

150
151
152
153
    @property
    def old_header(self):
        """Check if the header of this file has already been processed"""
        old_header = True
Carl Schaffer's avatar
Carl Schaffer committed
154
        if any([k in self.header for k in ["SOLARNET"]]):
155
156
            old_header = False
        return old_header
Carl Schaffer's avatar
Carl Schaffer committed
157

Carl Schaffer's avatar
Carl Schaffer committed
158
159
    @property
    def slit_orientation(self):
160
161
162
163
164
        """Get slit_orientation, defaults to 0. Not sure what this angle means...

        Returns:
            slit_orientation: float
        """
165
166
167
168

        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
169
170

    @property
171
172
    def number_of_maps(self):
        """
173
174

        Returns:
175
            nmaps: int number of maps in associated observation
176
        """
Carl Schaffer's avatar
Carl Schaffer committed
177
        return self._get_value(["NMAPS", "SERIES"])
Carl Schaffer's avatar
Carl Schaffer committed
178
179
180

    @property
    def map_index(self):
181
182
183
184
185
        """

        Returns:
            i_map: int map index within observation
        """
Carl Schaffer's avatar
Carl Schaffer committed
186
        return self._get_value(["IMAP", "ISERIE"])
Carl Schaffer's avatar
Carl Schaffer committed
187
188
189

    @property
    def n_steps(self):
190
191
192
193
194
195
196
197
198
199
200
201
202
        """
        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
Carl Schaffer's avatar
Carl Schaffer committed
203
        pattern = f"gris_*_*_*_{self.runnumber:03d}_{self.mapnumber:03d}_*.fits"
204
205
        matching_files = glob(join(self.dirname, pattern))
        max_steps = int(
Carl Schaffer's avatar
Carl Schaffer committed
206
            max([basename(m).split(".")[0].split("_")[-1] for m in matching_files])
207
208
209
210
211
212
        )

        if max_steps < nsteps:
            nsteps = max_steps

        return nsteps
Carl Schaffer's avatar
Carl Schaffer committed
213

214
215
    @property
    def obs_time(self):
216
217
218
219
220
        """

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

Carl Schaffer's avatar
Carl Schaffer committed
223
224
    @property
    def step_index(self):
225
226
227
228
229
        """

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

232
233
    @property
    def step_angle(self):
234
235
236
237
238
239
        """

        Returns:
            step_angle: scanning angle
        """

Carl Schaffer's avatar
Carl Schaffer committed
240
        return self._get_value(["STEPANGL"])
241
242
243
244

    @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
245
        return self._kw_mean(["FF1WLDSP", "FF2WLDSP"])
246
247
248
249
250

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

        # average for lower border:
Carl Schaffer's avatar
Carl Schaffer committed
254
255
        lower = self._kw_mean(["FF1WLOFF", "FF2WLOFF"])
        # add wl steps to calculate upper border
256
257
258
259
        upper = lower + step_size * n_steps_wl

        return lower, upper

260
261
262
263
264
265
266
267
268
    @property
    def derotator_in(self):
        """
        Determine whether the derotator is inserted

        Returns:
            bool: derotator in

        """
269
270
271
272
273
274
        try:
            rotval = self._get_value(["ROTCODE"])
        except KeyError:
            # if no derotator keywords are scontained, assume it doesn't exist yet
            return False

275
276
277
278
279
        if rotval == 0:
            return True
        else:
            return False

Carl Schaffer's avatar
Carl Schaffer committed
280
281
282
283
284
285
286
287
288
289
290
291
    @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")
292
        keyword = f"NAXIS{naxes}"
Carl Schaffer's avatar
Carl Schaffer committed
293
294
295
296
        states = self.query(keyword)

        return states

297
298
    @property
    def obs_mode(self):
299
        """Get observation mode"""
300
        self.parse()
301
302
303
304
        # extract obs mode keywords
        states = self.states
        n_maps = self.number_of_maps
        n_steps = self.n_steps
305

306
        mode = gris_obs_mode(states, n_maps, n_steps)
307
        return mode
Carl Schaffer's avatar
Carl Schaffer committed
308

Carl Schaffer's avatar
Carl Schaffer committed
309
310
311
    def parse(self):
        """extract data and metadata from file"""
        if self.is_parsed:
Carl Schaffer's avatar
Carl Schaffer committed
312
            return
313
314
315
316
317
318
319
320
321

        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
322
323

        # Read file
Carl Schaffer's avatar
Carl Schaffer committed
324
        try:
Carl Schaffer's avatar
Carl Schaffer committed
325
            header = self.header
326

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

Carl Schaffer's avatar
Carl Schaffer committed
330
            for key in header.keys():
331
332
333
                value = header[key]
                # check validity:
                value = value if is_valid(value) else ""
Carl Schaffer's avatar
Carl Schaffer committed
334
                if key in ["COMMENT", "HISTORY"]:
335
                    self.properties[key] = str(value)
Carl Schaffer's avatar
Carl Schaffer committed
336
                else:
337
338
                    self.properties[key] = value

Carl Schaffer's avatar
Carl Schaffer committed
339
            self.is_parsed = True
Carl Schaffer's avatar
Carl Schaffer committed
340
        except:
Carl Schaffer's avatar
linting    
Carl Schaffer committed
341
            exception = sys.exc_info()[0]
Carl Schaffer's avatar
Carl Schaffer committed
342
            self.errors.append("File %s" % self.filename + str(exception))
Carl Schaffer's avatar
Carl Schaffer committed
343
344

    def make_bson(self):
Carl Schaffer's avatar
Carl Schaffer committed
345
346
        """construct BSON document for mongodb storage"""
        return self.properties
Carl Schaffer's avatar
Carl Schaffer committed
347

348
    @property
349
    def _coords_from_wcs(self):
350
351
352
353
354
355
        """
        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
        """

Carl Schaffer's avatar
Carl Schaffer committed
356
357
358
        assert (
            not self.old_header
        ), f"Can't get coordinates from {self.path} ! Run header conversion first!"
359

360
361
362
        wcs = WCS(self.header)
        pixel_shape = wcs.pixel_shape
        n_y = pixel_shape[1]
363
        slit_list = np.zeros((n_y, self["NAXIS"]))
364
365
        slit_list[:, 1] = np.arange(n_y)
        slit_list[:, 0] = 0
366
        converted = np.array(wcs.pixel_to_world_values(slit_list))
367
        degrees = Angle(u.degree * converted[:, :2]).wrap_at(180 * u.degree)
368
        arcsec = degrees.to(u.arcsec)
369

370
371
        return arcsec

372
    def plot_slit(self):
373
        m = self.full_disk_map
374

375
        # retrieve different data sources
376
        coordinates = self.coords
377
378
379
380

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

        # determine fov limits for plotting
381

382
        plotting_padding = 300 * u.arcsec
383
384
385
386
        x_min = min(X) - plotting_padding
        x_max = max(X) + plotting_padding
        y_min = min(Y) - plotting_padding
        y_max = max(Y) + plotting_padding
387
388

        # Create submap of FOV and plot
Carl Schaffer's avatar
Carl Schaffer committed
389
390
391
392
        submap = m.submap(
            SkyCoord(x_min, y_min, frame=m.coordinate_frame),
            SkyCoord(x_max, y_max, frame=m.coordinate_frame),
        )
393
394
395
396
397
398
        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)
Carl Schaffer's avatar
Carl Schaffer committed
399
        ax.plot_coord(slit_wcs, "g-", label="WCS Coordinates")
400
401
402

        # Create Legend and annotations
        title = "Gris Slit {0}_{1:03d}_{2:03d}_{3:03d}"
Carl Schaffer's avatar
Carl Schaffer committed
403
404
405
406
407
408
        title = title.format(
            self.obs_time.strftime("%Y%m%d"),
            self.runnumber,
            self.mapnumber,
            self.step_index,
        )
409
410
        plt.title(title)

Carl Schaffer's avatar
Carl Schaffer committed
411
412
413
414
415
416
417
418
419
        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
        )
420
421
        plt.tight_layout(pad=3.5)
        _ = plt.legend()
422
423
424

        return fig, ax

Carl Schaffer's avatar
Carl Schaffer committed
425
    def __str__(self):
Carl Schaffer's avatar
Carl Schaffer committed
426
427
428
        outstring = (
            "GrisFitsFile %s from %s:\n\tRun: %s\n\tMap: %s\n\tSlitposition: %s\n"
        )
Carl Schaffer's avatar
linting    
Carl Schaffer committed
429
        outstring = outstring % (
Carl Schaffer's avatar
Carl Schaffer committed
430
            os.path.basename(self.filename),
Carl Schaffer's avatar
Carl Schaffer committed
431
432
            self.date,
            self.runnumber,
433
            self.map_index,
Carl Schaffer's avatar
Carl Schaffer committed
434
            self.slitposnumber,
Carl Schaffer's avatar
Carl Schaffer committed
435
        )
Carl Schaffer's avatar
linting    
Carl Schaffer committed
436
        return outstring
437
438


Carl Schaffer's avatar
Carl Schaffer committed
439
if __name__ == "__main__":
440
441
442
    gf = GrisFitsFile(sys.argv[1])
    gf.parse()
    print(gf)