# -*- coding: utf-8 -*-

"""

Created on Wed Apr 24 14:23:21 2024

@author: rein_pp

"""

import tl3_analysis_toolbox

#import sst_validation_2021

import xarray as xr

import matplotlib.pyplot as plt

import rioxarray

import rasterio

from rasterio.enums import Resampling

import pdb

import numpy as np

#from matplotlib.colors import LinearSegmentedColormap

from scipy.stats import linregress

#import mpl_scatter_density

import pandas as pd

from datetime import datetime

import rioxarray

from datetime import timedelta

import calendar

from sklearn.metrics import r2_score

import os

import warnings

import seaborn as sns

warnings.filterwarnings("ignore")

def buoy_date(row):

return datetime(int(row['Year']),int(row['Month']),int(row['Day']))

def plt_sst_buoy(data,var,ylabel,outfile):

fontsize=15

fig=plt.figure()

ax=fig.add_subplot(111)

ax.plot(data.sst_buoy,data[var],'.',markersize=1)

ax.plot([270,305],[270,305])

ax.set_xlim(270,305)

ax.set_ylim(270,305)

ax.set_aspect('equal')

ax.tick_params(labelsize=fontsize)

r2=r2_score(data.sst_buoy,data[var])

mad=np.mean(data.dif)

bias=np.mean(data[var]-data.sst_buoy)

n=len(data)

ax.text(290,280,'r2= '+str(round(r2,2)),fontsize=fontsize)

ax.text(290,277.5,'MAD= '+str(round(mad,2)),fontsize=fontsize)

ax.text(290,275,'Bias= '+str(round(bias,2)),fontsize=fontsize)

ax.text(290,272.5,'n= '+str(n),fontsize=fontsize)

ax.set_xlabel('NOAA Drifting Buoys SST [K]',fontsize=fontsize)

ax.set_ylabel(ylabel,fontsize=fontsize)

fig.tight_layout()

fig.savefig(outfile)

def buoy_validation_new():

aux=tl3_analysis_toolbox.aux()

prep=tl3_analysis_toolbox.l3_lst_sst_prep()

reproj=tl3_analysis_toolbox.reproj()

crop=tl3_analysis_toolbox.crop()

path='/nfs/IGARSS_2022/'

path_out='/nfs/IGARSS_2022/Validation_Revions_2024/buoy_outputs/'

path_buoy=path+'buoy_data/'

# Read in buoy data

print(datetime.now().strftime("%H:%M:%S")+' Read buoy data.')

data=pd.read_csv(path_buoy+'NOAA_Global_Lagrangian Drifter.csv')

data = data[(data.lat>26)&(data.lat<78)]

data = data[(data.lon>-41)&(data.lon<56)]

data=data[data.sea_surface_temperature<999]

data.sea_surface_temperature=data.sea_surface_temperature+273.15

data['time']=data.apply(buoy_date,axis=1)

for year in range(2007, 2014):

#print 'YEAR ',year

for month in range(1,13):

daynum=calendar.monthrange(year,month)[1]

for day in range(1,daynum+1):

print('Process '+str(datetime(year,month,day)))

tiles=['t01','t02','t03','t04']

for tile in tiles:

f_tl_dec=aux.get_l3_file_from_ida(str(year)+str(month).zfill(2)+str(day).zfill(2),

'Daily','SST','v01.01',tile)

if f_tl_dec is not None:

xds=xr.open_mfdataset(f_tl_dec,combine='nested',chunks={'x': 2300, 'y': 3250},concat_dim=['t'],

preprocess=prep.add_date)

data_val=pd.DataFrame()

data_dat=data[data.time==datetime(year,month,day)]

for buoy in data_dat.id.unique():

data_buoy=data_dat[data_dat.id==buoy]

lat_buoy=np.mean(data_buoy.lat)

lon_buoy=np.mean(data_buoy.lon)

lat_buoy,lon_buoy=reproj.coord_4326_to_3035(lat_buoy,lon_buoy)

stats_tl=crop.get_vars_by_coords(xds,['sst','tcwv','sat_zenith'],lat_buoy,lon_buoy)

lat_buoy=np.mean(data_buoy.lat)

lon_buoy=np.mean(data_buoy.lon)

sst_buoy=np.mean(data_buoy.sea_surface_temperature)

if np.isfinite(float(stats_tl['sst'].iloc[0])):

dat=pd.DataFrame({'sst_buoy':[sst_buoy],

'sst_tl':[float(stats_tl['sst'].iloc[0])],

'tcwv':[float(stats_tl['tcwv'].iloc[0])],

'sat_zenith':[float(stats_tl['sat_zenith'].iloc[0])],

'id':buoy,

'time':str(datetime(year,month,day))})

data_val=data_val.append(dat)

if len(data_val)>0:

# Filter some unrealistic TL SSTs

data_val=data_val[data_val.sst_tl<350]

# Write results to csv

outfile=path_out+str(datetime(year,month,day))[0:10]+'tl_buoy.csv'

data_val.to_csv(outfile)

def read_buoy_data(path_results):

files_buoy=[f for f in os.listdir(path_results) if 'buoy' in f]

dat_buoy=pd.DataFrame()

for file in files_buoy:

dat_buoy=dat_buoy.append(pd.read_csv(path_results+file))

dat_buoy['Date']=dat_buoy.time.apply(lambda x: datetime(int(x[0:4]),int(x[5:7]),int(x[8:10])))

dat_buoy['Month']=dat_buoy.Date.apply(lambda x: x.month)

dat_buoy['DOY']=dat_buoy['Date'].apply(lambda x: x.timetuple().tm_yday)

return dat_buoy

def filter_buoy_data(dat_buoy):

dat_buoy=dat_buoy[dat_buoy['sst_tl']>273]

