GrisFitsFile.py 18.3 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
20
21
import pytz
from astroplan import Observer
from astropy.coordinates import SkyCoord, Angle, get_sun
22
from astropy.wcs import WCS
Carl Schaffer's avatar
Carl Schaffer committed
23

Carl Schaffer's avatar
Carl Schaffer committed
24
from kis_tools.generic.fits import FitsFile
25
from kis_tools.gris.coord_correction_ml import get_coords_ml
Carl Schaffer's avatar
Carl Schaffer committed
26
from kis_tools.gris.util import telescope_to_hpc, expand_coords, calculate_slit_translation
27
from kis_tools.util.calculations import estimate_telescope_drift
28
from kis_tools.util.constants import gregor_coords, gris_stepsize_along_slit
Carl Schaffer's avatar
Carl Schaffer committed
29
from kis_tools.util.util import date_from_fn, gris_obs_mode
Carl Schaffer's avatar
Carl Schaffer committed
30

31
32
logger = logging.getLogger(__name__)

Carl Schaffer's avatar
Carl Schaffer committed
33

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

44
45
46
47
48
    Args:
      filename: filename of target fits file

    Returns:

Carl Schaffer's avatar
Carl Schaffer committed
49
    """
Carl Schaffer's avatar
Carl Schaffer committed
50

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

        self.filename = filename
61
        self.name = basename(filename)
62
        self.path = self.filename
63

64
        self.dirname = dirname(self.path)
Carl Schaffer's avatar
Carl Schaffer committed
65
66
67
68
        self.parse_filename()
        self.is_parsed = False

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

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

Carl Schaffer's avatar
linting    
Carl Schaffer committed
81
        # self.date = datetime.datetime(year,month,day,hh,mm,outstring)
Carl Schaffer's avatar
Carl Schaffer committed
82
        self.date = datetime.datetime(year, month, day, 0, 0, 0)
Carl Schaffer's avatar
Carl Schaffer committed
83

84
        # Extract Run-, Map-, and Stepnumber, as well as mode string
85
        self.runnumber = int(name[4])
86
        self.mapnumber = int(name[5])
87
        self.slitposnumber = int(name[6])
88
89
        self.modestring = name[3]

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

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

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

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

137
138
139
140
141
142
143
        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

144
        return x_telcenter, y_telcenter, x_std, y_std
145

Carl Schaffer's avatar
Carl Schaffer committed
146
    @property
147
    def _coords_from_simple_header(self):
148
149
150
151
152
153
154
155
156
157
158
159
160
        """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
161
        h = self.header
Carl Schaffer's avatar
Carl Schaffer committed
162
        slit_length_pixels = self.NAXIS2
Carl Schaffer's avatar
Carl Schaffer committed
163

164
        # Slit center in telescope coordinates
Carl Schaffer's avatar
Carl Schaffer committed
165
166
167
168
169
170
171
        x_slitcenter = self.SLITPOSX if self.header['SLITPOSX'] != "" else 0
        y_slitcenter = self.SLITPOSX if self.header['SLITPOSY'] != "" else 0

        coord_ref = slit_length_pixels // 2
        if self.SLITCNTR.strip() != 'scancntr':
            raise ValueError(
                f'Cannot determine reference pixel for coordinates! SLITCNTR is {self.SLITCNTR} instead of scancntr')
Carl Schaffer's avatar
Carl Schaffer committed
172

173
174
175
176
177
        # 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
178

179
        if std_delta_x > coord_model_stdx or std_delta_y > coord_model_stdy:
180
181
182
183
184
185
            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
186
187
        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
188

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

195
196
197
198
199
        X, Y = expand_coords(x_hpc, y_hpc,
                             slit_orientation=angle,
                             stepsize_slit=gris_stepsize_along_slit,
                             slit_length_pixels=slit_length_pixels)

