# Read the daily sea ice concentration from OISST
# Calculate sea ice area in the Arctic region:
#    - Northern Hemisphere
#    - at least 15% sea ice concentration
# Save results as a CSV file
#
# Script version 1.4
#
# Script works with:
#   Python version 3.10.13
#
#   Package version
#    numpy:  1.26.3
#    pandas:  2.2.0
#    netCDF4:  1.6.4
#
# Ver1.0.0: Created by Mariko Koseki, 22.05.2024
# Ver1.1.0: Updated by Mariko, 03.06 2024
# Ver1.2.0: Updated by Mariko, 10.06.2024
# Ver1.3.0: Updated by Mariko, 12.08.2024
# ver1.4: updated by mariko, 08.12.2025


import netCDF4
import numpy as np
import pandas as pd
import datetime
from datetime import timedelta
from dateutil.relativedelta import relativedelta
import calendar
import statistics
import math
import collections
import os
import glob
import sys
import platform


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

print('--- Package version ---')
print('numpy: ', np.__version__)
print('pandas: ', pd.__version__)
print('netCDF4: ', netCDF4.__version__)
print('')
#sys.exit()



##--- Input -----------------------------------------##
args = sys.argv
if len(args) == 2:
    date_oisst = args[1]
    if len(date_oisst) == 6:
        start_year = int(date_oisst[0:4])
        start_month = int(date_oisst[4:6])
        end_year = start_year
        end_month = start_month

    else:
        print('---How to use this script---')
        print('python read_oisst_ice.py <yyyymm>')
        print('or')
        print('python read_oisst_ice.py <yyyymm1> <yyyymm2>')
        print('<yyyymm1> = start month; <yyyymm2> = end month')
        print('')
        print('Example: python read_oisst_ice.py 202312')
        print('Example: python read_oisst_ice.py 199501 202312')
        print('')
        sys.exit()

elif len(args) == 3:
    date_oisst_start = args[1]
    date_oisst_end = args[2]
    if (len(date_oisst_start) == 6) & (len(date_oisst_end) == 6):
        start_year = int(date_oisst_start[0:4])
        start_month = int(date_oisst_start[4:6])
        end_year = int(date_oisst_end[0:4])
        end_month = int(date_oisst_end[4:6])

    else:
        print('---How to use this script---')
        print('python read_oisst_ice.py <yyyymm>')
        print('or')
        print('python read_oisst_ice.py <yyyymm1> <yyyymm2>')
        print('<yyyymm1> = start month; <yyyymm2> = end month')
        print('')
        print('Example: python read_oisst_ice.py 202312')
        print('Example: python read_oisst_ice.py 199501 202312')
        print('')
        sys.exit()


else:
    print('---How to use this script---')
    print('python download_oisst.py <yyyymm>')
    print('or')
    print('python download_oisst.py <yyyymm1> <yyyymm2>')
    print('<yyyymm1> = start month; <yyyymm2> = end month')
    print('')
    print('Example: python download_oisst.py 202312')
    print('Example: python download_oisst.py 199501 202312')
    print('')
    sys.exit()





#print(start_month, start_year)
#print(end_month, end_year)
#print('')



##--- Set version ---------------------------------------## #ver1.1.0
ver = 'v2.1'
#ver = 'v2'



##--- Set path ------------------------------------------## #ver1.1.0, #ver1.4
out_org = '/nird/datapeak/NS9873K/norcpm/validation/observation/NOAA/OISST/original/' + ver + '/'
out_ice_oisst = '/nird/datapeak/NS9873K/norcpm/validation/seaicepre/oisst/' + ver + '/'

grid_area = '/nird/datapeak/NS9873K/norcpm/validation/observation/NOAA/OISST/original/gridarea.nc'

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



##--- Create list of year-month -------------------------##
start_date = '{:0>4d}-{:0>2d}'.format(start_year,start_month)
#print(start_date)
end_date = datetime.datetime(end_year,end_month,1)
end_date = end_date + relativedelta(months=1)
end_date = end_date.strftime('%Y-%m')
#print(end_date)
ym_str_list = pd.date_range(start_date, end_date, freq='ME').strftime('%Y-%m').tolist()

