GrisRun.py 22.7 KB
Newer Older
Carl Schaffer's avatar
Carl Schaffer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
"""
Created on 2018-04-09
@author Carl Schaffer
@mail   carl.schaffer@leibniz-kis.de

class for retrieving run information for a set of GrisFitsFiles

If this file is called directly, it will attempt to create a GrisRun
instance for target passed as the first commandline argument. If no
argument is passed, it will attempt to process a dummy folder.
"""

import datetime
14
import re
15
16
from glob import glob  # get filenames in unix style
from os.path import join, dirname, abspath, isdir, basename
Carl Schaffer's avatar
Carl Schaffer committed
17

18
19
20
import astropy.io.fits as fitsio
import matplotlib.pyplot as plt
import numpy as np
21
22
23
from tqdm import tqdm

from kis_tools.util.locplot import make_loc_plot
24
25
from . import GrisFitsFile
from .util import get_observers
26
from ..generic import gris_fields
27
from ..util.util import reglob, gris_run_number, get_filter
Carl Schaffer's avatar
Carl Schaffer committed
28

Carl Schaffer's avatar
Carl Schaffer committed
29

Carl Schaffer's avatar
Carl Schaffer committed
30
31
class GrisRun(object):
    """Class to handle run-wise gris data
Carl Schaffer's avatar
Carl Schaffer committed
32

33
34
35
36
    Args:
      gris_fits_files: list of gris files to form one run

    Returns:
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
37
      ris_run: GrisRun instance
Carl Schaffer's avatar
Carl Schaffer committed
38
    """
Carl Schaffer's avatar
Carl Schaffer committed
39

40
    # Patterns for files:
41
42
    _l0pattern = r"(-\d\d){0,1}"
    _l1pattern = r"(?:-\d\d){0,1}c[cm]"
43

Carl Schaffer's avatar
Carl Schaffer committed
44
45
    # List of keywords to extract from header, values are passed
    # through self.query
46
47
48
49
50
51
52
    pars_native = {
        "BITPIX": "BITPIX",
        "NAXIS": "NAXIS",
        "BTYPE": "BTYPE",
        "DATE-OBS": "date",
        "TELESCOP": "TELESCOP",
        "CAMERA": "CAMERA",
53
        "WAVELENG": ["WAVELENG", "AWAVLNTH"],
54
55
56
57
        "STATES": "states",
        "AOSYSTEM": "AOSYSTEM",
        "AO_LOCK": "AO_LOCK",
        "NSUMEXP": "NSUMEXP",
58
        "SERIES": ["NSERIES", "NMAPS"],
59
        "IOBS": "IOBS",
60
        "RSUN_ARC": "RSUN_ARC",
61
62
        "STEPS": "NSTEPS",
        "STEPSIZE": "STEPSIZE",
Carl Schaffer's avatar
Carl Schaffer committed
63
64
        "EXPTIME": "XPOSURE",
        "TARGET": "OBS_TRGT",
65
        "OBS_MODE": "obs_mode",
66
    }
Carl Schaffer's avatar
Carl Schaffer committed
67

Carl Schaffer's avatar
Carl Schaffer committed
68
    # list of keywords where a range is deduced from the input files, values are passed through GrisFitsFile.query
69
70
71
72
73
74
75
76
    pars_minmax = {
        "UT": ["DATE-BEG"],
        "FRIEDR0": "ATMOS_R0",
        "AZIMUT": "AZIMUT",
        "ELEVATIO": "ELEV_ANG",
        "SLITPOSX": "SLITPOSX",
        "SLITPOSY": "SLITPOSY",
    }
Carl Schaffer's avatar
Carl Schaffer committed
77

Carl Schaffer's avatar
Carl Schaffer committed
78
    # Format string for storing revision dates in DB
Carl Schaffer's avatar
Carl Schaffer committed
79
80
    date_formatstring = "%Y-%m-%d %H:%M:%S"

Carl Schaffer's avatar
Carl Schaffer committed
81
    def __init__(self, gris_fits_files):
82
        # check inputs
Carl Schaffer's avatar
Carl Schaffer committed
83
        if not gris_fits_files:
84
            raise ValueError("Empty list of files passed to GrisRun")
Carl Schaffer's avatar
Carl Schaffer committed
85

