GrisFitsFile.py 17.7 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
from pathlib import Path
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
Carl Schaffer's avatar
Carl Schaffer committed
21

Carl Schaffer's avatar
Carl Schaffer committed
22
from kis_tools.generic.fits import FitsFile
23
24
25
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
26
from kis_tools.util.util import date_from_fn, gris_obs_mode
Carl Schaffer's avatar
Carl Schaffer committed
27

28
29
logger = logging.getLogger(__name__)

Carl Schaffer's avatar
Carl Schaffer committed
30

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

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

    Returns:

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

Carl Schaffer's avatar
Carl Schaffer committed
48
49
    def __init__(self, filename):
        """Constructor"""
50
51
        if isinstance(filename, Path):
            filename = str(filename.resolve())
52
        super().__init__(filename)
Carl Schaffer's avatar
Carl Schaffer committed
53
        assert exists(filename)
Carl Schaffer's avatar
Carl Schaffer committed
54
55
56
57
        self.properties = {}
        self.errors = []

        self.filename = filename
58
        self.name = basename(filename)
59
        self.path = self.filename
60

61
        self.dirname = dirname(self.path)
Carl Schaffer's avatar
Carl Schaffer committed
62
63
64
65
        self.parse_filename()
        self.is_parsed = False

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

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

Carl Schaffer's avatar
linting    
Carl Schaffer committed
78
        # self.date = datetime.datetime(year,month,day,hh,mm,outstring)
Carl Schaffer's avatar
Carl Schaffer committed
79
        self.date = datetime.datetime(year, month, day, 0, 0, 0)
Carl Schaffer's avatar
Carl Schaffer committed
80

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

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

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

99
100
101
102
103
104
105
106
107
108
109
110
111
    @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
112
113
            x_telcenter, x_valid, x_age = query_gdbs('xsuncenter', self.obs_time)
            y_telcenter, y_valid, y_age = query_gdbs('ysuncenter', self.obs_time)
114
115
116
117
118
119
120
121
122
123
            # 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
124
125
        except ValueError:
            gdbs_usable = False
126

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

134
135
136
137
138
139
140
        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

141
        return x_telcenter, y_telcenter, x_std, y_std
142

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

162
        # Slit center in telescope coordinates
163
164
        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
165

166
167
168
169
170
        # 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
171

172
        if std_delta_x > coord_model_stdx or std_delta_y > coord_model_stdy:
173
174
175
176
177
178
            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
179
180
        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
181

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

190
191
192
193
194
        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
195

196
197
198
        # 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
199
200
201
202
203

        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)

204
205
        # track coordinate uncertainty within class
        self.coord_uncertainty = (std_delta_x, std_delta_y)
206
        return coord_array
207

208
209
210
211
212
213
214
215
    @property
    def coords(self):
        try:
            coords = self._coords_from_wcs
        except AssertionError:
            coords = self._coords_from_simple_header
        return coords

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

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

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

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

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

248
249
250
251
252
253
254
        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:
255
            rotangle = self.ROTANGLE
256
257
258
259
            alpha = -2 * rotangle + (const_term / 2)
        angle = slitorie - alpha

        return angle
Carl Schaffer's avatar
Carl Schaffer committed
260
261

    @property
262
263
    def number_of_maps(self):
        """
264
265

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

    @property
    def map_index(self):
272
273
274
275
276
        """

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

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

        if max_steps < nsteps:
            nsteps = max_steps

        return nsteps
Carl Schaffer's avatar
Carl Schaffer committed
304

305
306
    @property
    def obs_time(self):
307
308
309
310
311
        """

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

Carl Schaffer's avatar
Carl Schaffer committed
314
315
    @property
    def step_index(self):
316
317
318
319
320
        """

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

323
324
325
326
327
328
329
    @property
    def p0_angle(self):
        """

        Returns:
            p0_angle: angle between telescope and solar north
        """
330
331
332
333
        angle = self._get_value(["P0ANGLE", 'SOLAR_P0'])
        if angle == '':
            raise ValueError(f'Header of {self} contains illegal p0_angle:"{angle}"')
        return angle
334

335
336
    @property
    def step_angle(self):
337
338
339
340
341
342
        """

        Returns:
            step_angle: scanning angle
        """

Carl Schaffer's avatar
Carl Schaffer committed
343
        return self._get_value(["STEPANGL"])
344
345
346
347

    @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
348
        return self._kw_mean(["FF1WLDSP", "FF2WLDSP"])
349
350
351
352
353

    @property
    def wavelength_region(self):
        """Get wavelength region covered by file. Returns tuple of (upper, lower)"""
        step_size = self.wavelength_step_size
354

355
        n_steps_wl = self.data.shape[-1] if self.old_header else self.NAXIS3
356
357

        # average for lower border:
Carl Schaffer's avatar
Carl Schaffer committed
358
359
        lower = self._kw_mean(["FF1WLOFF", "FF2WLOFF"])
        # add wl steps to calculate upper border
360
361
362
363
        upper = lower + step_size * n_steps_wl

        return lower, upper

364
365
366
367
368
369
370
371
372
    @property
    def derotator_in(self):
        """
        Determine whether the derotator is inserted

        Returns:
            bool: derotator in

        """
373
374
375
376
377
378
        try:
            rotval = self._get_value(["ROTCODE"])
        except KeyError:
            # if no derotator keywords are scontained, assume it doesn't exist yet
            return False

379
380
381
382
383
        if rotval == 0:
            return True
        else:
            return False

Carl Schaffer's avatar
Carl Schaffer committed
384
385
386
387
388
389
390
391
392
393
394
395
    @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")
396
        keyword = f"NAXIS{naxes}"
Carl Schaffer's avatar
Carl Schaffer committed
397
398
399
400
        states = self.query(keyword)

        return states

