meteortools.fileformats.ftpDetectInfo

  1#
  2# Load an FTPdetectInfo file - copied from RMS
  3#
  4# Copyright (C) 2018-2023 Mark McIntyre
  5
  6import os
  7import numpy as np
  8import configparser as crp
  9import json
 10import datetime
 11
 12try:
 13    from ..utils import date2JD, angleBetweenSphericalCoords
 14except Exception:
 15    from meteortools.utils import date2JD, angleBetweenSphericalCoords
 16
 17
 18def filterFTPforSpecificTime(ftpfile, dtstr):
 19    """ filter FTPdetect file for a specific event, by time, and copy it into a new file  
 20
 21    Arguments:  
 22        ftpfile - [string] full path to the ftpdetect file to filter  
 23        dtstr   - [string] date/time of required event in yyyymmdd_hhmmss format  
 24
 25    Returns:  
 26        tuple containing  
 27            full name of the new file containing just matching events  
 28            the number of matching events  
 29    """
 30    meteor_list = loadFTPDetectInfo(ftpfile, time_offsets=None, join_broken_meteors=True, locdata=None)
 31    refdt = datetime.datetime.strptime(dtstr, '%Y%m%d_%H%M%S')
 32    #print(refdt)
 33    new_met_list = []
 34    for met in meteor_list:
 35        #print(met.ff_name)
 36        dtpart = datetime.datetime.strptime(met.ff_name[10:25], '%Y%m%d_%H%M%S')
 37        tdiff = (refdt - dtpart).seconds
 38        #print(tdiff)
 39        if abs(tdiff) < 21:
 40            print('adding one entry')
 41            new_met_list.append(met)
 42    newname = writeNewFTPFile(ftpfile, new_met_list)
 43    return newname, len(new_met_list)
 44
 45
 46def writeNewFTPFile(srcname, metlist):
 47    """ creates a FTPDetect file from a list of MeteorObservation objects  
 48        
 49    Arguments: 
 50        srcname     - [string] full path to the original FTP file  
 51        metlist     - list of MeteorObservation objects  
 52
 53    Returns: 
 54        the full path to the created file 
 55
 56    """
 57    outdir, fname = os.path.split(srcname)
 58    newname = os.path.join(outdir, f'{fname}.old')
 59    try:
 60        os.rename(srcname, newname)
 61    except:
 62        pass
 63    if os.path.isfile(srcname):
 64        srcname = srcname[:-4] + '_new.txt'
 65    with open(srcname, 'w') as ftpf:
 66        _writeFTPHeader(ftpf, len(metlist), outdir, False)
 67        metno = 1
 68        ffname = ''
 69        for met in metlist:
 70            if ffname == met.ff_name:
 71                metno = metno + 1
 72            else:
 73                metno = 1
 74                ffname = met.ff_name
 75            _writeOneMeteor(ftpf, metno, met.station_id, met.time_data, len(met.frames), met.fps, met.frames, 
 76                np.degrees(met.ra_data), np.degrees(met.dec_data), 
 77                np.degrees(met.azim_data), np.degrees(met.elev_data),
 78                None, met.mag_data, False, met.x_data, met.y_data, met.ff_name)
 79    return srcname
 80
 81
 82def loadFTPDetectInfo(ftpdetectinfo_file_name, time_offsets=None,
 83        join_broken_meteors=True, locdata=None):
 84    """ Loads an FTPDEtect file into a list of MeteorObservation objects  
 85
 86    Arguments:  
 87        ftpdetectinfo_file_name: [str] Path to the FTPdetectinfo file.  
 88        stations: [dict] A dictionary where the keys are stations IDs, and values are lists of:  
 89            - latitude +N in radians  
 90            - longitude +E in radians  
 91            - height in meters  
 92        
 93
 94    Keyword arguments:  
 95        time_offsets: [dict] (key, value) pairs of (stations_id, time_offset) for every station. None by 
 96            default.
 97        join_broken_meteors: [bool] Join meteors broken across 2 FF files.
 98
 99
100    Return:  
101        meteor_list: [list] A list of MeteorObservation objects filled with data from the FTPdetectinfo file.  
102
103    """
104    stations={}
105    if locdata is None:
106        dirname, fname = os.path.split(ftpdetectinfo_file_name)
107        cfgfile = os.path.join(dirname, '.config')
108        cfg = crp.ConfigParser()
109        cfg.read(cfgfile)
110        try: 
111            lat = float(cfg['System']['latitude'].split()[0])
112            lon = float(cfg['System']['longitude'].split()[0])
113            height = float(cfg['System']['elevation'].split()[0])
114        except:
115            # try reading from platepars file
116            ppf = os.path.join(dirname, 'platepars_all_recalibrated.json')
117            if not os.path.isfile(ppf):
118                return []
119            js = json.load(open(ppf, 'r'))
120            if len(js) < 10:
121                return []
122            lat = js[list(js.keys())[0]]['lat']
123            lon = js[list(js.keys())[0]]['lon']
124            height = js[list(js.keys())[0]]['elev']
125        statid= fname.split('_')[1]
126    else:
127        statid = locdata['station_code']
128        lat = float(locdata['lat'])
129        lon = float(locdata['lon'])
130        height = float(locdata['elev'])
131
132    stations[statid] = [np.radians(lat), np.radians(lon), height*1000]
133    meteor_list = []
134
135    with open(ftpdetectinfo_file_name) as f:
136        # Skip the header
137        for i in range(11):
138            next(f)
139
140        current_meteor = None
141
142        bin_name = False
143        cal_name = False
144        meteor_header = False
145
146        for line in f:
147            # Skip comments
148            if line.startswith("#"):
149                continue
150
151            line = line.replace('\n', '').replace('\r', '')
152
153            # Skip the line if it is empty
154            if not line:
155                continue
156
157            if '-----' in line:
158                # Mark that the next line is the bin name
159                bin_name = True
160
161                # If the separator is read in, save the current meteor
162                if current_meteor is not None:
163                    current_meteor._finish()
164                    meteor_list.append(current_meteor)
165                continue
166
167            if bin_name:
168                bin_name = False
169
170                # Mark that the next line is the calibration file name
171                cal_name = True
172
173                # Save the name of the FF file
174                ff_name = line
175
176                # Extract the reference time from the FF bin file name
177                line = line.split('_')
178
179                # Count the number of string segments, and determine if it the old or new CAMS format
180                if len(line) == 6:
181                    sc = 1
182                else:
183                    sc = 0
184
185                ff_date = line[1 + sc]
186                ff_time = line[2 + sc]
187                milliseconds = line[3 + sc]
188
189                year = ff_date[:4]
190                month = ff_date[4:6]
191                day = ff_date[6:8]
192
193                hour = ff_time[:2]
194                minute = ff_time[2:4]
195                seconds = ff_time[4:6]
196
197                year, month, day, hour, minute, seconds, milliseconds = map(int, [year, month, day, hour, 
198                    minute, seconds, milliseconds])
199
200                # Calculate the reference JD time
201                jdt_ref = date2JD(year, month, day, hour, minute, seconds, milliseconds)
202                continue
203
204            if cal_name:
205                cal_name = False
206                # Mark that the next line is the meteor header
207                meteor_header = True
208                continue
209
210            if meteor_header:
211                meteor_header = False
212                line = line.split()
213
214                # Get the station ID and the FPS from the meteor header
215                station_id = line[0].strip()
216                fps = float(line[3])
217
218                # Try converting station ID to integer
219                try:
220                    station_id = int(station_id)
221                except:
222                    pass
223
224                # If the time offsets were given, apply the correction to the JD
225                if time_offsets is not None:
226                    if station_id in time_offsets:
227                        print('Applying time offset for station {:s} of {:.2f} s'.format(str(station_id),
228                            time_offsets[station_id]))
229
230                        jdt_ref += time_offsets[station_id]/86400.0
231                    else:
232                        print('Time offset for given station not found!')
233
234                # Get the station data
235                if station_id in stations:
236                    lat, lon, height = stations[station_id]
237                else:
238                    print('ERROR! No info for station ', station_id, ' found in CameraSites.txt file!')
239                    print('Exiting...')
240                    break
241                # Init a new meteor observation
242                current_meteor = MeteorObservation(jdt_ref, station_id, lat, lon, height, fps,
243                    ff_name=ff_name)
244                continue
245
246            # Read in the meteor observation point
247            if (current_meteor is not None) and (not bin_name) and (not cal_name) and (not meteor_header):
248
249                line = line.replace('\n', '').split()
250
251                # Read in the meteor frame, RA and Dec
252                frame_n = float(line[0])
253                x = float(line[1])
254                y = float(line[2])
255                ra = float(line[3])
256                dec = float(line[4])
257                azim = float(line[5])
258                elev = float(line[6])
259
260                # Read the visual magnitude, if present
261                if len(line) > 8:
262                    mag = line[8]
263                    if mag == 'inf':
264                        mag = None
265                    else:
266                        mag = float(mag)
267                else:
268                    mag = None
269
270                # Add the measurement point to the current meteor 
271                current_meteor._addPoint(frame_n, x, y, azim, elev, ra, dec, mag)
272
273        # Add the last meteor the the meteor list
274        if current_meteor is not None:
275            current_meteor._finish()
276            meteor_list.append(current_meteor)
277
278    # Concatenate observations across different FF files ###
279    if join_broken_meteors:
280
281        # Go through all meteors and compare the next observation
282        merged_indices = []
283        for i in range(len(meteor_list)):
284
285            # If the next observation was merged, skip it
286            if (i + 1) in merged_indices:
287                continue
288
289            # Get the current meteor observation
290            met1 = meteor_list[i]
291
292            if i >= (len(meteor_list) - 1):
293                break
294
295            # Get the next meteor observation
296            met2 = meteor_list[i + 1]
297            
298            # Compare only same station observations
299            if met1.station_id != met2.station_id:
300                continue
301
302            # Extract frame number
303            met1_frame_no = int(met1.ff_name.split("_")[-1].split('.')[0])
304            met2_frame_no = int(met2.ff_name.split("_")[-1].split('.')[0])
305
306            # Skip if the next FF is not exactly 256 frames later
307            if met2_frame_no != (met1_frame_no + 256):
308                continue
309
310
311            # Check for frame continouty
312            if (met1.frames[-1] < 254) or (met2.frames[0] > 2):
313                continue
314
315            # Check if the next frame is close to the predicted position 
316
317            # Compute angular distance between the last 2 points on the first FF
318            ang_dist = angleBetweenSphericalCoords(met1.dec_data[-2], met1.ra_data[-2], met1.dec_data[-1],
319                met1.ra_data[-1])
320
321            # Compute frame difference between the last frame on the 1st FF and the first frame on the 2nd FF
322            df = met2.frames[0] + (256 - met1.frames[-1])
323
324            # Skip the pair if the angular distance between the last and first frames is 2x larger than the 
325            #   frame difference times the expected separation
326            ang_dist_between = angleBetweenSphericalCoords(met1.dec_data[-1], met1.ra_data[-1],
327                met2.dec_data[0], met2.ra_data[0])
328
329            if ang_dist_between > 2*df*ang_dist:
330                continue
331
332            # If all checks have passed, merge observations ###
333            # Recompute the frames
334            frames = 256.0 + met2.frames
335
336            # Recompute the time data
337            time_data = frames/met1.fps
338
339            # Add the observations to first meteor object
340            met1.frames = np.append(met1.frames, frames)
341            met1.time_data = np.append(met1.time_data, time_data)
342            met1.x_data = np.append(met1.x_data, met2.x_data)
343            met1.y_data = np.append(met1.y_data, met2.y_data)
344            met1.azim_data = np.append(met1.azim_data, met2.azim_data)
345            met1.elev_data = np.append(met1.elev_data, met2.elev_data)
346            met1.ra_data = np.append(met1.ra_data, met2.ra_data)
347            met1.dec_data = np.append(met1.dec_data, met2.dec_data)
348            met1.mag_data = np.append(met1.mag_data, met2.mag_data)
349
350            # Merge the FF file name and create a list
351            if (met1.ff_name is not None) and (met2.ff_name is not None):
352                met1.ff_name = met1.ff_name + ',' + met2.ff_name
353
354            # Sort all observations by time
355            met1._finish()
356
357            # Indicate that the next observation is to be skipped
358            merged_indices.append(i + 1)
359
360        # Removed merged meteors from the list
361        meteor_list = [element for i, element in enumerate(meteor_list) if i not in merged_indices]
362
363    return meteor_list
364
365
366class MeteorObservation(object):
367    """ Container for meteor observations.  
368    """
369    def __init__(self, jdt_ref, station_id, latitude, longitude, height, fps, ff_name=None):
370        """ Construct the MeteorObservation object.  
371        Arguments:  
372            jdt_ref: [float] Reference Julian date when the relative time is t = 0s.  
373            station_id: [str] Station ID.  
374            latitude: [float] Latitude +N in radians.  
375            longitude: [float] Longitude +E in radians.  
376            height: [float] Elevation above sea level (MSL) in meters. 
377            fps: [float] Frames per second. 
378
379        Keyword arguments:  
380            ff_name: [str] Name of the originating FF file. 
381        """
382        self.jdt_ref = jdt_ref
383        self.station_id = station_id
384        self.latitude = latitude
385        self.longitude = longitude
386        self.height = height
387        self.fps = fps
388        self.ff_name = ff_name
389        self.frames = []
390        self.time_data = []
391        self.x_data = []
392        self.y_data = []
393        self.azim_data = []
394        self.elev_data = []
395        self.ra_data = []
396        self.dec_data = []
397        self.mag_data = []
398        self.abs_mag_data = []
399
400    def _addPoint(self, frame_n, x, y, azim, elev, ra, dec, mag):
401        """ Adds the measurement point to the meteor.
402
403        Arguments:
404            frame_n: [flaot] Frame number from the reference time.
405            x: [float] X image coordinate.
406            y: [float] X image coordinate.
407            azim: [float] Azimuth, J2000 in degrees.
408            elev: [float] Elevation angle, J2000 in degrees.
409            ra: [float] Right ascension, J2000 in degrees.
410            dec: [float] Declination, J2000 in degrees.
411            mag: [float] Visual magnitude.
412        """
413        self.frames.append(frame_n)
414
415        # Calculate the time in seconds w.r.t. to the reference JD
416        point_time = float(frame_n)/self.fps
417        self.time_data.append(point_time)
418        self.x_data.append(x)
419        self.y_data.append(y)
420
421        # Angular coordinates converted to radians
422        self.azim_data.append(np.radians(azim))
423        self.elev_data.append(np.radians(elev))
424        self.ra_data.append(np.radians(ra))
425        self.dec_data.append(np.radians(dec))
426        self.mag_data.append(mag)
427
428    def _finish(self):
429        """ When the initialization is done, convert data lists to numpy arrays. """
430        self.frames = np.array(self.frames)
431        self.time_data = np.array(self.time_data)
432        self.x_data = np.array(self.x_data)
433        self.y_data = np.array(self.y_data)
434        self.azim_data = np.array(self.azim_data)
435        self.elev_data = np.array(self.elev_data)
436        self.ra_data = np.array(self.ra_data)
437        self.dec_data = np.array(self.dec_data)
438        self.mag_data = np.array(self.mag_data)
439        # Sort by frame
440        temp_arr = np.c_[self.frames, self.time_data, self.x_data, self.y_data, self.azim_data, 
441        self.elev_data, self.ra_data, self.dec_data, self.mag_data]
442        temp_arr = temp_arr[np.argsort(temp_arr[:, 0])]
443        self.frames, self.time_data, self.x_data, self.y_data, self.azim_data, self.elev_data, self.ra_data, \
444            self.dec_data, self.mag_data = temp_arr.T
445
446
447def _writeFTPHeader(ftpf, metcount, fldr, ufo=True):
448    """
449    Internal function to create the header of the FTPDetect file  
450    """
451    l1 = 'Meteor Count = {:06d}\n'.format(metcount)
452    ftpf.write(l1)
453    ftpf.write('-----------------------------------------------------\n')
454    if ufo is True:
455        ftpf.write('Processed with UFOAnalyser\n')
456    else:
457        ftpf.write('Processed with RMS 1.0\n')
458    ftpf.write('-----------------------------------------------------\n')
459    l1 = 'FF  folder = {:s}\n'.format(fldr)
460    ftpf.write(l1)
461    l1 = 'CAL folder = {:s}\n'.format(fldr)
462    ftpf.write(l1)
463    ftpf.write('-----------------------------------------------------\n')
464    ftpf.write('FF  file processed\n')
465    ftpf.write('CAL file processed\n')
466    ftpf.write('Cam# Meteor# #Segments fps hnr mle bin Pix/fm Rho Phi\n')
467    ftpf.write('Per segment:  Frame# Col Row RA Dec Azim Elev Inten Mag\n')
468
469
470def _writeOneMeteor(ftpf, metno, sta, evttime, fcount, fps, fno, ra, dec, az, alt, b, mag, 
471        ufo=True, x=None, y=None, ffname = None):
472    """   Internal function to write one meteor event into the file in FTPDetectInfo style
473    """
474    ftpf.write('-------------------------------------------------------\n')
475    if ffname is None:
476        ms = '{:03d}'.format(int(evttime.microsecond / 1000))
477
478        fname = 'FF_' + sta + '_' + evttime.strftime('%Y%m%d_%H%M%S_') + ms + '_0000000.fits\n'
479    else:
480        fname = ffname + '\n'
481    ftpf.write(fname)
482
483    if ufo is True:
484        ftpf.write('UFO UKMON DATA Recalibrated on: ')
485    else:
486        ftpf.write('RMS data reprocessed on: ')
487    ftpf.write(datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S.%f UTC\n'))
488    li = f'{sta} {metno:04d} {fcount:04d} {fps:04.2f} 000.0 000.0  00.0 000.0 0000.0 0000.0\n'
489    ftpf.write(li)
490
491    for i in range(len(fno)):
492        #    204.4909 0422.57 0353.46 262.3574 +16.6355 267.7148 +23.0996 000120 3.41
493        bri = 0
494        if b is not None:
495            bri = int(b[i])
496        if ufo is True:
497            # UFO is timestamped as at the first detection
498            thisfn = fno[i] - fno[0]
499            thisx = 0
500            thisy = 0
501        else:
502            thisfn = fno[i]
503            thisx = x[i]
504            thisy = y[i]
505        li = f'{thisfn:.4f} {thisx:07.2f} {thisy:07.2f} '
506        li = li + f'{ra[i]:8.4f} {dec[i]:+7.4f} {az[i]:8.4f} '
507        li = li + f'{alt[i]:+7.4f} {bri:06d} {mag[i]:.02f}\n'
508        ftpf.write(li)