Carl Schaffer's avatar
Carl Schaffer committed
86
        # Setup data containers
Carl Schaffer's avatar
Carl Schaffer committed
87
88
89
90
91
        self.properties = {}
        self.files_broken = []
        self.errors = []
        self.is_parsed = False

92
93
94
        # check inputs:
        if not isinstance(gris_fits_files, list):
            print("Invalid input to GrisRun!")
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
95
96
97
98
            raise TypeError

        # check if inputs are already instances of GrisFitsFile, if
        # not, generate instances
99
        elif isinstance(gris_fits_files[0], str):
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
100
101
102
103
104
105
            gris_fits_files = [
                GrisFitsFile.GrisFitsFile(gff) for gff in gris_fits_files
            ]

        # store files sorted by filename
        self.files = sorted(gris_fits_files, key=lambda x: x.filename)
106

Carl Schaffer's avatar
Carl Schaffer committed
107
        # Extract run-wise parameters from first file of the run
Carl Schaffer's avatar
Carl Schaffer committed
108
109
110
        self.runnumber = self.files[0].runnumber
        self.properties["RUNNUMBER"] = self.runnumber
        self.date = self.files[0].date
111
        self.path = abspath(join(dirname(self.files[0].path), ".."))
Carl Schaffer's avatar
Carl Schaffer committed
112
        self.properties["DATE-OBS"] = self.date
113
        self.maps = list(
Carl Schaffer's avatar
Carl Schaffer committed
114
115
            set([gris_fitsfile.mapnumber for gris_fitsfile in self.files])
        )
Carl Schaffer's avatar
Carl Schaffer committed
116
        self.map_lengths = {}
117
        self.get_map_lengths()
Carl Schaffer's avatar
minor    
Carl Schaffer committed
118

Carl Schaffer's avatar
Carl Schaffer committed
119
        # Check data levels:
Carl Schaffer's avatar
Carl Schaffer committed
120
        self.check_data_levels()
Carl Schaffer's avatar
Carl Schaffer committed
121
        # Look for preview images and log files
Carl Schaffer's avatar
Carl Schaffer committed
122
123
124
125
126
        self.get_previews()

    def parse(self):
        """Performs all information gathering that requires opening
        individual fits files instead of only using info from their
127
128
129
130
131
132
133
        filenames

        Args:

        Returns:

        """
Carl Schaffer's avatar
Carl Schaffer committed
134
        if self.is_parsed:
Carl Schaffer's avatar
Carl Schaffer committed
135
            return self.is_parsed
136
137
138
139
140
141
142
143
144

        progress_bar = tqdm(
            sorted(self.files, key=lambda x: x.filename, reverse=True)
        )
        for fits_file in progress_bar:
            progress_bar.set_description(
                f"Parsing {basename(fits_file.filename)}"
            )
            fits_file.parse()
Carl Schaffer's avatar
Carl Schaffer committed
145
            if fits_file.errors:
Carl Schaffer's avatar
Carl Schaffer committed
146
                # Get errors from file,move file to broken list
147
                for error in fits_file.errors:
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
148
                    self.errors.append(error)
149
150
151
                self.files_broken.append(
                    self.files.pop(self.files.index(fits_file))
                )
Carl Schaffer's avatar
Carl Schaffer committed
152
        if not self.files:
Carl Schaffer's avatar
Carl Schaffer committed
153
            self.errors.append(
Carl Schaffer's avatar
Carl Schaffer committed
154
155
                "ERROR in %s: No fits files for run!" % (str(self))
            )
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
156
            return False
Carl Schaffer's avatar
Carl Schaffer committed
157

158
159
160
        # Extract observation properites from files
        pass

Carl Schaffer's avatar
Carl Schaffer committed
161
        # Extract run-wise parameters from files
162
163
164
        files = self.files
        file = files[0]

165
166
        xmin, xmax, ymin, ymax = self.calculate_bounding_box()