401
402
    @property
    def obs_mode(self):
403
        """Get observation mode"""
404
        self.parse()
405
406
407
408
        # extract obs mode keywords
        states = self.states
        n_maps = self.number_of_maps
        n_steps = self.n_steps
409

410
        mode = gris_obs_mode(states, n_maps, n_steps)
411
        return mode
Carl Schaffer's avatar
Carl Schaffer committed
412

Carl Schaffer's avatar
Carl Schaffer committed
413
414
415
    def parse(self):
        """extract data and metadata from file"""
        if self.is_parsed:
Carl Schaffer's avatar
Carl Schaffer committed
416
            return
417
418
419
420
421
422
423
424
425

        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
426
427

        # Read file
Carl Schaffer's avatar
Carl Schaffer committed
428
        try:
Carl Schaffer's avatar
Carl Schaffer committed
429
            header = self.header
430

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

Carl Schaffer's avatar
Carl Schaffer committed
434
            for key in header.keys():
435
436
437
                value = header[key]
                # check validity:
                value = value if is_valid(value) else ""
Carl Schaffer's avatar
Carl Schaffer committed
438
                if key in ["COMMENT", "HISTORY"]:
439
                    self.properties[key] = str(value)
Carl Schaffer's avatar
Carl Schaffer committed
440
                else:
441
442
                    self.properties[key] = value

Carl Schaffer's avatar
Carl Schaffer committed
443
            self.is_parsed = True
Carl Schaffer's avatar
Carl Schaffer committed
444
        except:
Carl Schaffer's avatar
linting    
Carl Schaffer committed
445
            exception = sys.exc_info()[0]
Carl Schaffer's avatar
Carl Schaffer committed
446
            self.errors.append("File %s" % self.filename + str(exception))
Carl Schaffer's avatar
Carl Schaffer committed
447
448

    def make_bson(self):
Carl Schaffer's avatar
Carl Schaffer committed
449
450
        """construct BSON document for mongodb storage"""
        return self.properties
Carl Schaffer's avatar
Carl Schaffer committed
451

452
    @property
453
    def _coords_from_wcs(self):
454
455
456
457
        """
        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
458
459
460

        Side-effect: Sets self.coord_uncertainty

461
462
        """

Carl Schaffer's avatar
Carl Schaffer committed
463
464
465
        assert (
            not self.old_header
        ), f"Can't get coordinates from {self.path} ! Run header conversion first!"
466

467
468
469
        wcs = WCS(self.header)
        pixel_shape = wcs.pixel_shape
        n_y = pixel_shape[1]
470
        slit_list = np.zeros((n_y, self["NAXIS"]))
471
472
        slit_list[:, 1] = np.arange(n_y)
        slit_list[:, 0] = 0
473
        converted = np.array(wcs.pixel_to_world_values(slit_list))
474
        degrees = Angle(u.degree * converted[:, :2]).wrap_at(180 * u.degree)
475
        arcsec = degrees.to(u.arcsec)
476

477
478
479
        # Track coord uncertainty within instance
        self.coord_uncertainty = (self.CSYER1, self.CSYER2)

480
481
        return arcsec

482
    def plot_slit(self):
483
        m = self.full_disk_map
484

485
        # retrieve different data sources
486
        coordinates = self.coords
487
488
489
490

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

        # determine fov limits for plotting
491

492
        plotting_padding = 300 * u.arcsec
493
494
495
496
        x_min = min(X) - plotting_padding
        x_max = max(X) + plotting_padding
        y_min = min(Y) - plotting_padding
        y_max = max(Y) + plotting_padding
497
498

        # Create submap of FOV and plot
Carl Schaffer's avatar
Carl Schaffer committed
499
500
501
502
        submap = m.submap(
            SkyCoord(x_min, y_min, frame=m.coordinate_frame),
            SkyCoord(x_max, y_max, frame=m.coordinate_frame),
        )
503
504
505
506
507
508
        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
509
        ax.plot_coord(slit_wcs, "g-", label="WCS Coordinates")
510
511
512

        # Create Legend and annotations
        title = "Gris Slit {0}_{1:03d}_{2:03d}_{3:03d}"
Carl Schaffer's avatar
Carl Schaffer committed
513
514
515
516
517
518
        title = title.format(
            self.obs_time.strftime("%Y%m%d"),
            self.runnumber,
            self.mapnumber,
            self.step_index,
        )
519
520
        plt.title(title)

Carl Schaffer's avatar
Carl Schaffer committed
521
522
523
524
525
526
527
528
529
        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
        )
530
531
        plt.tight_layout(pad=3.5)
        _ = plt.legend()
532
533
534

        return fig, ax

Carl Schaffer's avatar
Carl Schaffer committed
535
    def __str__(self):
Carl Schaffer's avatar
Carl Schaffer committed
536
537
538
        outstring = (
            "GrisFitsFile %s from %s:\n\tRun: %s\n\tMap: %s\n\tSlitposition: %s\n"
        )
Carl Schaffer's avatar
linting    
Carl Schaffer committed
539
        outstring = outstring % (
Carl Schaffer's avatar
Carl Schaffer committed
540
            os.path.basename(self.filename),
Carl Schaffer's avatar
Carl Schaffer committed
541
542
            self.date,
            self.runnumber,
543
            self.map_index,
Carl Schaffer's avatar
Carl Schaffer committed
544
            self.slitposnumber,
Carl Schaffer's avatar
Carl Schaffer committed
545
        )
Carl Schaffer's avatar
linting    
Carl Schaffer committed
546
        return outstring
547
548


Carl Schaffer's avatar
Carl Schaffer committed
549
if __name__ == "__main__":
550
551
552
    gf = GrisFitsFile(sys.argv[1])
    gf.parse()
    print(gf)