200
201
202
203
204
        coord_shift = calculate_slit_translation(slit_number=self.slitposnumber, stepsize_perp=self.STEPSIZE,
                                                 slit_orientation=angle)

        X = X + coord_shift[0]
        Y = Y + coord_shift[1]
205
206
207
208
209

        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)

210
211
        # track coordinate uncertainty within class
        self.coord_uncertainty = (std_delta_x, std_delta_y)
212
        return coord_array
213

214
215
216
217
218
219
220
221
    @property
    def coords(self):
        try:
            coords = self._coords_from_wcs
        except AssertionError:
            coords = self._coords_from_simple_header
        return coords

222
223
224
225
    @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
226
        if any([k in self.header for k in ["SOLARNET"]]):
227
228
            old_header = False
        return old_header
Carl Schaffer's avatar
Carl Schaffer committed
229

230
231
232
233
234
    @property
    def zenit(self):
        """
        Get Zenit angle (90° - elevation)
        Returns: zenit float in degrees
235
        Raises Value error if elevation is not contained in the header
236
237

        """
238
239
240
241
        elev = self._get_value(['ELEV_ANG', 'ELEVATIO'])
        if elev == '':
            raise ValueError(f'No elevation tracked for {self.__repr__()}')
        zenit = 90 - elev
242
243
        return zenit

Carl Schaffer's avatar
Carl Schaffer committed
244
245
    @property
    def slit_orientation(self):
246
247
        """Get slitorie:
         clockwise angle from the western section of the solar equator
248

249
        Attempts to correct for derotator and image rotation, usually correct within a couple 10s of degrees.
250
        Returns:
251
            slitorie: float (degrees)
252
        """
253

254
255
256
257
258
        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
259
260
261
262
263
264
265
266
267
268
269

        gregor_observer = Observer(location=gregor_coords)
        date_tz = pytz.utc.localize(self.obs_time)
        time = gregor_observer.datetime_to_astropy_time(date_tz)
        sun = get_sun(time)

        dec = sun.dec
        parallactic_angle = np.degrees(np.arcsin(np.sin(np.radians(azimut))*np.cos(gregor_coords.lat.to(u.rad))/np.cos(dec.to(u.rad))))
        self.parallactic_angle=parallactic_angle

        alpha = -1 * parallactic_angle.to(u.degree).value + azimut - zenit - p0_angle + const_term
270
        if self.derotator_in:
271
            rotangle = self.ROTANGLE
272
273
274
275
            alpha = -2 * rotangle + (const_term / 2)
        angle = slitorie - alpha

        return angle
Carl Schaffer's avatar
Carl Schaffer committed
276
277

    @property
278
279
    def number_of_maps(self):
        """
280
281

        Returns:
282
            nmaps: int number of maps in associated observation
283
        """
Carl Schaffer's avatar
Carl Schaffer committed
284
        return self._get_value(["NMAPS", "SERIES"])
Carl Schaffer's avatar
Carl Schaffer committed
285
286
287

    @property
    def map_index(self):
288
289
290
291
292
        """

        Returns:
            i_map: int map index within observation
        """
Carl Schaffer's avatar
Carl Schaffer committed
293
        return self._get_value(["IMAP", "ISERIE"])
Carl Schaffer's avatar
Carl Schaffer committed
294
295
296

    @property
    def n_steps(self):
297
298
299
300
301
302
303
304
305
306
307
308
309
        """
        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
310
        pattern = f"gris_*_*_*_{self.runnumber:03d}_{self.mapnumber:03d}_*.fits"
311
312
        matching_files = glob(join(self.dirname, pattern))
        max_steps = int(
Carl Schaffer's avatar
Carl Schaffer committed
313
            max([basename(m).split(".")[0].split("_")[-1] for m in matching_files])
314
315
316
317
318
319
        )

        if max_steps < nsteps:
            nsteps = max_steps

        return nsteps
Carl Schaffer's avatar
Carl Schaffer committed
320

321
322
    @property
    def obs_time(self):
323
324
325
326
327
        """

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