167
        observation_description = dict(
168
            OBS_NAME=self.obs_name,
169
170
171
172
            DATASET_NAME=f"gris_{self.date.strftime('%Y%m%d')}_{self.runnumber:03d}",
            DATAPRODUCT_TYPE="cube",
            CALIB_LEVEL=1,
            OBS_COLLECTION="gris_observations",
Carl Schaffer's avatar
Carl Schaffer committed
173
            DATE_BEG=min(self.query("obs_time")),
174
            DATE_END=max(self.query("obs_time")) + datetime.timedelta(seconds=file.query("XPOSURE") / 1000.),
Carl Schaffer's avatar
Carl Schaffer committed
175
176
            DATE_OBS=min(self.query("obs_time")),
            TIME_RESOLUTION=(max(self.query("obs_time")) - min(self.query("obs_time"))).seconds / file.query("NMAPS"),
Carl Schaffer's avatar
Carl Schaffer committed
177
            T_XEL=len(self.maps),
178
179
180
            EXPTIME=file.query("XPOSURE") / 1000.,
            ATMOS_R0_MIN=min(self.query("ATMOS_R0")),
            ATMOS_R0_MAX=max(self.query("ATMOS_R0")),
181
182
183
184
            HELIOPROJECTIVE_X_MIN=xmin,
            HELIOPROJECTIVE_X_MAX=xmax,
            HELIOPROJECTIVE_Y_MIN=ymin,
            HELIOPROJECTIVE_Y_MAX=ymax,
185
186
187
188
189
190
191
192
193
            MU=np.mean(self.query("mu")),
            THETA=np.mean(self.query("theta")),
            S_XEL1=file.query("NAXIS1"),
            S_XEL2=max(self.query("ISTEP")),
            S_RESOLUTION=file.query("CDELT1"),
            WAVELENGTH_MAX=max(self.query("AWAVMAX")) / 10.,
            WAVELENGTH_MIN=min(self.query("AWAVMIN")) / 10.,
            EM_XEL=file.query("NAXIS3"),
            EM_RESOLUTION=file.query("CDELT3"),
Carl Schaffer's avatar
Carl Schaffer committed
194
            POL_STATES="IQUV" if "NAXIS4" in file.header else "",
195
196
197
198
199
200
201
202
203
204
205
            POL_XEL=file.query("NAXIS4") if "NAXIS4" in file.header else 0,
            DIMENSIONS=file.query("NAXIS"),
            BTYPE=file.query("BTYPE"),
            TELESCOPE="GREGOR",
            INSTRUMENT="GRIS",
            TARGET=file.query("OBS_TRGT"),
            FILTER=get_filter(file.query("AWAVLNTH") / 10.),
            IOBS=file.query("IOBS"),
            STEPSIZE=file.query("STEPSIZE"),
            OBS_MODE=file.query("OBS_MODE")
        )
Carl Schaffer's avatar
Carl Schaffer committed
206

207
208
209
        # Checkthat all fields are contained in the observation description
        fields = gris_fields
        assert all([k in observation_description for k in fields.keys()])
Carl Schaffer's avatar
Carl Schaffer committed
210

211
        self.properties = observation_description
212
213

        # Set is parsed flag
Carl Schaffer's avatar
Carl Schaffer committed
214
215
        self.is_parsed = True

216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
    def query(self, attribute):
        """
        Get query attribute either from self or search constituents for that attribute
        Args:
            attribute: name of attribute to query

        Returns:
            res: value of given attribute if the current instance has the attribute, list of attributes for all
            constituents if there is no attribute in "self"
        """
        if hasattr(self, attribute):
            res = getattr(self, attribute)
        else:
            res = [f.query(attribute) for f in self.files]
        return res

232
    @property
233
234
235
236
237
238
    def verbose_name(self, date=None):
        """Make verbose name for GRIS runs, turns the date of the run,
        or a date passed as an argument into a string like 26aug14.001
        where the sectionafter the point is the runnumber of the
        corresponding run. If a specific date is passed, only the
        first section of the string is generated.
239

240
241
242
243
244
        Args:
            date (datetime): specific date defaults to None
        Returns:
            outstring (string): verbose date
        """
245
246
247
        if date is None:
            outstring = self.date.strftime("%d%b%y").lower()
            outstring += f".{int(self.runnumber):03d}"
248
        else:
249
            outstring = date.strftime("%d%b%y").lower()
250
251
252
253
254
        return outstring

    def matching_files(self, subfolder, pattern=None):
        query = self.verbose_name
        if not pattern:
255
            files = glob(join(self.path, subfolder, "*"))
