GrisFitsFile.py 13 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, Angle
19
from astropy.wcs import WCS
Carl Schaffer's avatar
Carl Schaffer committed
20

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

Carl Schaffer's avatar
Carl Schaffer committed
24

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

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

    Returns:

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

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

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

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

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

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

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

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

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

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

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

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

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

121
122
        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
123

Carl Schaffer's avatar
Carl Schaffer committed
124
125
        message = (
            "Using FITS file coordinates, they are probably wrong by at least 100 arcsec!"
126
            " Also the slit angle is definitely off!"
Carl Schaffer's avatar
Carl Schaffer committed
127
        )
128
129
130
131
132
133
134
        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
135

136
137
138
139
140
141
142
143
    @property
    def coords(self):
        try:
            coords = self._coords_from_wcs
        except AssertionError:
            coords = self._coords_from_simple_header
        return coords

144
145
146
147
    @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
148
        if any([k in self.header for k in ["SOLARNET"]]):
149
150
            old_header = False
        return old_header
Carl Schaffer's avatar
Carl Schaffer committed
151

Carl Schaffer's avatar
Carl Schaffer committed
152
153
    @property
    def slit_orientation(self):
154
155
156
157
158
        """Get slit_orientation, defaults to 0. Not sure what this angle means...

        Returns:
            slit_orientation: float
        """
159
160
161
162

        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
163
164

    @property
165
166
    def number_of_maps(self):
        """
167
168

        Returns:
169
            nmaps: int number of maps in associated observation
170
        """
Carl Schaffer's avatar
Carl Schaffer committed
171
        return self._get_value(["NMAPS", "SERIES"])
Carl Schaffer's avatar
Carl Schaffer committed
172
173
174

    @property
    def map_index(self):
175
176
177
178
179
        """

        Returns:
            i_map: int map index within observation
        """
Carl Schaffer's avatar
Carl Schaffer committed
180
        return self._get_value(["IMAP", "ISERIE"])
Carl Schaffer's avatar
Carl Schaffer committed
181
182
183

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

        if max_steps < nsteps:
            nsteps = max_steps

        return nsteps
Carl Schaffer's avatar
Carl Schaffer committed
207

208
209
    @property
    def obs_time(self):
210
211
212
213
214
        """

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

Carl Schaffer's avatar
Carl Schaffer committed
217
218
    @property
    def step_index(self):
219
220
221
222
223
        """

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

226
227
    @property
    def step_angle(self):
228
229
230
231
232
233
        """

        Returns:
            step_angle: scanning angle
        """

Carl Schaffer's avatar
Carl Schaffer committed
234
        return self._get_value(["STEPANGL"])
235
236
237
238

    @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
239
        return self._kw_mean(["FF1WLDSP", "FF2WLDSP"])
240
241
242
243
244

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

        # average for lower border:
Carl Schaffer's avatar
Carl Schaffer committed
248
249
        lower = self._kw_mean(["FF1WLOFF", "FF2WLOFF"])
        # add wl steps to calculate upper border
250
251
252
253
        upper = lower + step_size * n_steps_wl

        return lower, upper

254
255
256
257
258
259
260
261
262
    @property
    def derotator_in(self):
        """
        Determine whether the derotator is inserted

        Returns:
            bool: derotator in

        """
263
264
265
266
267
268
        try:
            rotval = self._get_value(["ROTCODE"])
        except KeyError:
            # if no derotator keywords are scontained, assume it doesn't exist yet
            return False

269
270
271
272
273
        if rotval == 0:
            return True
        else:
            return False

Carl Schaffer's avatar
Carl Schaffer committed
274
275
276
277
278
279
280
281
282
283
284
285
    @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")
286
        keyword = f"NAXIS{naxes}"
Carl Schaffer's avatar
Carl Schaffer committed
287
288
289
290
        states = self.query(keyword)

        return states

291
292
    @property
    def obs_mode(self):
293
        """Get observation mode"""
294
        self.parse()
295
296
297
298
        # extract obs mode keywords
        states = self.states
        n_maps = self.number_of_maps
        n_steps = self.n_steps
299

300
        mode = gris_obs_mode(states, n_maps, n_steps)
301
        return mode
Carl Schaffer's avatar
Carl Schaffer committed
302

Carl Schaffer's avatar
Carl Schaffer committed
303
304
305
    def parse(self):
        """extract data and metadata from file"""
        if self.is_parsed:
Carl Schaffer's avatar
Carl Schaffer committed
306
            return
