#! /usr/bin/env python3
"""proj_ll2ra_ascii — ASCII variant of proj_ll2ra.

Python port of csh proj_ll2ra_ascii.csh (D. Sandwell 2007, M. Wei 2008).

Usage:  proj_ll2ra_ascii trans.dat phase_ll.xyz phase_ra.xyz

Inputs:
  trans.dat        — file generated by llt_grid2rat (r a topo lon lat)
  phase_ll.xyz     — ASCII text file of (lon lat value) triples
Output:
  phase_ra.xyz     — ASCII text file of (range azimuth value) triples
"""
import os
import sys
from gmtsar_lib import run


def proj_ll2ra_ascii():
    if len(sys.argv) < 4:
        sys.exit(
            "Usage: proj_ll2ra_ascii trans.dat phase_ll.xyz phase_ra.xyz\n"
            "  trans.dat    - file generated by llt_grid2rat (r a topo lon lat)\n"
            "  phase_ll.xyz - input ASCII file of (lon lat value)\n"
            "  phase_ra.xyz - output ASCII file of (range azimuth value)"
        )
    trans, phase_ll, phase_ra = sys.argv[1], sys.argv[2], sys.argv[3]
    V = "" if os.path.isfile(os.path.expanduser("~/.quiet")) else "-V"

    run(f"gmt gmtconvert {trans} -o3,4,0 -bi5d -bo3f > llr")
    run(f"gmt gmtconvert {trans} -o3,4,1 -bi5d -bo3f > lla")
    # 0.01 deg ≈ 1.1 km region step; 0.005 deg interpolation
    run(f"gmt surface llr `gmt gmtinfo {phase_ll} -I0.01` "
        f"-bi3f -I0.005 -T.50 -Gllr.grd {V}")
    run(f"gmt surface lla `gmt gmtinfo {phase_ll} -I0.01` "
        f"-bi3f -I0.005 -T.50 -Glla.grd {V}")
    run(f"gmt grdtrack {phase_ll} -nl -Gllr.grd > llpr")
    run("gmt grdtrack llpr -nl -Glla.grd > llpra")
    # Output (range, azimuth, value) columns
    run(f"awk '{{print $4,$5,$3}}' < llpra > {phase_ra}")
    run("rm -f ll*")


if __name__ == "__main__":
    proj_ll2ra_ascii()
