GrisFitsFile.py 17.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
import logging
9
10
import os  # path and file manipulations
import sys  # error handling
11
from glob import glob
12
from os.path import basename, exists, join, dirname
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

Carl Schaffer's avatar
Carl Schaffer committed
21
from kis_tools.generic.fits import FitsFile
22
23
24
from kis_tools.gris.coord_correction_ml import get_coords_ml
from kis_tools.gris.util import gregor_parallactic_angle, telescope_to_hpc
from kis_tools.util.calculations import estimate_telescope_drift
Carl Schaffer's avatar
Carl Schaffer committed
25
from kis_tools.util.util import date_from_fn, gris_obs_mode
Carl Schaffer's avatar
Carl Schaffer committed
26

27
28
logger = logging.getLogger(__name__)

Carl Schaffer's avatar
Carl Schaffer committed
29

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

40
41
42
43
44
    Args:
      filename: filename of target fits file

    Returns:

Carl Schaffer's avatar
Carl Schaffer committed
45
    """
Carl Schaffer's avatar
Carl Schaffer committed
46

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

        self.filename = filename
55
        self.name = basename(filename)
56
        self.path = self.filename
57

58
        self.dirname = dirname(self.path)
Carl Schaffer's avatar
Carl Schaffer committed
59
60
61
62
        self.parse_filename()
        self.is_parsed = False

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

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

Carl Schaffer's avatar
linting    
Carl Schaffer committed
75
        # self.date = datetime.datetime(year,month,day,hh,mm,outstring)
Carl Schaffer's avatar
Carl Schaffer committed
76
        self.date = datetime.datetime(year, month, day, 0, 0, 0)
Carl Schaffer's avatar
Carl Schaffer committed
77

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

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

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

96
97
98
99
100
101
102
103
104
105
106
107
108
    @property
    def telescope_centering(self):
        """
        Determine state of telescope centering at the time of observation, queries GDBS
        for centering data

        Returns:
            outstring: x_offset, y_offset, x_uncertainty, y_uncertainty all in arcseconds
        """
        # Query centering from GDBS
        try:
            from kis_tools.gdbs.query_gdbs import query_gdbs
            gdbs_usable = True
109
110
            x_telcenter, x_valid, x_age = query_gdbs('xsuncenter', self.obs_time)
            y_telcenter, y_valid, y_age = query_gdbs('ysuncenter', self.obs_time)
111
112
113
114
115
116
117
118
119
120
            # Determine uncertainties for centering values
            centering_age = max(x_age, y_age)
            drift = estimate_telescope_drift(centering_age.seconds / 3600)
            x_std, y_std = drift, drift

            # if the centering is older than 24 hours, discard
            if centering_age > datetime.timedelta(hours=24):
                gdbs_usable = False
        except IOError:
            gdbs_usable = False
121
122
        except ValueError:
            gdbs_usable = False
123

124
        # Catch instances where no numbers but error strings are written to GDBS
Carl Schaffer's avatar
bugfix    
Carl Schaffer committed
125
126
127
128
129
        if gdbs_usable:
            try:
                x_telcenter, y_telcenter = float(x_telcenter), float(y_telcenter)
            except ValueError:
                gdbs_usable = False
130

131
132
133
134
135
136
137
        if not gdbs_usable:
            # Determined from offset between uncentered and cross_correlated data
            # STD-DEV of distances in X and Y again see gris coord study for details
            uncertainty_no_center_x = 117.1
            uncertainty_no_center_y = 182.4
            x_telcenter, y_telcenter, x_std, y_std = 0, 0, uncertainty_no_center_x, uncertainty_no_center_y

138
        return x_telcenter, y_telcenter, x_std, y_std
139

Carl Schaffer's avatar
Carl Schaffer committed
140
    @property
141
    def _coords_from_simple_header(self):
142
143
144
145
146
147
148
149
150
151
152
153
154
        """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
155
        h = self.header
156
        slit_length_pixels = h["NAXIS2"] if self.old_header else h['NAXIS1']
Carl Schaffer's avatar
Carl Schaffer committed
157
158
        center_pixel = slit_length_pixels // 2

159
        # Slit center in telescope coordinates
160
161
        x_slitcenter = h["SLITPOSX"] if h["SLITPOSX"] != "" else 0
        y_slitcenter = h["SLITPOSY"] if h["SLITPOSY"] != "" else 0
Carl Schaffer's avatar
Carl Schaffer committed
162

