dougal-software/bin/p190.py

#!/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