256
257
258
259
260
        else:
            query += pattern
            files = reglob(join(self.path, subfolder), query)
        return files

261
262
263
264
265
    @property
    def obs_name(self):
        """Observation name unique within all data for the instrument"""
        return f"{self.date.strftime('%Y%m%d')}_{self.runnumber:03d}"

266
267
    @property
    def files_l0(self):
268
        return self.matching_files("level0", self._l0pattern)
269
270
271

    @property
    def files_l1(self):
272
        return self.matching_files("level1", self._l1pattern)
273
274
275

    @property
    def files_l1_split(self):
276
        files = glob(join(self.path, "level1_split", "*.fits"))
277
        res = list(
Carl Schaffer's avatar
Carl Schaffer committed
278
279
            filter(lambda x: gris_run_number(x) == self.runnumber, files)
        )
280
281
282
283
        return res

    @property
    def files_l1_split_mod(self):
284
        files = glob(join(self.path, "level1_split_modified", "*.fits"))
285
        res = list(
Carl Schaffer's avatar
Carl Schaffer committed
286
287
            filter(lambda x: gris_run_number(x) == self.runnumber, files)
        )
288
289
290
291
        return res

    @property
    def files_l2(self):
292
        files = glob(join(self.path, "level2", "*"))
293
294
295
296
297
        res = list(filter(lambda x: self.verbose_name in basename(x), files))
        return res

    @property
    def files_l3(self):
298
        files = glob(join(self.path, "level3", "*"))
299
300
301
302
303
        res = list(filter(lambda x: self.verbose_name in basename(x), files))
        return res

    @property
    def context_files(self):
304
        return self.matching_files("context_data")
Carl Schaffer's avatar
Carl Schaffer committed
305

Carl Schaffer's avatar
Carl Schaffer committed
306
307
308
309
    @property
    def context_folder(self):
        return join(self.path, 'context_data')

310
311
    @property
    def map_saves(self):
312
        res = list(filter(lambda x: x[-1] == "m", self.files_l1))
313
        return res
Carl Schaffer's avatar
Carl Schaffer committed
314

315
316
    @property
    def coord_saves(self):
317
        res = list(filter(lambda x: "coord_mu" in x, self.context_files))
318
319
320
321
        return res

    @property
    def mask_files(self):
Carl Schaffer's avatar
Carl Schaffer committed
322
        res = list(filter(lambda x: "pen" in x and ".sav" in x, self.files_l2))
323
324
        return res

Carl Schaffer's avatar
Carl Schaffer committed
325
326
    @property
    def map_files(self):
327
        res = list(filter(lambda x: "maps" in x, self.files_l2))
Carl Schaffer's avatar
Carl Schaffer committed
328
329
330
331
        return res

    @property
    def spec_files(self):
332
        res = list(filter(lambda x: "cor_spec" in x, self.files_l2))
Carl Schaffer's avatar
Carl Schaffer committed
333
334
335
336
        return res

    @property
    def stokes_files(self):
337
        res = list(filter(lambda x: "stok" in x, self.files_l2))
Carl Schaffer's avatar
Carl Schaffer committed
338
339
        return res

340
341
    @property
    def loc_previews(self):
342
        def filterfunc(x):
343
            return re.search(f"HMI_{self.verbose_name}\\.png|{self.runnumber}_location.png", x)
344

345
346
347
        res = list(filter(filterfunc, self.context_files))
        return res

348
349
350
351
352
    @property
    def spatial_uncertainties(self):
        """Average Uncertainties as described by WCS"""
        return np.mean(self.query("CSYER1")), np.mean(self.query("CSYER2"))

353
354
    @property
    def map_previews(self):
355
356
357
        def filterfunc(x):
            return re.match(f"{self.verbose_name}\\.png", basename(x))

358
359
360
361
362
        res = list(filter(filterfunc, self.context_files))
        return res

    @property
    def logfiles(self):
363
        files = glob(join(self.path, basename(self.path) + ".txt"))
364
365
        return files

366
367
368
369
370
371
    @property
    def observers(self):
        if self.logfiles:
            observers = get_observers(self.logfiles[0])
        return observers

372
373
    @property
    def was_aborted(self):
374
        # determine number of maps found
