#!/usr/bin/python3 """ P1/90 parsing functions. """ import math import re from datetime import datetime, timedelta, timezone from parse_fwr import parse_fwr def parse_p190_header (string): """Parse a generic P1/90 header record. Returns a dictionary of fields. """ names = [ "record_type", "header_type", "header_type_modifier", "description", "data" ] record = parse_fwr(string, [1, 2, 2, 27, 48]) return dict(zip(names, record)) def parse_p190_type1 (string): """Parse a P1/90 Type 1 non receiver group record.""" names = [ "record_type", "line_name", "spare1", "vessel_id", "source_id", "tailbuoy_id", "point_number", "latitude", "longitude", "easting", "northing", "water_depth", "doy", "time", "spare2" ] record = parse_fwr(string, [1, 12, 3, 1, 1, 1, 6, 10, 11, 9, 9, 6, 3, 6, 1]) return dict(zip(names, record)) def parse_p190_rcv_group (string): """Parse a P1/90 Type 1 receiver group record.""" names = [ "record_type", "rcv_group_number1", "easting1", "northing1", "depth1", "rcv_group_number2", "easting2", "northing2", "depth2", "rcv_group_number3", "easting3", "northing3", "depth3", "streamer_id" ] record = parse_fwr(string, [1, 4, 9, 9, 4, 4, 9, 9, 4, 4, 9, 9, 4, 1]) return dict(zip(names, record)) def parse_line (string): type = string[0] if string[:3] == "EOF": return if type == "H": return parse_p190_header(string) elif type == "R": return parse_p190_rcv_group(string) else: return parse_p190_type1(string) def p190_type(type, records): return [ r for r in records if r["record_type"] == type ] def p190_header(code, records): return [ h for h in p190_type("H", records) if h["header_type"]+h["header_type_modifier"] == code ] def object_id(record): return "".join([ record[k] for k in ["record_type", "vessel_id", "source_id", "tailbuoy_id", "rcv_group_number1", "rcv_group_number2", "rcv_group_number3"] if k in record]) def normalise_record(record): for key in record: # These are floats if key in ("easting", "northing", "water_depth", "easting1", "northing1", "depth1", "easting2", "northing2", "depth2", "easting3", "northing3", "depth3"): try: record[key] = float(record[key]) except ValueError: pass # These are ints elif key in ("point_number", "rcv_group_number1", "rcv_group_number2", "rcv_group_number3"): try: record[key] = int(record[key]) except ValueError: pass # These are geographic coordinates in DMS, # which we convert to decimal and return as a float elif key in ("latitude", "longitude"): try: record[key] = dms(record[key]) except ValueError: pass # These are probably strings elif "strip" in dir(record[key]): record[key] = record[key].strip() return record def normalise(records): for record in records: normalise_record(record) return records def from_file(path, only_records=None, shot_range=None, with_objrefs=False): records = [] with open(path) as fd: cnt = 0 line = fd.readline() while line: cnt = cnt + 1 if line == "EOF": break record = parse_line(line) if record is not None: if only_records: if not record["record_type"] in only_records: continue if shot_range: shot = int(record["point_number"]) if not (shot >= min(shot_range) and shot <= max(shot_range)): continue if with_objrefs: record["objref"] = object_id(record) records.append(record) line = fd.readline() return records def apply_tstamps(recordset, tstamp=None, fix_bad_seconds=False): #print("tstamp", tstamp, type(tstamp)) if type(tstamp) is int: year = tstamp elif type(tstamp) is str: ts_rx = re.compile(tstamp) date_value = next(iter(p190_header("0200", recordset)), {"data": ""})["data"].strip() year = datetime.strptime(date_value, tstamp).year #print("d,y", date_value, year) elif tstamp is None: # Take a guess. If "doy" of the first record is less than the current Julian day, # use the present year, else go back one year. year = datetime.utcnow().year today = int(datetime.strftime(datetime.utcnow(), "%j")) doy = next((int(r["doy"]) in r for r in recordset if "doy" in r), today) if doy > today: year -= 1 else: raise ValueError("Cannot handle tstamp of type %s" % type(tstamp)) prev = dict() for record in recordset: if "doy" in record and "time" in record: while True: doy = record["doy"] time = record["time"].replace(" ", "0") # Fix for some dodgy systems if fix_bad_seconds and time.endswith("60"): time = time[:4]+"59" expr = "{} {} {} UTC".format(year, doy, time) ts = datetime.strptime(expr, "%Y %j %H%M%S %Z").replace(tzinfo=timezone.utc) if object_id(record) in prev and doy < prev[object_id(record)]: # We have gone through 31 Dec and "doy" rolled back to 1 year += 1 else: record["tstamp"] = ts prev[object_id(record)] = doy break return recordset def dms(value): # 591544.61N hemisphere = 1 if value[-1] in "NnEe" else -1 seconds = float(value[-6:-1]) minutes = int(value[-8:-6]) degrees = int(value[:-8]) return (degrees + minutes/60 + seconds/3600) * hemisphere def tod(record): if "tstamp" in record and record["tstamp"].timestamp: return record["tstamp"].timestamp() else: time = record["time"] d = int(record["doy"]) h = int(time[:2]) m = int(time[2:4]) s = float(time[4:]) return d*86400 + h*3600 + m*60 + s def duration(record0, record1): ts0 = tod(record0) ts1 = tod(record1) return timedelta(seconds=ts1-ts0) def distance(record0, record1): x0, y0 = float(record0["easting"]), float(record0["northing"]) x1, y1 = float(record1["easting"]), float(record1["northing"]) return math.sqrt((x1-x0)**2 + (y1-y0)**2) def azimuth(record0, record1): x0, y0 = float(record0["easting"]), float(record0["northing"]) x1, y1 = float(record1["easting"]), float(record1["northing"]) return math.degrees(math.atan2(x1-x0, y1-y0)) % 360 def speed(record0, record1, knots=False): scale = 3600/1852 if knots else 1 t0 = tod(record0) t1 = tod(record1) return (distance(record0, record1) / math.fabs(t1-t0)) * scale