307
308
309
310
311
312
313
314
315

        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
316
317

        # Read file
Carl Schaffer's avatar
Carl Schaffer committed
318
        try:
Carl Schaffer's avatar
Carl Schaffer committed
319
            header = self.header
320

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

Carl Schaffer's avatar
Carl Schaffer committed
324
            for key in header.keys():
325
326
327
                value = header[key]
                # check validity:
                value = value if is_valid(value) else ""
Carl Schaffer's avatar
Carl Schaffer committed
328
                if key in ["COMMENT", "HISTORY"]:
329
                    self.properties[key] = str(value)
Carl Schaffer's avatar
Carl Schaffer committed
330
                else:
331
332
                    self.properties[key] = value

Carl Schaffer's avatar
Carl Schaffer committed
333
            self.is_parsed = True
Carl Schaffer's avatar
Carl Schaffer committed
334
        except:
Carl Schaffer's avatar
linting    
Carl Schaffer committed
335
            exception = sys.exc_info()[0]
Carl Schaffer's avatar
Carl Schaffer committed
336
            self.errors.append("File %s" % self.filename + str(exception))
Carl Schaffer's avatar
Carl Schaffer committed
337
338

    def make_bson(self):
Carl Schaffer's avatar
Carl Schaffer committed
339
340
        """construct BSON document for mongodb storage"""
        return self.properties
Carl Schaffer's avatar
Carl Schaffer committed
341

342
    @property
343
    def _coords_from_wcs(self):
344
345
346
347
348
349
        """
        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
350
351
352
        assert (
            not self.old_header
        ), f"Can't get coordinates from {self.path} ! Run header conversion first!"
353

354
355
356
        wcs = WCS(self.header)
        pixel_shape = wcs.pixel_shape
        n_y = pixel_shape[1]
357
        slit_list = np.zeros((n_y, self["NAXIS"]))
358
359
        slit_list[:, 1] = np.arange(n_y)
        slit_list[:, 0] = 0
360
        converted = np.array(wcs.pixel_to_world_values(slit_list))
361
        degrees = Angle(u.degree * converted[:, :2]).wrap_at(180 * u.degree)
362
        arcsec = degrees.to(u.arcsec)
363

364
365
        return arcsec

366
    def plot_slit(self):
367
        m = self.full_disk_map
368

369
        # retrieve different data sources
370
        coordinates = self.coords
371
372
373
374

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

        # determine fov limits for plotting
375

376
        plotting_padding = 300 * u.arcsec
377
378
379
380
        x_min = min(X) - plotting_padding
        x_max = max(X) + plotting_padding
        y_min = min(Y) - plotting_padding
        y_max = max(Y) + plotting_padding
381
382

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

        # Create Legend and annotations
        title = "Gris Slit {0}_{1:03d}_{2:03d}_{3:03d}"
Carl Schaffer's avatar
Carl Schaffer committed
397
398
399
400
401
402
        title = title.format(
            self.obs_time.strftime("%Y%m%d"),
            self.runnumber,
            self.mapnumber,
            self.step_index,
        )
403
404
        plt.title(title)

Carl Schaffer's avatar
Carl Schaffer committed
405
406
407
408
409
410
411
412
413
        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
        )
414
415
        plt.tight_layout(pad=3.5)
        _ = plt.legend()
416
417
418

        return fig, ax

Carl Schaffer's avatar
Carl Schaffer committed
419
    def __str__(self):
Carl Schaffer's avatar
Carl Schaffer committed
420
421
422
        outstring = (
            "GrisFitsFile %s from %s:\n\tRun: %s\n\tMap: %s\n\tSlitposition: %s\n"
        )
Carl Schaffer's avatar
linting    
Carl Schaffer committed
423
        outstring = outstring % (
Carl Schaffer's avatar
Carl Schaffer committed
424
            os.path.basename(self.filename),
Carl Schaffer's avatar
Carl Schaffer committed
425
426
            self.date,
            self.runnumber,
427
            self.map_index,
Carl Schaffer's avatar
Carl Schaffer committed
428
            self.slitposnumber,
Carl Schaffer's avatar
Carl Schaffer committed
429
        )
Carl Schaffer's avatar
linting    
Carl Schaffer committed
430
        return outstring
431
432


Carl Schaffer's avatar
Carl Schaffer committed
433
if __name__ == "__main__":
434
435
436
    gf = GrisFitsFile(sys.argv[1])
    gf.parse()
    print(gf)