375
376
        map_lengths = self.map_lengths
        n_maps = len(map_lengths)
377
378
379
380
381
382
383
384

        if n_maps == 0:
            return None

        # check header for planned number of maps
        planned_maps = 0
        planned_steps = 0
        header = fitsio.getheader(self.files[0].path, ignore_missing_end=True)
385
386
        planned_maps = int(header["NMAPS"])
        planned_steps = int(header["NSTEPS"])
387
388
389
390
391
392
393

        if n_maps != planned_maps:
            return True

        len_first_map = map_lengths[sorted(map_lengths.keys())[0]]
        if len_first_map != planned_steps:
            return True
Carl Schaffer's avatar
Carl Schaffer committed
394
        was_aborted = any([v != len_first_map for v in map_lengths.values()])
395
396
        return was_aborted

397
    def get_map_lengths(self):
398
        """Check number of slit positions and check if this number is as
399
        specified by the header
400
401
402
403
404

        Args:

        Returns:

405
        """
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
406
        for mapnumber in self.maps:
407
            slit_numbers = [
408
409
                f.slitposnumber for f in self.files if f.mapnumber == mapnumber
            ]
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
410
            self.map_lengths[mapnumber] = max(slit_numbers)
411

Carl Schaffer's avatar
Carl Schaffer committed
412
    def fix_times(self):
413
        """Fix all time values for consistent format in database before
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
414
        generating a database entry
415

Carl Schaffer's avatar
Linting    
Carl Schaffer committed
416
         """
Carl Schaffer's avatar
Carl Schaffer committed
417
        # parse datetime objects from strings
418
        keys = ["UT-MIN", "UT-MAX"]
419
        patterns = ["%Y-%m-%dT%H:%M:%S.%f", "%Y-%m-%dT%H:%M:%S"]
Carl Schaffer's avatar
Carl Schaffer committed
420

421
422
        for key in keys:
            value = self.properties[key]
423
424
425
426
427
428
429
            for p in patterns:
                try:
                    time = datetime.datetime.strptime(value, p)
                    break
                except ValueError:
                    continue

Carl Schaffer's avatar
Linting    
Carl Schaffer committed
430
            self.properties[key] = time
Carl Schaffer's avatar
Carl Schaffer committed
431
432
433

    def get_native(self):
        """Extract runwise parameter from filelist"""
Carl Schaffer's avatar
Carl Schaffer committed
434
        # Check for a file that has been successfully parsed
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
435
436
        for gris_fitsfile in self.files:
            if gris_fitsfile.is_parsed is True:
Carl Schaffer's avatar
Carl Schaffer committed
437
438
                break
        else:
Carl Schaffer's avatar
Carl Schaffer committed
439
            self.errors.append(
Carl Schaffer's avatar
Carl Schaffer committed
440
441
                "ERROR in %s: No file successfully parsed!" % (str(self))
            )
Carl Schaffer's avatar
Carl Schaffer committed
442
443
            return -1

Carl Schaffer's avatar
Carl Schaffer committed
444
        # Try retrieving native keys from file, log error is unsuccessful
445
        first_file = self.files[0]
Carl Schaffer's avatar
Carl Schaffer committed
446

447
        for key, parameter in self.pars_native.items():
Carl Schaffer's avatar
Carl Schaffer committed
448
            try:
449
                self.properties[key] = first_file.query(parameter)
Carl Schaffer's avatar
Carl Schaffer committed
450
451
            except KeyError:
                message = "ERROR in %s: Key %s not available in fits file!"
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
452
                message = message % (str(self), parameter)
Carl Schaffer's avatar
Carl Schaffer committed
453
                self.errors.append(message)
Carl Schaffer's avatar
Carl Schaffer committed
454

455
456
457
458
459
        # get naxisi values
        for i in range(1, self.properties["NAXIS"] + 1):
            key = f"NAXIS{i}"
            self.properties[key] = first_file.query(key)

Carl Schaffer's avatar
Carl Schaffer committed
460
461
462
        # fix type and value for wavelength
        self.properties["WAVELENG"] = int(self.properties["WAVELENG"] / 10)

463
    def get_minmax(self, key, parameter):
