#!/usr/bin/env python3
"""
python drop-in replacement for the IDL remap.pro routine

Copyright (C) 2017-2018 University of Wisconsin Space Science and Engineering Center

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""


import sys
import numpy as np
import math

import subprocess

image_file = 'image.dat'
lat_file = 'lat.dat'
lon_file = 'lon.dat'

_, image_type, ncol, nrow, rowsperscan, res, missing, ul_lon, lr_lon, lr_lat, ul_lat = sys.argv
ncol=int(ncol)
nrow=int(nrow)
rowsperscan=int(rowsperscan)
res=float(res)
missing=float(missing)
ul_lon=float(ul_lon)
lr_lon=float(lr_lon)
lr_lat=float(lr_lat)
ul_lat=float(ul_lat)

print(ncol, nrow, rowsperscan, res, missing, ul_lon, lr_lon, lr_lat, ul_lat)

tag = 'image_proj'
xsize = 0
ysize = 0

lat = np.fromfile(lat_file)
lon = np.fromfile(lon_file)

if ul_lon >= 0.0 and lr_lon <= 0.0:
    lr_lon = lr_lon + 360.0

# set center of projection
c_lat = (lr_lat + ul_lat) * 0.5
c_lon = (lr_lon + ul_lon) * 0.5

# Create Map Projection Parameter (MPP) file
mpp_file = 'remap.mpp'
proj_name = 'Cylindrical Equidistant'
fmpp = open(mpp_file, 'w')
fmpp.write('Map Projection: %s\n' % proj_name)
fmpp.write('Map Reference Latitude: %.5f\n' % c_lat)
fmpp.write('Map Reference Longitude: %.5f\n' % c_lon)
fmpp.write('Map lr lat: %.5f\n' % lr_lat)
fmpp.write('Map lr lon: %.5f\n' % lr_lon)
fmpp.write('Map ul lat: %.5f\n' % ul_lat)
fmpp.write('Map ul lon: %.5f\n' % ul_lon)
fmpp.close()

# Compute XY coordinates for upper left corner of projection
ll2xy_ul_cmd = "echo %s %s | ll2xy %s" % (str(ul_lat), str(ul_lon), mpp_file)
print(ll2xy_ul_cmd)
r = subprocess.Popen(ll2xy_ul_cmd, stdout=subprocess.PIPE, shell=True)
if not r:
    print("echo", str(ul_lat), str(ul_lon), "|", "ll2xy", mpp_file)
    raise RuntimeError("Error: Non-zero exit status from ll2xy in cspp_remap")
r.wait()
for line in r.stdout:
    zlat, zlon, ul_x, ul_y, status = [float(s) for s in line.split()]

# Compute XY coordinates for lower right corner of projection
ll2xy_lr_cmd = "echo %s %s | ll2xy %s" % (str(lr_lat), str(lr_lon), mpp_file)
print(ll2xy_lr_cmd)
r = subprocess.Popen(ll2xy_lr_cmd, stdout=subprocess.PIPE, shell=True)
if not r:
    raise RuntimeError("Error: Non-zero exit status from ll2xy in cspp_remap")
r.wait()
for line in r.stdout:
    zlat, zlon, lr_x, lr_y, status = [float(s) for s in line.split()]

# Compute width and height of projection grid in meters
width_x = lr_x - ul_x
height_y = ul_y - lr_y

# Compute size of projection grid in pixels
xsize = math.ceil(width_x / res)
ysize = math.ceil(height_y / res)

# Compute map origin (center of grid)
xcen = (xsize * 0.5) - 0.5
ycen = (ysize * 0.5) - 0.5

# Create Grid Parameter Definition (GPD) file
gpd_file = 'remap.gpd'
fgpd = open(gpd_file, 'w')
fgpd.write('Grid MPP File: %s\n' % mpp_file)
fgpd.write("Grid Map Units per Cell: %.3f\n" % res)
fgpd.write("Grid Width: %.1f\n" % xsize)
fgpd.write("Grid Height: %.1f\n" % ysize)
fgpd.write("Grid Map Origin Column: %.1f\n" % xcen)
fgpd.write("Grid Map Origin Row: %.1f\n" % ycen)
fgpd.close()

# Run ll2cr (convert latitude-longitude pairs to column-row pairs)
nscan = nrow / rowsperscan

ll2cr_cmd = "ll2cr -v -f %5d %5d %5d %s %s %s %s" % (ncol, nscan, rowsperscan, lat_file, lon_file, gpd_file, tag)
print(ll2cr_cmd)
r = subprocess.Popen(ll2cr_cmd, stdout=subprocess.PIPE, shell=True)
r.wait()
if not r:
    raise RuntimeError("Error: Non-zero exit status from ll2cr in cspp_remap")

scan_first = 0
col_file = "%s_%s_%05d_%05d_%05d_%02d.img" % (tag, "cols", ncol, nscan, scan_first, rowsperscan)
row_file = "%s_%s_%05d_%05d_%05d_%02d.img" % (tag, "rows", ncol, nscan, scan_first, rowsperscan)

print("col_file:", col_file)
print("row_file:", row_file)

fornav_cmd = "fornav 1 -v -t %s -f %8d -F 0 -D 40 %5d %5d %5d %s %s %s %5d %5d %s" % (image_type, missing, ncol, nscan, rowsperscan, col_file, row_file, image_file, xsize, ysize, tag + ".dat")
print(fornav_cmd)
r = subprocess.Popen(fornav_cmd, stdout=subprocess.PIPE, shell=True)
r.wait()
if not r:
    raise RuntimeError("Error: Non-zero exit status from fornav in cspp_remap")