#print(ym_str_list)




##--- Set year, month ------------------------------------##
for ym in ym_str_list:
    year_int = int(ym[0:4])
    month_int = int(ym[5:])
    print('')
    print('Compute sea ice area: ', year_int, '/', month_int)
    print('')


    strdt = datetime.date(year_int, month_int, 1) #1st of this month
    enddt = datetime.date(year_int, month_int, calendar.monthrange(year_int, month_int)[1]) # last date of this month
    days_num = (enddt - strdt).days + 1

    ymd_str_list = []
    for i in range(days_num):
        ymd_str_list.append(strdt + timedelta(days=i))


    print('start date: ', ymd_str_list[0])
    print('end date: ', ymd_str_list[-1])

    #print(ymd_str_list)



    ##--- Calculate sea ice area/extent in the Arctic region ----------## #ver1.2.0
    print('')
    print('Start computing sea ice area...')
    print('')

    all_sea_ice_area = [] #ver1.2.0
    all_sea_ice_extent = [] #ver1.2.0
    all_time = []
    for ymd in ymd_str_list:
        YEAR = ymd.strftime('%Y')
        MONTH = ymd.strftime('%m')
        DATE = ymd.strftime('%d')


        if ver == 'v2.1': #ver1.1.0
            nc_oisst = out_org + 'oisst-avhrr-v02r01.' + YEAR + MONTH + DATE + '.nc'
            if not os.path.exists(nc_oisst): #ver1.3.0
                nc_oisst = out_org + 'oisst-avhrr-v02r01.' + YEAR + MONTH + DATE + '_preliminary.nc'
        else:
            nc_oisst = out_org + 'avhrr-only-v2.' + YEAR + MONTH + DATE + '.nc'

        print('fine name: ', nc_oisst)

        if not os.path.exists(nc_oisst):
            print('no such a directory!!!')
            continue

        else:
            print('Data exists!')


            ## Read netCDF files ##
            nc = netCDF4.Dataset(nc_oisst, 'r') #oisst
            gc = netCDF4.Dataset(grid_area, 'r') #grid cell area



            ## Check length of variables ##
            nlat = len(nc.dimensions['lat']) # Latitude: 720
            nlon = len(nc.dimensions['lon']) # Longitude: 1440 
            #print(nlat, nlon)
        
            nlat_gc = len(gc.dimensions['lat']) # Latitude: 720
            nlon_gc = len(gc.dimensions['lon']) # Longitude: 1440
            #print(nlat_gc, nlon_gc)


        
            ## Read variables ##
            ### time, lon, lat ###
            time = int(nc.variables['time'][:]) # Center time of the day, days since 1978-01-01 12:00:00
            timeunits = nc.variables['time'].getncattr('units')
            dtime = datetime.date(1978,1,1) + datetime.timedelta(days=time)
        
            lat = gc.variables['lat'][:][:] # Latitude
            lon = gc.variables['lon'][:][:] # Longitude
            '''
            print(lat); print(lon); print(dtime)
            '''

            lat_gc = gc.variables['lat'][:][:] # Latitude
            lon_gc = gc.variables['lon'][:][:] # Longitude
            '''
            print(lat_gc); print(lon_gc)
            '''

            ### Sea ice concentration ###
            '''
            short ice(time, zlev, lat, lon) ;
	            ice:long_name = "Sea ice concentration" ;
		    ice:units = "%" ;
                    ice:_FillValue = -999s ;
		    ice:add_offset = 0.f ;
		    ice:scale_factor = 0.01f ;
		    ice:valid_min = 0s ;
		    ice:valid_max = 100s ;
            '''
            ice = np.squeeze(nc.variables['ice'][:][:][:][:]) # Sea ice concentration, units: %, ice(time, zlev, lat, lon)
            #add_offset = nc.variables['ice'].add_offset; scale_factor = nc.variables['ice'].scale_factor
            #ice_real = (ice * scale_factor) + add_offset
            

            ### Grid cell area ###
            '''
            double cell_area(lat, lon) ;
		    cell_area:standard_name = "area" ;
		    cell_area:long_name = "area of grid cell" ;
		    cell_area:units = "m2" ;
            '''
            cell_area = gc.variables['cell_area'][:][:] # Grid Cell Area (m2)
            area = cell_area/1000000 # Grid-Cell Area (km2)


            ### Close netCDF files ###
            nc.close()
            gc.close()


            ## Convert variables into Numpy ndarray ##
            lon_ar = np.array(lon)
            lat_ar = np.array(lat)
            area_ar = np.array(area)
            ice_ar = np.array(ice)


            ## Replace missing values with NaN ##
            '''
            Missing value:
                ice:_FillValue = -999s
            '''
            ice_ar[np.where(ice_ar <= -999.0)] = np.nan

            #print(ice_ar)
            print('Number of NaN: ', np.count_nonzero(np.isnan(ice_ar)))



            ## Calculate sea ice area/extent ##
            '''
            The cumulative area of all polar grid cells of the Northern Hemisphere
                that have at least 15% sea ice concentration
            '''
            ### Check index of the Northern Hemispher ###
            lat_nor_ind = np.where(lat_ar > 0)[0][0]
            lat_sou_ind = np.where(lat_ar <= 0)[0][-1]
            '''
            print(lat_nor_ind)
            print(lat_sou_ind)
            '''

            ### Extract sea ice and grid area only in the Northern Hemispher ###
            ice_nor = ice_ar[lat_nor_ind:, :]
            ice_sou = ice_ar[0:lat_sou_ind+1, :]
            area_nor = area_ar[lat_nor_ind:, :]
            #print('Number of NaN in sea ice (Northern Hemisphere): ', np.count_nonzero(np.isnan(ice_nor)))

            ### Extract grid cells of the Northern Hemisphere that have at least 15% sea ice concentration ###
            #print('Number of grid cells (>= 15%): ', np.count_nonzero(ice_nor >= 0.15))
            #print('Number of grid cells (< 15%): ', np.count_nonzero(ice_nor < 0.15))
            

            ice_nor[np.where(ice_nor < 0.15)] = np.nan
            #print('Number of NaN in sea ice (Northern Hemisphere) + <15%: ', np.count_nonzero(np.isnan(ice_nor)))

           

            ### Calculate sea ice area ###
            sea_ice_area = ice_nor * area_nor
            #print('Number of NaN in sea ice area: ', np.count_nonzero(np.isnan(sea_ice_area)))


            ### Calculate sum of sea ice area ###
            sum_sea_ice_area = np.nansum(sea_ice_area)
            print('Sea ice area [km2]: ', sum_sea_ice_area)


            ### Calculate sea ice extent ### #ver1.2.0
            #print(np.isnan(ice_nor))
            area_nor[np.isnan(ice_nor)] = np.nan
            #print('Number of NaN in sea ice (Northern Hemisphere) + <15%: ', np.count_nonzero(np.isnan(area_nor)))
            sum_sea_ice_extent =  np.nansum(area_nor)
            print('Seaice extent [km2]: ', sum_sea_ice_extent)
            

            ### Add results into list ###
            all_sea_ice_area.append(sum_sea_ice_area)
            all_sea_ice_extent.append(sum_sea_ice_extent)
            all_time.append(dtime)


    #print('')
    #print(all_sea_ice_area)
    #print(all_time)



    ##--- Stack 1-D arrays as columns into a 2-D array -------------## #ver1.2.0
    np_arr = np.column_stack([all_time, all_sea_ice_area, all_sea_ice_extent])



    ##--- Save results as CSV --------------------------------------## #ver1.2.0
    ### Convert all results into Pandas DataFrame ###
    df = pd.DataFrame(np_arr, columns=['dtime', 'sea_ice_area', 'sea_ice_extent'])


    ### Sort dataset by time ###
    df = df.sort_values('dtime');


    ### Save DataFrame as CSV ###
    df.to_csv(out_ice_oisst + 'seaice_oisst_' + YEAR + MONTH + '.csv', index=False)
    print('Saved sea ice area!')
    print('')


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