Carl Schaffer's avatar
Carl Schaffer committed
464
465
        """get parameter range from fits files

466
        Args:
Carl Schaffer's avatar
Carl Schaffer committed
467
          key:
468
          parameter: key to process, valid values are any values that
469
470
471
472
473
474
475
        are included in  GrisFitsFile.properties.keys()
          newkey: key under which to store the values in self
        default is identical to original key for each
        original key a set of new keys appended with
        '-MIN' and '-MAX' are added

        Returns:
Carl Schaffer's avatar
Carl Schaffer committed
476
477

        """
478

Carl Schaffer's avatar
Carl Schaffer committed
479
        try:
480
            data = [
481
                gris_fitsfile.query(parameter) for gris_fitsfile in self.files
482
            ]
483

484
            # drop empty strings which are inserted instead of invalid values
485
            data = [d for d in data if d != ""]
486
487
488
489
490
491
492
493
494

            defaults = {
                "AZIMUT": -1.0,
                "ELEVATIO": -1,
                "SLITPOSX": 0,
                "SLITPOSY": 0,
                "FRIEDR0": 0,
            }

Carl Schaffer's avatar
Carl Schaffer committed
495
            if not data:
496
497
498
499
500
                print(f"No valid entries found in header for {key}")
                if key in defaults:
                    print(f"Defaulting to {defaults[key]}")
                    data.append(defaults[key])

501
            # get min and max
502
503
            self.properties[key + "-MAX"] = max(data)
            self.properties[key + "-MIN"] = min(data)
504

Carl Schaffer's avatar
Linting    
Carl Schaffer committed
505
        except KeyError:
Carl Schaffer's avatar
Carl Schaffer committed
506
            self.errors.append(
507
508
                "Missing Parameter: %s not found for %s\n"
                % (parameter, str(self))
Carl Schaffer's avatar
Carl Schaffer committed
509
            )
Carl Schaffer's avatar
Carl Schaffer committed
510
511
512
513
514

    def get_previews(self):
        """look for preview files and store their paths in
        self.previews as 'key':path

515
516
517
518
        Args:

        Returns:

Carl Schaffer's avatar
Carl Schaffer committed
519
520
        """
        run = self.runnumber
Carl Schaffer's avatar
Carl Schaffer committed
521

Carl Schaffer's avatar
Carl Schaffer committed
522
        # retrieve splits_level1_folder from first fits file
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
523
        folder = dirname(self.files[0].filename)
Carl Schaffer's avatar
Carl Schaffer committed
524
525

        # navigate to root folder:
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
526
        folder = abspath(join(folder, "../"))
527
        mapfolder = folder.replace("/gris/", "/gris_maps/")
Carl Schaffer's avatar
Carl Schaffer committed
528

Carl Schaffer's avatar
minor    
Carl Schaffer committed
529
        preview_map = []
530
531
        preview_loc = self.loc_previews
        preview_log = self.logfiles
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
532

Carl Schaffer's avatar
Carl Schaffer committed
533
        preview_map = glob(
534
            join(self.context_folder, f"???????.{run:03}.png")
Carl Schaffer's avatar
Carl Schaffer committed
535
        )
Carl Schaffer's avatar
Carl Schaffer committed
536
        if not preview_map:
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
537
            preview_map = glob(
Carl Schaffer's avatar
Carl Schaffer committed
538
539
                join(mapfolder, "context_data", f"*{run:03}_*.gif")
            )
Carl Schaffer's avatar
Carl Schaffer committed
540
        if not preview_map:
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
541
            preview_map = glob(
Carl Schaffer's avatar
Carl Schaffer committed
542
543
                join(mapfolder, "context_data", f"*{run:03}_???.png")
            )
Carl Schaffer's avatar
Carl Schaffer committed
544

545
        self.previews = {}
Carl Schaffer's avatar
Carl Schaffer committed
546
        if preview_loc:
547
            self.previews["PREVIEWLOC"] = preview_loc[0]
Carl Schaffer's avatar
Carl Schaffer committed
548
        if preview_log:
549
            self.previews["PREVIEWLOG"] = preview_log[0]
Carl Schaffer's avatar
Carl Schaffer committed
550
        if preview_map:
551
            self.previews["PREVIEWMAP"] = preview_map[0]
Carl Schaffer's avatar
Carl Schaffer committed
552

