#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Sep 15 20:10:33 2021

@author: Camaron.George
"""
import os
import glob
import numpy as np
import netCDF4 as nc
from scipy.interpolate import griddata

path = '/scratch2/NCEPDEV/ohd/Camaron.George/'
grid = 'EastGulf.gr3'
out = 'TideConstEastGulf.nc'

gridFile = path+'Scripts/Tide/'+grid
lon = []
lat = []
elems = []
bnodes = []
with open(gridFile) as f:
    next(f)
    line = f.readline()
    ne = int(line.split()[0])
    nn = int(line.split()[1])
    for i in range(nn):
        line = f.readline()
        lon.append(float(line.split()[1]))
        lat.append(float(line.split()[2]))
    for i in range(ne):
        f.readline()
    next(f)
    next(f)
    line = f.readline()
    nb = int(line.split()[0])
    for j in range(nb):
        bnodes.append(int(f.readline()))

lo = [lon[b-1] for b in bnodes]
la = [lat[b-1] for b in bnodes]

consfiles = glob.glob(path+'TidalConst/*')
consfiles.sort() #'k1','k2','m2','n2','o1','p1','q1','s2'
pha = np.zeros((len(lo),len(consfiles)))
amp = np.zeros((len(lo),len(consfiles)))

col = 0
for file in consfiles:
    data = nc.Dataset(file,'r')
    if col == 0:
        x = data.variables['lon'][:]
        y = data.variables['lat'][:]
        x,y = np.meshgrid(x,y)
        I = np.where(x > 180.0)
        x[I] = x[I]-360.0
        I = np.where((x < max(lo)+1) & (x > min(lo)-1) & (y < max(la)+1) & (y > min(la)-1))
        x = x[I]
        y = y[I]
        
    a = data.variables['amplitude'][:]
    a = a[I]
    p = data.variables['phase'][:]
    p = p[I]
    
    test = griddata((x,y),a,(lo,la),method='nearest')
    amp[:,col] = griddata((x,y),a,(lo,la),method='nearest')
    pha[:,col] = griddata((x,y),p,(lo,la),method='nearest')
    col+=1
    
if os.path.exists(path+'Scripts/Tide/'+out):
    os.remove(path+'Scripts/Tide/'+out)
    
# open a netCDF file to write
ncout = nc.Dataset(path+'Scripts/Tide/'+out,'w',format='NETCDF4')

# define axis size
ncout.createDimension('nOpenBndNodes',None)
ncout.createDimension('nConstituents',len(consfiles))

# create time axis
ncamp = ncout.createVariable('amp','f8',('nOpenBndNodes','nConstituents',))
ncpha = ncout.createVariable('pha','f8',('nOpenBndNodes','nConstituents',))

# copy axis from original dataset
ncamp[:] = amp
ncpha[:] = pha

ncout.close()