163
164
165
166
167
        # Transform to HPC (shift by telescope centering, rotate by p0 angle)
        delta_x, delta_y, std_delta_x, std_delta_y = self.telescope_centering

        # If uncertainties are too large, use ML model to correct the coordinates
        from kis_tools.gris.coord_correction_ml import coord_model_stdx, coord_model_stdy
168

169
        if std_delta_x > coord_model_stdx or std_delta_y > coord_model_stdy:
170
171
172
173
174
175
            try:
                x_hpc, y_hpc, std_delta_x, std_delta_y = get_coords_ml(self)
            except ValueError as e:
                logger.warning(f"Error in finding ML coords for {self.__repr__()}:{e}")
                # Point to slitpos with 1 solar disk uncertainty if coord retrieval didn't work
                x_hpc, y_hpc, std_delta_x, std_delta_y = x_slitcenter, y_slitcenter, 1000, 1000
176
177
        else:
            x_hpc, y_hpc = telescope_to_hpc(x_slitcenter, y_slitcenter, self.p0_angle, delta_x, delta_y)
Carl Schaffer's avatar
Carl Schaffer committed
178

179
        # Setup grid of pixels
Carl Schaffer's avatar
Carl Schaffer committed
180
181
        delta_pix = [i - center_pixel for i in range(slit_length_pixels)]
        delta_pix = np.array(delta_pix)
182
        arcsec_per_pixel_x = h["STEPSIZE"]  # assuming square pixels
Carl Schaffer's avatar
Carl Schaffer committed
183
184
185
        arcsec_per_pixel_y = (
            h["STEPSIZE"] if h["STEPSIZE"] else 0.135
        )  # assuming square pixels
Carl Schaffer's avatar
Carl Schaffer committed
186

187
188
189
190
191
        try:
            angle = self.slit_orientation
        except ValueError as e:
            logger.warning(f"Could not determine Slit orientation! \n {e}")
            angle = 0
Carl Schaffer's avatar
Carl Schaffer committed
192

193
194
195
        # generate coordinates for full grid of slit pixels
        X = x_hpc - np.sin(angle) * delta_pix * arcsec_per_pixel_x
        Y = y_hpc - np.cos(angle) * delta_pix * arcsec_per_pixel_y
196
197
198
199
200

        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)

201
202
        # track coordinate uncertainty within class
        self.coord_uncertainty = (std_delta_x, std_delta_y)
203
        return coord_array
204

205
206
207
208
209
210
211
212
    @property
    def coords(self):
        try:
            coords = self._coords_from_wcs
        except AssertionError:
            coords = self._coords_from_simple_header
        return coords

213
214
215
216
    @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
217
        if any([k in self.header for k in ["SOLARNET"]]):
218
219
            old_header = False
        return old_header
Carl Schaffer's avatar
Carl Schaffer committed
220

221
222
223
224
225
    @property
    def zenit(self):
        """
        Get Zenit angle (90° - elevation)
        Returns: zenit float in degrees
226
        Raises Value error if elevation is not contained in the header
227
228

        """
229
230
231
232
        elev = self._get_value(['ELEV_ANG', 'ELEVATIO'])
        if elev == '':
            raise ValueError(f'No elevation tracked for {self.__repr__()}')
        zenit = 90 - elev
233
234
        return zenit

Carl Schaffer's avatar
Carl Schaffer committed
235
236
    @property
    def slit_orientation(self):
237
238
        """Get slitorie:
         clockwise angle from the western section of the solar equator
239

240
        Attempts to correct for derotator and image rotation, usually correct within a couple 10s of degrees.
241
        Returns:
242
            slitorie: float (degrees)
243
        """
244

245
246
247
248
249
250
251
        p0_angle = self.p0_angle
        azimut = self._get_value(['AZIMUT'])
        zenit = self.zenit
        slitorie = self.SLITORIE
        const_term = -130  # determined by trial and error to match results from cross correlation with HMI
        alpha = -1 * gregor_parallactic_angle(self.obs_time) + azimut - zenit - p0_angle + const_term
        if self.derotator_in:
252
            rotangle = self.ROTANGLE
253
254
255
256
            alpha = -2 * rotangle + (const_term / 2)
        angle = slitorie - alpha

        return angle
Carl Schaffer's avatar
Carl Schaffer committed
257
258

    @property
259
260
    def number_of_maps(self):
        """
261
262

        Returns:
263
            nmaps: int number of maps in associated observation
264
        """