Carl Schaffer's avatar
Carl Schaffer committed
330
331
    @property
    def step_index(self):
332
333
334
335
336
        """

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

339
340
341
342
343
344
345
    @property
    def p0_angle(self):
        """

        Returns:
            p0_angle: angle between telescope and solar north
        """
346
347
348
349
        angle = self._get_value(["P0ANGLE", 'SOLAR_P0'])
        if angle == '':
            raise ValueError(f'Header of {self} contains illegal p0_angle:"{angle}"')
        return angle
350

351
352
    @property
    def step_angle(self):
353
354
355
356
357
358
        """

        Returns:
            step_angle: scanning angle
        """

Carl Schaffer's avatar
Carl Schaffer committed
359
        return self._get_value(["STEPANGL"])
360
361
362
363

    @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
364
        return self._kw_mean(["FF1WLDSP", "FF2WLDSP"])
365
366
367
368
369

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

371
        n_steps_wl = self.data.shape[-1] if self.old_header else self.NAXIS3
372
373

        # average for lower border:
Carl Schaffer's avatar
Carl Schaffer committed
374
375
        lower = self._kw_mean(["FF1WLOFF", "FF2WLOFF"])
        # add wl steps to calculate upper border
376
377
378
379
        upper = lower + step_size * n_steps_wl

        return lower, upper

380
381
382
383
384
385
386
387
388
    @property
    def derotator_in(self):
        """
        Determine whether the derotator is inserted

        Returns:
            bool: derotator in

        """
389
390
391
392
393
394
        try:
            rotval = self._get_value(["ROTCODE"])
        except KeyError:
            # if no derotator keywords are scontained, assume it doesn't exist yet
            return False

395
396
397
398
399
        if rotval == 0:
            return True
        else:
            return False

Carl Schaffer's avatar
Carl Schaffer committed
400
401
402
403
404
405
406
407
408
409
410
411
    @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")
412
        keyword = f"NAXIS{naxes}"
Carl Schaffer's avatar
Carl Schaffer committed
413
414
415
416
        states = self.query(keyword)

        return states

417
418
    @property
    def obs_mode(self):
419
        """Get observation mode"""
420
        self.parse()
421
422
423
424
        # extract obs mode keywords
        states = self.states
        n_maps = self.number_of_maps
        n_steps = self.n_steps
425

426
        mode = gris_obs_mode(states, n_maps, n_steps)
427
        return mode
Carl Schaffer's avatar
Carl Schaffer committed
428

Carl Schaffer's avatar
Carl Schaffer committed
429
430
431
    def parse(self):
        """extract data and metadata from file"""
        if self.is_parsed:
Carl Schaffer's avatar
Carl Schaffer committed
432
            return
433
434
435
436
437
438
439
440
441

        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
442
443

        # Read file
Carl Schaffer's avatar
Carl Schaffer committed
444
        try:
Carl Schaffer's avatar
Carl Schaffer committed
445
            header = self.header
446

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

Carl Schaffer's avatar
Carl Schaffer committed
450
            for key in header.keys():
451
452
453
                value = header[key]
                # check validity:
                value = value if is_valid(value) else ""
Carl Schaffer's avatar
Carl Schaffer committed
454
                if key in ["COMMENT", "HISTORY"]:
455
                    self.properties[key] = str(value)
Carl Schaffer's avatar
Carl Schaffer committed
456
                else:
457
458
                    self.properties[key] = value

Carl Schaffer's avatar
Carl Schaffer committed
459
            self.is_parsed = True
Carl Schaffer's avatar
Carl Schaffer committed
460
        except:
Carl Schaffer's avatar
linting    
Carl Schaffer committed
461
            exception = sys.exc_info()[0]
Carl Schaffer's avatar
Carl Schaffer committed
462
            self.errors.append("File %s" % self.filename + str(exception))
Carl Schaffer's avatar
Carl Schaffer committed
463
464

    def make_bson(self):
Carl Schaffer's avatar
Carl Schaffer committed
465
466
        """construct BSON document for mongodb storage"""
        return self.properties
Carl Schaffer's avatar
Carl Schaffer committed
467

468
    @property
469
    def _coords_from_wcs(self):
470
471
472
473
        """
        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