Carl Schaffer's avatar
Carl Schaffer committed
553
554
555
556
557
    def make_bson(self):
        """construct BSON entry for run. parameters RUNNUMBER and NFILES are
        always inserted. Apart from that, all parameters from the argument are
        appended

558
        Returns:
Carl Schaffer's avatar
Carl Schaffer committed
559
560

        """
561
        return self.properties.copy()
Carl Schaffer's avatar
Carl Schaffer committed
562

Carl Schaffer's avatar
Carl Schaffer committed
563
    def check_data_levels(self):
564
        """Check which data levels are present in parent directory
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
565

566
        todo : so far there is no checking wheter a specific run is
Carl Schaffer's avatar
Carl Schaffer committed
567
568
569
        contained in all folders, we assume that if the folder is
        present in the day, the corresponding run is contained in that
        day. this could be improved!
570
571
572
573
574

        Args:

        Returns:

Carl Schaffer's avatar
Carl Schaffer committed
575
        """
Carl Schaffer's avatar
Carl Schaffer committed
576
        # Keys under which to store the result
577
        keys = ["DATALVL0", "DATALVL1", "DATALVL2", "DATALV1S"]
Carl Schaffer's avatar
Carl Schaffer committed
578
        # folders to check for, order MUST match 'keys' above!
579
        folders = ["level0", "level1", "level2", "level1_split"]
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
580
581
582
        path = dirname(self.files[0].filename)
        path = abspath(path)
        directory = join(path, "..")
Carl Schaffer's avatar
Carl Schaffer committed
583
584

        # check if given folder is present and set run flags
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
585
586
587
        for k, folder in zip(keys, folders):
            path = join(directory, folder)
            folder_present = isdir(path)
Carl Schaffer's avatar
Carl Schaffer committed
588
            self.properties[k] = folder_present
Carl Schaffer's avatar
Carl Schaffer committed
589
590

    def __str__(self):
591
592
        outstring = f"GrisRun#{self.runnumber:3d}  Date: {self.date.strftime('%Y-%m-%d')},"
        outstring += f" {len(self.maps):3d} maps, {len(self.files):3d} split files"
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
593
        return outstring
Carl Schaffer's avatar
Carl Schaffer committed
594

Carl Schaffer's avatar
Carl Schaffer committed
595
596
597
598
    def __repr__(self):
        return str(self)

    def get_cube(self):
599
        """Retrieve data cube for run from slit files"""
Carl Schaffer's avatar
Carl Schaffer committed
600
        # Determine data shape and initialize numpy cube
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
601
        gris_fitsfile = fitsio.open(
Carl Schaffer's avatar
Carl Schaffer committed
602
603
            self.files[0].filename, ignore_missing_end=True
        )
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
604
605
606
607
        shape = gris_fitsfile[0].data.shape
        gris_fitsfile.close()
        n_x = len(self.files)
        n_y = shape[1]
Carl Schaffer's avatar
Carl Schaffer committed
608
609
        nlambda = shape[-1]
        npol = shape[0]
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
610
        cube = np.ndarray([npol, n_x, n_y, nlambda])
Carl Schaffer's avatar
Carl Schaffer committed
611
612

        # get data into cube:
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
613
        for i, fitsfile in enumerate(self.files):
Carl Schaffer's avatar
Carl Schaffer committed
614
615
616
            data = fitsio.open(fitsfile.filename, ignore_missing_end=True)[
                0
            ].data
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
617
            cube[:, i, :, :] = data
618
            data.close()
Carl Schaffer's avatar
Carl Schaffer committed
619
620
621
622
        self.cube = cube

    def get_mean_specs(self):
        """Calculate mean spectra for Stokes I,Q,U and V"""
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
623
        # check if cube has been made
624
        if not hasattr(self, "cube"):
Carl Schaffer's avatar
Carl Schaffer committed
625
626
627
            print("Missing cube, can't plot!")
            return

Carl Schaffer's avatar
Linting    
Carl Schaffer committed
628
        # average over spatial dimensions and slit positions
Carl Schaffer's avatar
Carl Schaffer committed
629
630
        specs = np.ndarray([self.cube.shape[0], self.cube.shape[-1]])
        for i in range(specs.shape[0]):
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
631
            specs[i, :] = self.cube[i, :, :, :].mean(axis=0).mean(axis=0)