Carl Schaffer's avatar
Carl Schaffer committed
265
        return self._get_value(["NMAPS", "SERIES"])
Carl Schaffer's avatar
Carl Schaffer committed
266
267
268

    @property
    def map_index(self):
269
270
271
272
273
        """

        Returns:
            i_map: int map index within observation
        """
Carl Schaffer's avatar
Carl Schaffer committed
274
        return self._get_value(["IMAP", "ISERIE"])
Carl Schaffer's avatar
Carl Schaffer committed
275
276
277

    @property
    def n_steps(self):
278
279
280
281
282
283
284
285
286
287
288
289
290
        """
        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
291
        pattern = f"gris_*_*_*_{self.runnumber:03d}_{self.mapnumber:03d}_*.fits"
292
293
        matching_files = glob(join(self.dirname, pattern))
        max_steps = int(
Carl Schaffer's avatar
Carl Schaffer committed
294
            max([basename(m).split(".")[0].split("_")[-1] for m in matching_files])
295
296
297
298
299
300
        )

        if max_steps < nsteps:
            nsteps = max_steps

        return nsteps
Carl Schaffer's avatar
Carl Schaffer committed
301

302
303
    @property
    def obs_time(self):
304
305
306
307
308
        """

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

Carl Schaffer's avatar
Carl Schaffer committed
311
312
    @property
    def step_index(self):
313
314
315
316
317
        """

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

320
321
322
323
324
325
326
327
328
    @property
    def p0_angle(self):
        """

        Returns:
            p0_angle: angle between telescope and solar north
        """
        return self._get_value(["P0ANGLE", 'SOLAR_P0'])

329
330
    @property
    def step_angle(self):
331
332
333
334
335
336
        """

        Returns:
            step_angle: scanning angle
        """

Carl Schaffer's avatar
Carl Schaffer committed
337
        return self._get_value(["STEPANGL"])
338
339
340
341

    @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
342
        return self._kw_mean(["FF1WLDSP", "FF2WLDSP"])
343
344
345
346
347

    @property
    def wavelength_region(self):
        """Get wavelength region covered by file. Returns tuple of (upper, lower)"""
        step_size = self.wavelength_step_size
348
        n_steps_wl = self.data.shape[-1] if self.old_header else self.NAXIS3
349
350

        # average for lower border:
Carl Schaffer's avatar
Carl Schaffer committed
351
352
        lower = self._kw_mean(["FF1WLOFF", "FF2WLOFF"])
        # add wl steps to calculate upper border
353
354
355
356
        upper = lower + step_size * n_steps_wl

        return lower, upper

357
358
359
360
361
362
363
364
365
    @property
    def derotator_in(self):
        """
        Determine whether the derotator is inserted

        Returns:
            bool: derotator in

        """
366
367
368
369
370
371
        try:
            rotval = self._get_value(["ROTCODE"])
        except KeyError:
            # if no derotator keywords are scontained, assume it doesn't exist yet
            return False

372
373
374
375
376
        if rotval == 0:
            return True
        else:
            return False

Carl Schaffer's avatar
Carl Schaffer committed
377
378
379
380
381
382
383
384
385
386
387
388
    @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")
389
        keyword = f"NAXIS{naxes}"
Carl Schaffer's avatar
Carl Schaffer committed
390
391
392
393
        states = self.query(keyword)

        return states

394
395
    @property
    def obs_mode(self):
396
        """Get observation mode"""
397
        self.parse()
398
399
400
401
        # extract obs mode keywords
        states = self.states
        n_maps = self.number_of_maps
        n_steps = self.n_steps
402

403
        mode = gris_obs_mode(states, n_maps, n_steps)
404
        return mode
Carl Schaffer's avatar
Carl Schaffer committed
405

Carl Schaffer's avatar
Carl Schaffer committed
406
407
408
    def parse(self):
        """extract data and metadata from file"""
        if self.is_parsed:
Carl Schaffer's avatar
Carl Schaffer committed
409
            return
410
411
412
413
414
415
416
417
418

        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
419
420

        # Read file
Carl Schaffer's avatar
Carl Schaffer committed
421
        try:
Carl Schaffer's avatar
Carl Schaffer committed
422
            header = self.header
423

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

Carl Schaffer's avatar
Carl Schaffer committed
427
            for key in header.keys():
428
429
430
                value = header[key]
                # check validity:
                value = value if is_valid(value) else ""
Carl Schaffer's avatar
Carl Schaffer committed
431
                if key in ["COMMENT", "HISTORY"]:
432
                    self.properties[key] = str(value)
Carl Schaffer's avatar
Carl Schaffer committed
433
                else:
434
435
                    self.properties[key] = value

Carl Schaffer's avatar
Carl Schaffer committed
436
            self.is_parsed = True
Carl Schaffer's avatar
Carl Schaffer committed
437
        except:
Carl Schaffer's avatar
linting    
Carl Schaffer committed
438
            exception = sys.exc_info()[0]
Carl Schaffer's avatar
Carl Schaffer committed
439
            self.errors.append("File %s" % self.filename + str(exception))
Carl Schaffer's avatar
Carl Schaffer committed
440
441

    def make_bson(self):
Carl Schaffer's avatar
Carl Schaffer committed
442
443
        """construct BSON document for mongodb storage"""
        return self.properties
Carl Schaffer's avatar
Carl Schaffer committed
444

445
    @property
446
    def _coords_from_wcs(self):
447
448
449
450
451
452
        """
        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
