GrisRun.py 22.4 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
24
from ..generic import gris_fields
from . import GrisFitsFile
from .util import get_observers
from ..util.util import reglob, gris_run_number, get_filter
25
from tqdm import tqdm
Carl Schaffer's avatar
Carl Schaffer committed
26

Carl Schaffer's avatar
Carl Schaffer committed
27

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

31
32
33
34
    Args:
      gris_fits_files: list of gris files to form one run

    Returns:
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
35
      ris_run: GrisRun instance
Carl Schaffer's avatar
Carl Schaffer committed
36
    """
Carl Schaffer's avatar
Carl Schaffer committed
37

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

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

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

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

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

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

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

        # check if inputs are already instances of GrisFitsFile, if
        # not, generate instances
97
        elif isinstance(gris_fits_files[0], str):
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
98
99
100
101
102
103
            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)
104

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

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

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

        Args:

        Returns:

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

        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
143
            if fits_file.errors:
Carl Schaffer's avatar
Carl Schaffer committed
144
                # Get errors from file,move file to broken list
145
                for error in fits_file.errors:
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
146
                    self.errors.append(error)
147
148
149
                self.files_broken.append(
                    self.files.pop(self.files.index(fits_file))
                )
Carl Schaffer's avatar
Carl Schaffer committed
150
        if not self.files:
Carl Schaffer's avatar
Carl Schaffer committed
151
            self.errors.append(
Carl Schaffer's avatar
Carl Schaffer committed
152
153
                "ERROR in %s: No fits files for run!" % (str(self))
            )
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
154
            return False
Carl Schaffer's avatar
Carl Schaffer committed
155

156
157
158
        # Extract observation properites from files
        pass

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

163
164
        xmin, xmax, ymin, ymax = self.calculate_bounding_box()

165
        observation_description = dict(
166
            OBS_NAME=self.obs_name,
167
168
169
170
            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
171
            DATE_BEG=min(self.query("obs_time")),
172
            DATE_END=max(self.query("obs_time")) + datetime.timedelta(seconds=file.query("XPOSURE") / 1000.),
Carl Schaffer's avatar
Carl Schaffer committed
173
174
            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
175
            T_XEL=len(self.maps),
176
177
178
            EXPTIME=file.query("XPOSURE") / 1000.,
            ATMOS_R0_MIN=min(self.query("ATMOS_R0")),
            ATMOS_R0_MAX=max(self.query("ATMOS_R0")),
179
180
181
182
            HELIOPROJECTIVE_X_MIN=xmin,
            HELIOPROJECTIVE_X_MAX=xmax,
            HELIOPROJECTIVE_Y_MIN=ymin,
            HELIOPROJECTIVE_Y_MAX=ymax,
183
184
185
186
187
188
189
190
191
            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
192
            POL_STATES="IQUV" if "NAXIS4" in file.header else "",
193
194
195
196
197
198
199
200
201
202
203
            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
204

205
206
207
        # 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
208

209
        self.properties = observation_description
210
211

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

214
215
216
217
218
219
220
221
222
223
    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"
        """
Carl Schaffer's avatar
Carl Schaffer committed
224
        print(attribute)
225
226
227
228
229
230
        if hasattr(self, attribute):
            res = getattr(self, attribute)
        else:
            res = [f.query(attribute) for f in self.files]
        return res

231
    @property
232
233
234
235
236
237
    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.
238

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

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

260
261
262
263
264
    @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}"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        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)
384
385
        planned_maps = int(header["NMAPS"])
        planned_steps = int(header["NSTEPS"])
386
387
388
389
390
391
392

        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
393
        was_aborted = any([v != len_first_map for v in map_lengths.values()])
394
395
        return was_aborted

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

        Args:

        Returns:

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

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

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

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

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

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

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

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

454
455
456
457
458
        # 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
459
460
461
        # fix type and value for wavelength
        self.properties["WAVELENG"] = int(self.properties["WAVELENG"] / 10)

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

465
        Args:
Carl Schaffer's avatar
Carl Schaffer committed
466
          key:
467
          parameter: key to process, valid values are any values that
468
469
470
471
472
473
474
        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
475
476

        """
477

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

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

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

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

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

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

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

514
515
516
517
        Args:

        Returns:

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

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

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

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

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

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

Carl Schaffer's avatar
Carl Schaffer committed
552
553
554
555
556
    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

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

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

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

565
        todo : so far there is no checking wheter a specific run is
Carl Schaffer's avatar
Carl Schaffer committed
566
567
568
        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!
569
570
571
572
573

        Args:

        Returns:

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

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

    def __str__(self):
590
591
        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
592
        return outstring
Carl Schaffer's avatar
Carl Schaffer committed
593

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

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

        # get data into cube:
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
612
        for i, fitsfile in enumerate(self.files):
Carl Schaffer's avatar
Carl Schaffer committed
613
614
615
            data = fitsio.open(fitsfile.filename, ignore_missing_end=True)[
                0
            ].data
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
616
            cube[:, i, :, :] = data
617
            data.close()
Carl Schaffer's avatar
Carl Schaffer committed
618
619
620
621
        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
622
        # check if cube has been made
623
        if not hasattr(self, "cube"):
Carl Schaffer's avatar
Carl Schaffer committed
624
625
626
            print("Missing cube, can't plot!")
            return

Carl Schaffer's avatar
Linting    
Carl Schaffer committed
627
        # average over spatial dimensions and slit positions
Carl Schaffer's avatar
Carl Schaffer committed
628
629
        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
630
            specs[i, :] = self.cube[i, :, :, :].mean(axis=0).mean(axis=0)
Carl Schaffer's avatar
Carl Schaffer committed
631
632
633
        self.specs = specs

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

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

    def plot_specs(self):
651
        """plot average spectra for given run. requres average spectra and cube
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
652
        """
Carl Schaffer's avatar
Carl Schaffer committed
653
        # check if mean spectra have been calculated
654
655
        if not hasattr(self, "specs"):
            print("Missing specs, not plotting!")
Carl Schaffer's avatar
Carl Schaffer committed
656
            return
Carl Schaffer's avatar
Linting    
Carl Schaffer committed
657

Carl Schaffer's avatar
Carl Schaffer committed
658
        for i in range(4):
659
            plt.plot(self.specs[i, :] / np.linalg.norm(self.specs[i, :]))
Carl Schaffer's avatar
Carl Schaffer committed
660
661
        plt.ion()
        plt.show()
662

663
664
665
666
667
668
669
670
671
672
    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)
673
674
        first_coords = map_files[0]._coords_from_wcs
        last_coords = map_files[-1]._coords_from_wcs
675
676
677
678
679
680
        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

681

Carl Schaffer's avatar
Carl Schaffer committed
682
# noinspection PyMethodOverriding
683
684
685
686
687
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)
688
        self.date = datetime.datetime.strptime(fn.split(".")[0], "%d%b%y")
689
690
691
        self.runnumber = runnumber
        self.map_lengths = {}
        self.files = []
692
        self.path = abspath(join(dirname(path), ".."))
693
694
        self.maps = []

Carl Schaffer's avatar
Carl Schaffer committed
695
    # noinspection PyMethodOverriding
696
697
    @property
    def verbose_name(self):
698
699
        outstring = self.date.strftime("%d%b%y").lower()
        outstring += f".{int(self.runnumber):03d}"
700
        return outstring
701
702
703
704
705


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