# Plot sea ice area (OISST)
#
#
#
# Script version 1.1
#
# Python version: 3.10
#
#
# Ver1.0: Created by Mariko Koseki, 18.06.2025
# ver1.1: updated by Mariko, 09.01.2026




import numpy as np
import pandas as pd
import datetime
from datetime import timedelta
from dateutil.relativedelta import relativedelta
import statistics
import math
import collections
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import matplotlib.ticker as ptick
import os
import glob
import sys
import platform
import calendar


print('--- Python version ---')
print(platform.python_version())

print('--- Package version ---')
print('numpy: ', np.__version__)
print('pandas: ', pd.__version__)
print('matplotlib: ', mpl.__version__)




##--- Set path ------------------------------------------##
## Input ##
'''
Set path to input directory ###
'''
seaice_dir = '/nird/datapeak/NS9873K/norcpm/validation/seaicepre/' #ver1.1
### OISST ###
out_oisst = seaice_dir + 'oisst/merge/'


## Output ##
'''
Set path to output file
'''
fig_dir = '/nird/datapeak/NS9873K/www/users/mariko/fig/seaice/OISST/' #ver1.1
out_dir = '/nird/datapeak/NS9873K/norcpm/validation/seaicepre/oisst/all/' #ver1.1

if not os.path.exists(fig_dir):
    os.makedirs(fig_dir)

if not os.path.exists(out_dir):
    os.makedirs(out_dir)




##--- Read CSV files ---------------------------------------##
## search for csv files ##
file_oisst = glob.glob(out_oisst + '*.csv')
file_oisst_sorted = sorted(file_oisst)
print('')
print(file_oisst_sorted)
print('')




df_all = pd.DataFrame() # initialize an empty dataframe
for filename in file_oisst_sorted:
    df = pd.read_csv(filename)

    
    ## Change 'object' to 'datetime' ##
    df['dtime'] = pd.to_datetime(df['dtime'])

    yyyy = os.path.splitext(os.path.basename(filename))[0][19:23]

    if yyyy == '1982':
        var = df.loc[:,['dtime', 'sea_ice_area_interpolated']]
        df_all = var.rename({'sea_ice_area_interpolated': yyyy}, axis='columns')
        print(df_all)
    else:
        var = df.loc[:,['sea_ice_area_interpolated']]
        var_re = var.rename({'sea_ice_area_interpolated': yyyy}, axis='columns')
        df_all = pd.concat([df_all, var_re], axis=1)


print(df_all)

##--- Save DataFrame as CSV --------------------------------------##
file_name_csv = out_dir + 'seaice_area_interpolated_oisst_all.csv'
df_all.to_csv(file_name_csv, index=False)

print('')
print('Saved CSV file!')
print('')




##-- Plot Arctic sea ice area ------------------------------------##
print('')
print('Plot Arctic sea ice area')

df_all = df_all.set_index('dtime')

plt.figure(figsize=(20, 12))
#ax = df_all.plot(linewidth = 0.5)
ax = df_all.plot(y='1982', color='blue', linewidth = 0.5, label='1980s')
ax = df_all.plot(y='1983', color='blue', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1984', color='blue', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1985', color='blue', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1986', color='blue', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1987', color='blue', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1988', color='blue', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1989', color='blue', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1990', color='green', linewidth = 0.5, label='1990s', ax=plt.gca())
ax = df_all.plot(y='1991', color='green', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1992', color='green', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1993', color='green', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1994', color='green', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1995', color='green', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1996', color='green', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1997', color='green', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1998', color='green', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='1999', color='green', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2000', color='yellow', linewidth = 0.5, label='2000s', ax=plt.gca())
ax = df_all.plot(y='2001', color='yellow', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2002', color='yellow', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2003', color='yellow', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2004', color='yellow', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2005', color='yellow', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2006', color='yellow', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2007', color='yellow', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2008', color='yellow', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2009', color='yellow', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2010', color='orange', linewidth = 0.5, label='2010s', ax=plt.gca())
ax = df_all.plot(y='2011', color='orange', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2012', color='orange', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2013', color='orange', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2014', color='orange', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2015', color='orange', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2016', color='orange', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2017', color='orange', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2018', color='orange', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2019', color='orange', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2020', color='magenta', linewidth = 0.5, label='2020s', ax=plt.gca())
ax = df_all.plot(y='2021', color='magenta', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2022', color='magenta', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2023', color='magenta', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2024', color='magenta', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2025', color='magenta', linewidth = 0.5, label='', ax=plt.gca())
ax = df_all.plot(y='2026', color='k', linewidth = 1, label='2026', ax=plt.gca())





## Set xlim, ylim ##
span = [datetime.datetime(1982, 1, 1), datetime.datetime(1982, 12, 31)]
#span = [datetime.datetime(1992, 1, 1), datetime.datetime(1992, 12, 31)]
ax.set_xlim(span)
ax.set_ylim([0, 16000000])


### Set X axis ###
ax.set_xlabel('')
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=1)) # interval: 1 month
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b'))


### Set Y axis ###
ax.yaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax.ticklabel_format(style="sci", axis="y", scilimits=(6,6)) # 1*10^6


## Add grid ##
ax.minorticks_on()
ax.grid(color = "gray", linestyle = "--")


## Add legend ##
plt.legend(frameon=False, fontsize = 'small')


## Add labels, title ##
ax.set_ylabel('Sea Ice Area [$\mathrm{km^2}$]', fontsize = 12)

fig_title = 'Daily Arctic sea ice area (OISST)'
plt.suptitle(fig_title, fontsize = 15)



plt.show()
plt.savefig(fig_dir + 'seaice_area_all_oisst.png', dpi=300, format='png')
plt.close('all')

print('')
print('Completed!!')