#####################################################################################
#Code to run mode on the previously created PP vs ERO netCDF file and then
#gather/calculate simple 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_obs = working_dir+'latlon_vday'+str(validday)+'_ERO'+ero_category+'_simpleobj_obs'+'.csv'
latlon_mod = working_dir+'latlon_vday'+str(validday)+'_ERO'+ero_category+'_simpleobj_mod'+'.csv'
print(latlon_obs)
print(latlon_mod)

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)

#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_obs)==True):
	erase_latlon_content = open(latlon_obs,'r+')
	erase_latlon_content.seek(0)
	erase_latlon_content.truncate()

if (os.path.exists(latlon_mod)==True):
        erase_latlon_content = open(latlon_mod,'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)
	
	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
	
			#Define hour string	
			#valid_hr = filename[124:126]
			valid_hr = filename[filename.find('vhr')+3:filename.find('vhr')+5]
	 
			file = open(mode_stuff_dir+'mode_240000L_'+yrmonday_tom+'_120000V_240000A_obj.txt','r')
			#file = open(mode_stuff_dir+'mode_240000L_20190527_120000V_240000A_obj.txt','r')
			lines = file.readlines()
			for line in lines[1:]:
				values = line.strip().split()
				object_id = values[22]
				object_cat = values[23][2:5]
				lat = values[26]
				lon = values[27]
				area = values[31]
				if 'F' in object_id and (not('NA') in lat) and (not('NA') in lon) and not 'CF' in object_id:
					lat_mod=values[26]
					lon_mod=values[27]
					date_mod=values[6]
					area_mod=values[31]
					object_cat_mod=values[23][2:5]
					if (np.any(object_cat_mod != '000')):
						matched_status_mod = float(True)	 
					if (np.any(object_cat_mod == '000')):
						matched_status_mod = float(False)
					#sanity check
					#if matched_status_mod == False:
					#	print('ERO not matched')
					#	print(object_id)
					with open(latlon_mod, mode='a') as lines:
						latlon_lines_mod = csv.writer(lines, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL)
						latlon_lines_mod.writerow([str(date_mod),valid_hr,str(int(matched_status_mod)),\
							str(lat_mod),str(lon_mod),str(area_mod)])
					file.close()
				
				elif 'O' in object_id and (not('NA') in lat) and (not('NA') in lon) and not 'CO' in object_id:
					date_obs=values[6]
					lat_obs=values[26]
					lon_obs=values[27]
					area_obs=values[31]
					object_cat_obs=values[23][2:5]
					if (np.any(object_cat_obs != '000')):
						matched_status_obs = float(True)
					if (np.any(object_cat_obs == '000')):
						matched_status_obs = float(False)
					#sanity check
					#if matched_status_obs == False:	
					#	print('observation not matched')
					#	print(object_id)
					with open(latlon_obs, mode='a') as lines: 
						latlon_lines_obs = csv.writer(lines, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL)
						latlon_lines_obs.writerow([str(date_obs),valid_hr,str(int(matched_status_obs)),\
							str(lat_obs),str(lon_obs),str(area_obs)])
					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')