# Read the daily observed time series of sea ice area/extent in the Arctic
# Convert extracted sea ice data into CSV for each year
#
# Script version 2.0.1
#
# Script works with:
#   Python version 3.10.13
#
#   Package version
#    numpy: 1.26.3 
#    pandas:  2.2.0


# Ver1.0.0: Created by Mariko Koseki, 30.08.2023
# Ver1.0.1: Updated by Mariko, 02.11.2023
# Ver1.1.0: Updated by Mariko, 02.02.2024
# Ver1.1.1: Updated by Mariko, 13.03.2024
# Ver2.0.0: Updated by Mariko, 06.06.2024
# Ver2.0.1: Updated by Mariko, 06.08.2024


import pandas as pd
import numpy as np
import calendar
#import datetime
import sys
import ssl # Ver1.0.1
import os # Ver1.1.0
import platform # Ver1.1.0



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

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


##--- Input -----------------------------------------## #ver2.0.0
args = sys.argv
if len(args) == 2:
    year_start = args[1]
    if len(year_start) == 4:
        yyyy1 = int(year_start[0:4])
        yyyy2 = yyyy1

    else:
        print('---How to use this script---')
        print('python read_obs_ice.py <yyyy>')
        print('or')
        print('python read_obs_ice.py <yyyy1> <yyyy2>')
        print('<yyyy1> = start year; <yyyy2> = end year')
        print('')
        print('Example: python read_obs_ice.py 2023')
        print('Example: python read_obs_ice.py 1995 2023')
        print('')
        sys.exit()



elif len(args) == 3:
    year_start = args[1]
    year_end = args[2]
    if (len(year_start) == 4) & (len(year_end) == 4):
        yyyy1 = int(year_start[0:4])
        yyyy2 = int(year_end[0:4])

    else:
        print('')
        print('---How to use this script---')
        print('python read_obs_ice.py <yyyy>')
        print('or')
        print('python read_obs_ice.py <yyyy1> <yyyy2>')
        print('<yyyy1> = start year; <yyyy2> = end year')
        print('')
        print('Example: python read_obs_ice.py 2023')
        print('Example: python read_obs_ice.py 1995 2023')
        print('')
        sys.exit()


else:
    print('')
    print('---How to use this script---')
    print('python read_obs_ice.py <yyyy>')
    print('or')
    print('python read_obs_ice.py <yyyy1> <yyyy2>')
    print('<yyyy1> = start year; <yyyy2> = end year')
    print('')
    print('Example: python read_obs_ice.py 2023')
    print('Example: python read_obs_ice.py 1995 2023')
    print('')
    sys.exit()





##--- Read sea ice text data from NERSC website -------------------------------##
## Set URL ##
sea_ice_txt = 'https://web.nersc.no/WebData/arctic-roos.org/observation/DailyArcticIceAreaExtent.txt'


## Set path to output file ## # Ver1.1.0
'''
Set 'out_obs' (path to csv file)
'''
seaice_dir = '/nird/projects/NS9873K/norcpm/validation/seaicepre/' # Ver1.1.1
out_obs = seaice_dir + 'obs/'

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



## Convert dataset into Pandas DataFrame ##
ssl._create_default_https_context = ssl._create_unverified_context # Ver1.0.1
df = pd.read_csv(sea_ice_txt, sep="\s+", dtype = {'YEAR':'float', 'DOY':'int', 'AREA':'float', 'EXTENT':'float'}) # DOY: Day of Year




##--- Create list of year ------------------------------------------------------## #ver2.0.0
yyyy_int_list = []
for yyyy in range(yyyy1, yyyy2+1):
    yyyy_int_list.append(yyyy)

print(yyyy_int_list)



##--- Data treatment 1 -------------------------------------------------##
## Change 'YEAR' float to int ##
df['YEAR'] = np.floor(df['YEAR']).astype(int)



## Treat 'YEAR' ##
'''
If YEAR == 2005, 2007, 2008, 2009, 2010, 2011, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023
Data treatment is needed
'''
df.loc[(df['DOY'] == 365) & (df['YEAR'] >= 2005), 'YEAR'] = df['YEAR'] - 1


'''
## Check leap year ##
for y in range(2002,2030):
    if calendar.isleap(y):
        print(y, ' is a leap year')
'''




##--- Extract sea ice data for each year ---------------------------------------------## #ver2.0.0
for YEAR in yyyy_int_list:
    df_copy = df.copy()
    df_copy = df_copy[(df_copy['YEAR'] == YEAR)]



    ##--- Data treatment 2-------------------------------------------------##
    ## Drop duplicates except for the last occurrence ## #ver2.0.0
    df_copy = df_copy.drop_duplicates(subset='DOY', keep='last')


    ## Sort DataFrame by DOY ## #ver2.0.0
    df_copy = df_copy.sort_values('DOY')


    ## Covert DOY into Month and Date and add it into DataFrame ## #ver2.0.0
    df_copy['MMDD'] = pd.to_datetime(df_copy['DOY'], format='%j').dt.strftime('%m-%d')


    '''
    ## Remove 29/02 (DOY = 60) ##
    # ??? No need to remove leap day????
    ## Check leap year ##
    for y in range(2002,2030):
        if calendar.isleap(y):
            print(y, ' is a leap year')


    drop_ind_2004 = ice_data_2004.index[ice_data_2004['DOY'] == 60]; ice_data_2004 = ice_data_2004.drop(drop_ind_2004)
    #print(ice_data_2020) # for debugging
    drop_ind_2008 = ice_data_2008.index[ice_data_2008['DOY'] == 60]; ice_data_2008 = ice_data_2008.drop(drop_ind_2008)
    drop_ind_2012 = ice_data_2012.index[ice_data_2012['DOY'] == 60]; ice_data_2012 = ice_data_2012.drop(drop_ind_2012)
    drop_ind_2016 = ice_data_2016.index[ice_data_2016['DOY'] == 60]; ice_data_2016 = ice_data_2016.drop(drop_ind_2016)
    drop_ind_2020 = ice_data_2020.index[ice_data_2020['DOY'] == 60]; ice_data_2020 = ice_data_2020.drop(drop_ind_2020)
    drop_ind_2024 = ice_data_2024.index[ice_data_2024['DOY'] == 60]; ice_data_2024 = ice_data_2024.drop(drop_ind_2024)
    print(ice_data_2024.head(60)) # for debugging
    #sys.exit()
    '''


    ## Reset index ## #ver2.0.0
    df_copy = df_copy.reset_index()

    



    ##--- Save DataFrame as CSV --------------------------------------## #ver2.0.0
    df_copy.to_csv(out_obs + 'sea_ice_obs_' + str(YEAR) + '.csv', index=False)
    print('Saved sea ice area!')
    print('')
    print('-----------------------')



print('')
print("completed")