Carl Schaffer's avatar
Carl Schaffer committed
632
633
634
        self.specs = specs

    def plot_toti(self):
635
636
        """ plot self.cube totI channel

Carl Schaffer's avatar
Linting    
Carl Schaffer committed
637
        """
638
639
        if not hasattr(self, "cube"):
            print("Missing cube, not plotting!")
Carl Schaffer's avatar
Carl Schaffer committed
640
            return
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
641
        toplot = self.cube[0, :, :, :].mean(axis=2)
Carl Schaffer's avatar
Carl Schaffer committed
642
643
        plt.imshow(
            np.rot90(toplot),
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
644
645
            vmax=np.percentile(toplot, 95),
            vmin=np.percentile(toplot, 10),
646
            cmap="gray",
Carl Schaffer's avatar
Carl Schaffer committed
647
        )
Carl Schaffer's avatar
Carl Schaffer committed
648
649
650
        plt.ion()
        plt.show()

651
652
    def plot_location(self):
        coords = self.calculate_bounding_box()
Carl Schaffer's avatar
Carl Schaffer committed
653
654
        date = min(self.query("obs_time"))
        fig, ax = make_loc_plot(coords[:2], coords[2:], date, uncertainties=self.spatial_uncertainties)
655
656
657
658
        plt.ion()
        plt.show()
        return fig,ax

Carl Schaffer's avatar
Carl Schaffer committed
659
    def plot_specs(self):
660
        """plot average spectra for given run. requres average spectra and cube
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
661
        """
Carl Schaffer's avatar
Carl Schaffer committed
662
        # check if mean spectra have been calculated
663
664
        if not hasattr(self, "specs"):
            print("Missing specs, not plotting!")
Carl Schaffer's avatar
Carl Schaffer committed
665
            return
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
666

Carl Schaffer's avatar
Carl Schaffer committed
667
        for i in range(4):
668
            plt.plot(self.specs[i, :] / np.linalg.norm(self.specs[i, :]))
Carl Schaffer's avatar
Carl Schaffer committed
669
670
        plt.ion()
        plt.show()
671

672
673
674
675
676
677
678
679
680
681
    def calculate_bounding_box(self):
        """
        Determine the bounding box for the run. The box is calculated as the area covered by the first map in the
        observation. If there are subsequent maps, they are ignored.

        Returns:
            xmin,xmax,ymin,ymax : Helioprojective coordinates in arcsec
        """
        first_map_id = self.maps[0]
        map_files = sorted(list(filter(lambda x: x.mapnumber == first_map_id, self.files)), key=lambda x: x.filename)
682
683
        first_coords = map_files[0]._coords_from_wcs
        last_coords = map_files[-1]._coords_from_wcs
684
685
686
687
688
689
        x_vals = np.concatenate([first_coords[:, 0], last_coords[:, 0]]).value
        y_vals = np.concatenate([first_coords[:, 1], last_coords[:, 1]]).value
        xmin, xmax = min(x_vals), max(x_vals)
        ymin, ymax = min(y_vals), max(y_vals)
        return xmin, xmax, ymin, ymax

690

Carl Schaffer's avatar
Carl Schaffer committed
691
# noinspection PyMethodOverriding
692
693
694
695
696
class EmptyGrisRun(GrisRun):
    def __init__(self, path):
        """run constructed from a path to a single l1 or l0 filename"""
        runnumber = gris_run_number(path)
        fn = basename(path)
697
        self.date = datetime.datetime.strptime(fn.split(".")[0], "%d%b%y")
698
699
700
        self.runnumber = runnumber
        self.map_lengths = {}
        self.files = []
701
        self.path = abspath(join(dirname(path), ".."))
702
703
        self.maps = []

Carl Schaffer's avatar
Carl Schaffer committed
704
    # noinspection PyMethodOverriding
705
706
    @property
    def verbose_name(self):
707
708
        outstring = self.date.strftime("%d%b%y").lower()
        outstring += f".{int(self.runnumber):03d}"
709
        return outstring
710
711
712
713
714


if __name__ == "__main__":
    files = glob("/dat/sdc/gris/20150426/level1_split/*_001_???_*.fits")
    gr = GrisRun(files)
715
    props = gr.parse()