#! /usr/bin/env python3
"""proj_ll2ra — project a GRD file from lon/lat into range/azimuth.

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

Usage:  proj_ll2ra trans.dat phase_ll.grd phase_ra.grd [-I32/64]

Inputs:
  trans.dat        — file generated by llt_grid2rat (r a topo lon lat)
  phase_ll.grd     — GRD file in lon/lat coordinates
  -I32/64          — optional output sampling (range/azimuth px); default -I32/64
Output:
  phase_ra.grd     — same data in radar (range/azimuth) coordinates
"""
import os
import sys
from gmtsar_lib import run


def proj_ll2ra():
    if len(sys.argv) < 4:
        sys.exit(
            "Usage: proj_ll2ra trans.dat phase_ll.grd phase_ra.grd [-I32/64]\n"
            "  trans.dat    - file generated by llt_grid2rat (r a topo lon lat)\n"
            "  phase_ll.grd - input GRD in lon/lat coordinates\n"
            "  phase_ra.grd - output GRD in radar coordinates\n"
            "  -I32/64      - optional output sampling (range/azimuth)"
        )
    trans, phase_ll, phase_ra = sys.argv[1], sys.argv[2], sys.argv[3]
    II = sys.argv[4] if len(sys.argv) >= 5 else "-I32/64"
    V = "" if os.path.isfile(os.path.expanduser("~/.quiet")) else "-V"

    run("gmt set FORMAT_GEO_OUT D")
    run(f"gmt grd2xyz {phase_ll} -s -bo3f -fg > llp")
    # Build (lon,lat,range) and (lon,lat,azimuth) tables from trans.dat
    run(f"gmt gmtconvert {trans} -o3,4,0 -bi5d -bo3f > llr")
    run(f"gmt gmtconvert {trans} -o3,4,1 -bi5d -bo3f > lla")
    # Interpolate to a regular grid; -I0.083... = 5 arc-min, -I0.000833... = 3 arc-sec
    run(f"gmt surface llr `gmt gmtinfo llp -I0.08333333333 -bi3f` "
        f"-bi3f -I.00083333333333 -T.50 -Gllr.grd {V}")
    run(f"gmt surface lla `gmt gmtinfo llp -I0.08333333333 -bi3f` "
        f"-bi3f -I.00083333333333 -T.50 -Glla.grd {V}")
    # Sample range/azimuth at each input lon/lat point
    run("gmt grdtrack llp -nl -Gllr.grd -bi3f -bo4f > llpr")
    run("gmt grdtrack llpr -nl -Glla.grd -bi4f -bo5f > llpra")
    # Reorder columns to (range, azimuth, phase) and grid them
    run("gmt gmtconvert llpra -bi5f -bo3f -o3,4,2 > rap")
    run(f"gmt xyz2grd rap `gmt gmtinfo rap -I32/64 -bi3f` {II} -r -G{phase_ra} -bi3f")
    run("rm -f ll* rap")


if __name__ == "__main__":
    proj_ll2ra()
