#####################################################################################
#Very simple script that tests the rose plot capability. MJE. 20220504.
#
######################################################################################

#!/usr/bin/python
import os
import os.path
#from netCDF4 import Dataset
import numpy as np
import datetime
import csv
import math
import haversine
from windrose import WindroseAxes
import glob
import pygrib
from matplotlib import pyplot as plt
import subprocess

#######################LIST OF VARIABLES########################################################
MET_PATH     = '/opt/MET/METPlus4_0/bin/'
DATA_PATH    = '/export/hpc-lw-dtbdev5/merickson/code/python/lapenta/data/'
working_dir  = '/export/hpc-lw-dtbdev5/merickson/code/python/lapenta/work/'
validday     = 3
ero_category = 'HIGH'
beg_date     = datetime.datetime(2015,1,1,12,0,0)
end_date     = datetime.datetime(2021,5,31,12,0,0)
################################################################################################

os.environ["LD_LIBRARY_PATH"] ='/opt/MET/METPlus4_0/external_libs/lib'
os.environ["LIBS"] = '/opt/MET/METPlus4_0/external_libs'

lat_obs = [40,40,40,40]
lon_obs = [-75,-75,-75,-75]
lat_mod = [38,38,42,42]
lon_mod = [-77,-73,-77,-73]

bin = [0,1,2,3,4,5,6,7,8,9,10]

magnitude = []
angle     = []
if len(lat_obs) > 0 and len(lon_obs) > 0 and len(lat_mod) > 0 and len(lon_mod) > 0:
	for line in range(0,len(lat_obs)):

		#Latitude model displacement (model south is negative)
		if float(lat_mod[line]) > float(lat_obs[line]):
			ysign = 1
		else:
			ysign = -1
		#longitude model displacement (model west is negative)
		if float(lon_mod[line]) > float(lon_obs[line]):
			xsign = 1
		else:
			xsign = -1
		#Calculate x and y vector distance separately
		ydiff = ysign * (haversine.distance((float(lat_mod[line]), float(lon_mod[line])), \
			(float(lat_obs[line]), float(lon_mod[line]))))
		xdiff = xsign * (haversine.distance((float(lat_mod[line]) , float(lon_mod[line])), \
			(float(lat_mod[line]), float(lon_obs[line]))))

		#Calculate displacement in degrees
		ydiff_lat = float(lat_mod[line]) - float(lat_obs[line])
		xdiff_lon = float(lon_mod[line]) - float(lon_obs[line])

		#Taking x/y calculate magnitude and angle (0 to 360 degrees)
		magnitude = np.append(magnitude,np.sqrt(ydiff**2 + xdiff**2))
		angle     = np.append(angle,(math.atan2(xdiff,ydiff)/math.pi*180)+180)
			
#windrose plotting son southern plains
ax = WindroseAxes.from_ax()
ax.bar(angle[1:3], magnitude[1:3], bins=bin,blowto=True)
ax.set_legend()
plt.savefig(working_dir+'test.png',bbox_inches='tight',dpi=150)
plt.close()

subprocess.call('scp '+working_dir+'test.png '+'hpc@vm-lnx-rzdm06:/home/people/hpc/ftp/erickson/lapenta',shell=True)