453
454
455
        assert (
            not self.old_header
        ), f"Can't get coordinates from {self.path} ! Run header conversion first!"
456

457
458
459
        wcs = WCS(self.header)
        pixel_shape = wcs.pixel_shape
        n_y = pixel_shape[1]
460
        slit_list = np.zeros((n_y, self["NAXIS"]))
461
462
        slit_list[:, 1] = np.arange(n_y)
        slit_list[:, 0] = 0
463
        converted = np.array(wcs.pixel_to_world_values(slit_list))
464
        degrees = Angle(u.degree * converted[:, :2]).wrap_at(180 * u.degree)
465
        arcsec = degrees.to(u.arcsec)
466

467
468
469
        # Track coord uncertainty within instance
        self.coord_uncertainty = (self.CSYER1, self.CSYER2)

470
471
        return arcsec

472
    def plot_slit(self):
473
        m = self.full_disk_map
474

475
        # retrieve different data sources
476
        coordinates = self.coords
477
478
479
480

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

        # determine fov limits for plotting
481

482
        plotting_padding = 300 * u.arcsec
483
484
485
486
        x_min = min(X) - plotting_padding
        x_max = max(X) + plotting_padding
        y_min = min(Y) - plotting_padding
        y_max = max(Y) + plotting_padding
487
488

        # Create submap of FOV and plot
Carl Schaffer's avatar
Carl Schaffer committed
489
490
491
492
        submap = m.submap(
            SkyCoord(x_min, y_min, frame=m.coordinate_frame),
            SkyCoord(x_max, y_max, frame=m.coordinate_frame),
        )
493
494
495
496
497
498
        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
499
        ax.plot_coord(slit_wcs, "g-", label="WCS Coordinates")
500
501
502

        # Create Legend and annotations
        title = "Gris Slit {0}_{1:03d}_{2:03d}_{3:03d}"
Carl Schaffer's avatar
Carl Schaffer committed
503
504
505
506
507
508
        title = title.format(
            self.obs_time.strftime("%Y%m%d"),
            self.runnumber,
            self.mapnumber,
            self.step_index,
        )
509
510
        plt.title(title)

Carl Schaffer's avatar
Carl Schaffer committed
511
512
513
514
515
516
517
518
519
        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
        )
520
521
        plt.tight_layout(pad=3.5)
        _ = plt.legend()
522
523
524

        return fig, ax

Carl Schaffer's avatar
Carl Schaffer committed
525
    def __str__(self):
Carl Schaffer's avatar
Carl Schaffer committed
526
527
528
        outstring = (
            "GrisFitsFile %s from %s:\n\tRun: %s\n\tMap: %s\n\tSlitposition: %s\n"
        )
Carl Schaffer's avatar
linting    
Carl Schaffer committed
529
        outstring = outstring % (
Carl Schaffer's avatar
Carl Schaffer committed
530
            os.path.basename(self.filename),
Carl Schaffer's avatar
Carl Schaffer committed
531
532
            self.date,
            self.runnumber,
533
            self.map_index,
Carl Schaffer's avatar
Carl Schaffer committed
534
            self.slitposnumber,
Carl Schaffer's avatar
Carl Schaffer committed
535
        )
Carl Schaffer's avatar
linting    
Carl Schaffer committed
536
        return outstring
537
538


Carl Schaffer's avatar
Carl Schaffer committed
539
if __name__ == "__main__":
540
541
542
    gf = GrisFitsFile(sys.argv[1])
    gf.parse()
    print(gf)