474
475
476

        Side-effect: Sets self.coord_uncertainty

477
478
        """

Carl Schaffer's avatar
Carl Schaffer committed
479
480
481
        assert (
            not self.old_header
        ), f"Can't get coordinates from {self.path} ! Run header conversion first!"
482

483
484
485
        wcs = WCS(self.header)
        pixel_shape = wcs.pixel_shape
        n_y = pixel_shape[1]
486
        slit_list = np.zeros((n_y, self["NAXIS"]))
487
488
        slit_list[:, 1] = np.arange(n_y)
        slit_list[:, 0] = 0
489
        converted = np.array(wcs.pixel_to_world_values(slit_list))
490
        degrees = Angle(u.degree * converted[:, :2]).wrap_at(180 * u.degree)
491
        arcsec = degrees.to(u.arcsec)
492

493
494
495
        # Track coord uncertainty within instance
        self.coord_uncertainty = (self.CSYER1, self.CSYER2)

496
497
        return arcsec

498
    def plot_slit(self):
499
        m = self.full_disk_map
500

501
        # retrieve different data sources
502
        coordinates = self.coords
503
504
505
506

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

        # determine fov limits for plotting
507

508
        plotting_padding = 300 * u.arcsec
509
510
511
512
        x_min = min(X) - plotting_padding
        x_max = max(X) + plotting_padding
        y_min = min(Y) - plotting_padding
        y_max = max(Y) + plotting_padding
513
514

        # Create submap of FOV and plot
Carl Schaffer's avatar
Carl Schaffer committed
515
516
517
518
        submap = m.submap(
            SkyCoord(x_min, y_min, frame=m.coordinate_frame),
            SkyCoord(x_max, y_max, frame=m.coordinate_frame),
        )
519
520
521
522
523
524
        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
525
        ax.plot_coord(slit_wcs, "g-", label="WCS Coordinates")
526
527
528

        # Create Legend and annotations
        title = "Gris Slit {0}_{1:03d}_{2:03d}_{3:03d}"
Carl Schaffer's avatar
Carl Schaffer committed
529
530
531
532
533
534
        title = title.format(
            self.obs_time.strftime("%Y%m%d"),
            self.runnumber,
            self.mapnumber,
            self.step_index,
        )
535
536
        plt.title(title)

Carl Schaffer's avatar
Carl Schaffer committed
537
538
539
540
541
542
543
544
545
        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
        )
546
547
        plt.tight_layout(pad=3.5)
        _ = plt.legend()
548
549
550

        return fig, ax

Carl Schaffer's avatar
Carl Schaffer committed
551
    def __str__(self):
Carl Schaffer's avatar
Carl Schaffer committed
552
553
554
        outstring = (
            "GrisFitsFile %s from %s:\n\tRun: %s\n\tMap: %s\n\tSlitposition: %s\n"
        )
Carl Schaffer's avatar
linting    
Carl Schaffer committed
555
        outstring = outstring % (
Carl Schaffer's avatar
Carl Schaffer committed
556
            os.path.basename(self.filename),
Carl Schaffer's avatar
Carl Schaffer committed
557
558
            self.date,
            self.runnumber,
559
            self.map_index,
Carl Schaffer's avatar
Carl Schaffer committed
560
            self.slitposnumber,
Carl Schaffer's avatar
Carl Schaffer committed
561
        )
Carl Schaffer's avatar
linting    
Carl Schaffer committed
562
        return outstring
563
564


Carl Schaffer's avatar
Carl Schaffer committed
565
if __name__ == "__main__":
566
567
568
    gf = GrisFitsFile(sys.argv[1])
    gf.parse()
    print(gf)