dat_buoy['dif']=abs(dat_buoy['sst_tl']-dat_buoy['sst_buoy'])

data_new=pd.DataFrame()

for idi in dat_buoy.id.unique():

data_sub=dat_buoy[dat_buoy.id==idi]

data_sub['mean_dif']=abs(data_sub['sst_tl']-np.mean(data_sub['sst_tl']))

std=np.std(data_sub.sst_buoy)

#data_sub=data_sub[data_sub.dif<2*std]

data_sub=data_sub[data_sub.mean_dif<2*std]

data_new=data_new.append(data_sub)

return data_new

def confidence_plot_sst_dif(data,title,ylabel,outfile):

fontsize=24

fig=plt.figure()

fig.set_figwidth(20)

fig.set_figheight(6)

ax= fig.add_subplot(111)

# Median plot

p=sns.lineplot(data=data, x='date',y='q2',ax=ax,marker='o',color='blue',

legend='brief',label='Median',ms=1)

# confidence intervals

f1=ax.fill_between(x=data['date'],

y1=data['q1'],

y2=data['q3'], alpha=0.2,color='blue')

# legends

l1=ax.legend(loc='upper left', bbox_to_anchor=(0.8, 0.5, 0.5, 0.5),frameon=False,

fontsize=fontsize)

l2=ax.legend(handles=[f1], labels=['5-95% Quantile'],loc='upper left',

bbox_to_anchor=(0.8, 0.42, 0.5, 0.5),frameon=False,

fontsize=fontsize)

# axes

start=datetime(1990,1,1)

stop=datetime(2017,1,1)

ax.set_ylim(-5,5)

ax.set_xlim(start,stop)

ax.add_artist(l1)

ax.add_artist(l2)

ax.set_title(title,fontsize=fontsize)

ax.set_xlabel("Time",fontsize=fontsize)

ax.set_ylabel(ylabel,fontsize=fontsize)

plt.xticks(fontsize=fontsize)

plt.yticks(fontsize=fontsize)

# zero line

ax.plot([start,stop],[0,0],color='black')

mad=np.mean(data.mads)

md=np.mean(data.mds)

ax.text(datetime(1990,3,1),4,'MAD :'+str(round(mad,2)),fontsize=fontsize)

ax.text(datetime(1990,3,1),3,'Bias :'+str(round(md,2)),fontsize=fontsize)

plt.tight_layout()

fig.savefig(outfile)

def calc_quantiles(path,files,var):

q1=[]

q2=[]

q3=[]

date=[]

ns=[]

mds=[]

mads=[]

i=0

for file in files:

print(str(i)+'/'+str(len(files)))

dat=pd.read_csv(path+file)

dat=dat[(np.isfinite(dat.sst_dif))&(dat.sst_max>275)]

std=np.std(dat.sst_max)

dat['mean_dif']=abs(np.mean(dat['sst_max'])-dat['sst_max'])

dat=dat[dat['mean_dif']<2*std]

if len(dat)>0:

ns.append(len(dat))

q1.append(np.quantile(dat[var],0.05))

q2.append(np.quantile(dat[var],0.5))

q3.append(np.quantile(dat[var],0.95))

mds.append(np.mean(dat[var]))

mads.append(np.mean(abs(dat[var])))

#if np.quantile(dat[var],0.5)< -5:

#pdb.set_trace()

if file[0:6]=='Baltic':

date.append(datetime(int(file[11:15]),int(file[16:18]),int(file[19:21])))

if file[0:5]=='North':

date.append(datetime(int(file[10:14]),int(file[15:17]),int(file[18:20])))

i=i+1

data=pd.DataFrame({'date':date,'q1':q1,'q2':q2,'q3':q3,'ns':ns,'mds':mds,'mads':mads})

return data

if __name__ == '__main__':

#buoy_validation_new()

path_plot='E:/Publications/SST_analysis/Submission/ReviewRound1/Figures/Supplementary_Validation/'

path_res='E:/SST_Analysis/Results/'

path_val='E:/SST_Analysis/Validationfiles/'

# Plot buoy data

path_results='E:/Publications/SST_analysis/Submission/ReviewRound1/buoy_outputs/'

dat_buoy=read_buoy_data(path_results)

dat_buoy=filter_buoy_data(dat_buoy)

outfile=path_plot+'bouy_scatter_tl.png'

plt_sst_buoy(dat_buoy,'sst_tl','TIMELINE SST[K]',outfile)

pdb.set_trace()

# Plot CCI monthly max comparison

# Baltic Sea

res_file=path_res+'baltic_dif_quantiles.csv'

files=[file for file in os.listdir(path_val) if file[0:6]=='Baltic']

var='sst_dif'

#data=calc_quantiles(path_val,files,var)

#data.to_csv(res_file)

data=pd.read_csv(res_file)

data['date']=data.date.apply(lambda x:datetime(int(x[0:4]),int(x[5:7]),int(x[8:10])))

outfile=path_plot+'baltic_confidence_dif.png'

confidence_plot_sst_dif(data, 'Baltic Sea SST difference', 'TIMELINE - CCI SST[K]', outfile)

# North Sea

files=[file for file in os.listdir(path_val) if file[0:5]=='North']

res_file=path_res+'north_dif_quantiles.csv'

var='sst_dif'

#data=calc_quantiles(path_val,files,var)

#data.to_csv(res_file)

data=pd.read_csv(res_file)

data['date']=data.date.apply(lambda x:datetime(int(x[0:4]),int(x[5:7]),int(x[8:10])))

outfile=path_plot+'north_confidence_dif.png'

confidence_plot_sst_dif(data, 'North Sea SST difference', 'TIMELINE - CCI SST[K]', outfile)