#####################################################################################
#Code to run mode on the previously created PP vs ERO netCDF file and then
#gather/calculate piared object statistics and save as a csv file.
#
#-Run runMODE.py first to create the netCDF files (modify mode_config for MODE sensitivity studies)
#-Run createCSV.py to create a csv of paired objects
#-Run createCSV_simpleobj.py to create a csv of simple objects
#-Run plotting.py to create non-grid plots for paired objects
#-Run plotting_grid.py to create grid plots for paired objects
#-Run plotting_grid_simpleobj.py to create grid plots for simple objects (needs work)
#
#Created in the summer/fall of 2020 by CLC, editted by MJE on 20210601.
#Additional modifications by AJ throughout summer of 2021.
#Further documented by MJE. 20210921.
#
######################################################################################

#!/usr/bin/python
import os
import os.path
#from netCDF4 import Dataset
import numpy as np
import datetime
import csv
import math
import haversine
import glob
import pygrib

#######################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     = 1
ero_category = 'SLGT'
beg_date     = datetime.datetime(2016,6,1,12,0,0)                           #Start date
end_date     = datetime.datetime(2016,6,15,12,0,0)                           #End date
################################################################################################

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

#Define the output file
latlon = working_dir+'latlon_vday'+str(validday)+'_ERO'+ero_category+'.csv'
print(latlon)

mode_stuff_dir = working_dir+'MODEobjects/validday'+str(validday)+'_'+ero_category+'/'

#Set proper directories and files
DATA_PATH_EX = DATA_PATH+'ERO_verif_day'+str(validday)+'_ALL_noMRGL//'
config_file  = working_dir+'mode_config_'+str(ero_category)
print(DATA_PATH_EX)
print(working_dir)

#Remove any old CSV file
#os.system('rm -rf '+mode_stuff_dir+'mode*')

#Make necessary directories
try:
	os.mkdir(working_dir)
except:
	pass
try:
        os.mkdir(mode_stuff_dir)
except:
        pass
try:
	os.mkdir(working_dir+'MODEobjects')
except:
	pass
try:
        os.mkdir(working_dir+'figures')
except:
        pass

#Convert datetime to julian dates
beg_date_jul = pygrib.datetime_to_julian(beg_date)
end_date_jul = pygrib.datetime_to_julian(end_date)

#remove old contents in the old  latlon.csv file
if (os.path.exists(latlon)==True):
	erase_latlon_content = open(latlon,'r+')
	erase_latlon_content.seek(0)
	erase_latlon_content.truncate()

for dates in range(int(round(beg_date_jul)),int(round(end_date_jul))+1): #through the dates

	#Create datetime element for day being loaded
	curdate      = pygrib.julian_to_datetime(dates)
	tomdate      = pygrib.julian_to_datetime(dates+1)
	yrmonday_cur = '{:04d}'.format(curdate.year)+'{:02d}'.format(curdate.month)+'{:02d}'.format(curdate.day)
	yrmonday_tom = '{:04d}'.format(tomdate.year)+'{:02d}'.format(tomdate.month)+'{:02d}'.format(tomdate.day)

	print(DATA_PATH_EX+'grid_stat_PP_ALL_ERO_s'+yrmonday_cur+'*.nc.gz')
	
	for filename in glob.glob(DATA_PATH_EX+'grid_stat_PP_ALL_ERO_s'+yrmonday_cur+'*.nc.gz'):
		if ('vhr09' in filename): #Only grab 09 UTC for validday === 1

			print(filename)	
			#Define hour string	
			#valid_hr = filename[124:126]
			valid_hr = filename[filename.find('vhr')+3:filename.find('vhr')+5]
	 
			#Lat and Lon
			lat_mod=[]
			lat_obs=[]
			lon_mod=[]
			lon_obs=[]
			date_mod=[]

			area_mod=[]
			area_obs=[]
			
			#print(mode_stuff_dir+'/mode_240000L_'+yrmonday_tom+'_120000V_240000A_obj.txt')
			print(mode_stuff_dir+'mode_240000L_'+yrmonday_tom+'_120000V_240000A_obj.txt')	
			file = open(mode_stuff_dir+'mode_240000L_'+yrmonday_tom+'_120000V_240000A_obj.txt','r')
			lines = file.readlines()
			for line in lines:
				values = line.strip().split()
				object_id = values[22]
				lat = values[26]
				lon = values[27]
				area = values[31]
				if 'CF' in object_id and (not('NA') in lat) and (not('NA') in lon):
					lat_mod=np.append(lat_mod,values[26])
					lon_mod=np.append(lon_mod,values[27])
					date_mod=np.append(date_mod,values[6])
					area_mod=np.append(area_mod,values[31])
				elif 'CO' in object_id and (not('NA') in lat) and (not('NA') in lon):
					lat_obs=np.append(lat_obs,values[26])
					lon_obs=np.append(lon_obs,values[27])
					area_obs=np.append(area_obs,values[31])
					
			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)):

				#	print(mode_stuff_dir+'/mode_240000L_'+yrmonday_tom+'_120000V_240000A_obj.txt')

					#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.sqrt(ydiff**2 + xdiff**2)
					angle = (math.atan2(xdiff,ydiff)/math.pi*180)+180
			
					#print(angle)
					with open(latlon, mode='a') as lines: 
						latlon_lines = csv.writer(lines, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL)

						latlon_lines.writerow([str(date_mod[line]),valid_hr,str(lat_mod[line]),str(lon_mod[line]),str(area_mod[line]), \
							str(lat_obs[line]),str(lon_obs[line]),str(area_obs[line]),str(magnitude),str(angle),str(ydiff_lat),\
                                                        str(xdiff_lon)])
					file.close()        

###############EXTRA CODE########################################################
##This is an extra snippit of code that can be used to read in a netCDF file
#from netCDF4 import Dataset
#f        = Dataset(working_dir+'MODEobjects/mode_240000L_'+yrmonday_tom+'_120000V_240000A_obj.nc', "a", format="NETCDF4")
#temp     = f.variables['fcst_obj_id'][:]
#print(np.unique(temp))
#temp     = f.variables['fcst_clus_id'][:]
#print(np.unique(temp))
#f.close()

##This is an extra snippit that pushes the postscript to the web for viewing at: https://ftp.wpc.ncep.noaa.gov/erickson/lapenta/
#os.system('scp '+working_dir+'MODEobjects/mode_240000L_'+yrmonday_tom+'_120000V_240000A.ps hpc@vm-lnx-rzdm01:/home/people/hpc/ftp/erickson/lapenta')

##This is an extra snippit to calculate a rose plot for a sample and port it to the web for viewing
#from windrose import WindroseAxes
#import matplotlib.pyplot as plt
#ax = WindroseAxes.from_ax()
#ax.bar(np.array([angle]), np.array([magnitude]), bins=[0,10,20,50,75,100,125,150,200,400,500],blowto=True)
#ax.set_legend()
#plt.savefig(working_dir+'figures/test.png')
#plt.close()
#os.system('scp '+working_dir+'figures/test.png hpc@vm-lnx-rzdm01:/home/people/hpc/ftp/erickson/lapenta')