# Read the daily sea ice concentration from OSTIA
# 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.1
#
# Script works with:
#   Python version 3.10
#
#   Package version
#    numpy:  
#    pandas:  
#    netCDF4:  
#
# Ver1.0: Created by Mariko Koseki, 18.06.2025
# ver1.1, updated by Mariko, 09.01.2026



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_input = args[1]
    if len(date_input) == 6:
        start_year = int(date_input[0:4])
        start_month = int(date_input[4:6])
        end_year = start_year
        end_month = start_month

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

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


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


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





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





##--- Set dataset ID ---------------------------##
'''
Dataset name: Global Ocean OSTIA Sea Surface Temperature and Sea Ice Reprocessed
Product ID: SST_GLO_SST_L4_REP_OBSERVATIONS_010_011
Dataset ID: METOFFICE-GLO-SST-L4-REP-OBS-SST
Spatial resolution: 0.05° × 0.05°
Period: 1 Oct 1981 to 31 May 2022

Dataset name: Global Ocean OSTIA Sea Surface Temperature and Sea Ice Analysis
Product ID: SST_GLO_SST_L4_NRT_OBSERVATIONS_010_001
Dataset ID: METOFFICE-GLO-SST-L4-NRT-OBS-SST-V2
Spatial resolution: 0.05° × 0.05°
Period: 1 Jan 2007 to present
'''

dataset = 'nrt'
#dataset = 'rep'


if dataset == 'rep':
    datasetID = 'METOFFICE-GLO-SST-L4-REP-OBS-SST'
    data_path = 'SST_GLO_SST_L4_REP_OBSERVATIONS_010_011/METOFFICE-GLO-SST-L4-REP-OBS-SST_202003/'
    filename='120000-UKMO-L4_GHRSST-SSTfnd-OSTIA-GLOB_REP-v02.0-fv02.0.nc'

elif dataset == 'nrt':
    datasetID = 'METOFFICE-GLO-SST-L4-NRT-OBS-SST-V2'
    data_path = 'SST_GLO_SST_L4_NRT_OBSERVATIONS_010_001/METOFFICE-GLO-SST-L4-NRT-OBS-SST-V2/'
    filename='120000-UKMO-L4_GHRSST-SSTfnd-OSTIA-GLOB-v02.0-fv02.0.nc'

print('ID: ',datasetID)
print('')



##--- Set path ------------------------------------------##
out_org = '/nird/datapeak/NS9039K/data/external/observation/OSTIA/original/' #ver1.1
out_ice = '/nird/datapeak/NS9873K/norcpm/validation/seaicepre/ostia/monthly/'#ver1.1

grid_area = '/nird/datapeak/NS9039K/data/external/observation/OSTIA/gridarea/gridarea_005deg.nc' #ver1.1

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



##--- 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 ----------##
    print('')
    print('Start computing sea ice area...')
    print('')

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


        nc_ostia = out_org + data_path + YEAR + '/' + MONTH +'/' +YEAR + MONTH + DATE + filename
        

        print('fine name: ', nc_ostia)
        print('')

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

        else:
            print('Data exists!')
            print('')


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



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


        
            ## Read variables ##
            ### time, lon, lat ###
            time = int(nc.variables['time'][:]) #seconds since 1981-01-01 00:00:00
            timeunits = nc.variables['time'].getncattr('units')
            dtime = datetime.date(1981,1,1) + datetime.timedelta(seconds=time)
            print(dtime)
            print('')
            
        
            lat = gc.variables['lat'][:] # Latitude
            lon = gc.variables['lon'][:] # Longitude
            #print(lat); print(lon)
            

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

            ### Sea ice concentration ###
            '''
            byte sea_ice_fraction(time, lat, lon) ;
		    sea_ice_fraction:long_name = "sea ice area fraction" ;
		    sea_ice_fraction:standard_name = "sea_ice_area_fraction" ;
     		sea_ice_fraction:units = "1" ;
	    	sea_ice_fraction:coordinates = "lon lat" ;
    		sea_ice_fraction:_FillValue = -128b ;
	    	sea_ice_fraction:add_offset = 0.f ;
		    sea_ice_fraction:scale_factor = 0.01f ;
     		sea_ice_fraction:valid_min = 0b ;
	    	sea_ice_fraction:valid_max = 100b ;
		    sea_ice_fraction:source = "EUMETSAT OSI-SAF" ;
		    sea_ice_fraction:comment = " Sea ice area fraction" ;
            '''
            ice = np.squeeze(nc.variables['sea_ice_fraction'][:][:][:][:]) # sea ice area fraction, units: %, ice(time, lat, lon)
            #add_offset = nc.variables['sea_ice_fraction'].add_offset; scale_factor = nc.variables['sea_ice_fraction'].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 ##
            '''
            sea_ice_fraction:_FillValue = -128b
            '''
            ice_ar[np.where(ice_ar <= -128)] = 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 ##
            #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 -------------##
    np_arr = np.column_stack([all_time, all_sea_ice_area, all_sea_ice_extent])

    
    ### 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 results as CSV --------------------------------------##
    ### Save DataFrame as CSV ###
    df.to_csv(out_ice + 'seaice_ostia_' + YEAR + MONTH + '.csv', index=False)
    print('Saved sea ice area!')
    print('')


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