# Compute monthly mean of the IO-SST-V2.1 from NOAA # Calculate horizontal interpolation: 0.25 degree -> 1 degree # # # Script version 1.2 # # Script works with: # Python version 3.11.7 # # Package version # numpy: 1.26.3 # pandas: 2.1.4 # Based on CDO version: 2.4.0 # Based on NCO version: 5.2.3 # #--------------------------------------------- # This script needs "python-cdo" and "pynco" packages # The "python-cdo" and "pynco" packages can be installed: # # $ conda install -c conda-forge python-cdo pynco # #--------------------------------------------- # # # Ver1.0.0: Created by Mariko Koseki, 16.04.2024 # Ver1.1.0: Updated by Mariko, 05.06.2024 # ver1.2: updated by Mariko, 08.12.2025 ##--- Import modules ---------------------------## import numpy as np import pandas as pd import math import datetime from datetime import timedelta #from dateutil.relativedelta import relativedelta import os import glob import sys import platform import calendar import urllib.request import urllib.error import shutil from cdo import * cdo = Cdo() from nco import Nco nco = Nco() from nco.custom import Atted, Rename print('--- Python version ---') print(platform.python_version()) print('--- Package version ---') print('numpy: ', np.__version__) print('pandas: ', pd.__version__) #print('python-cdo version: ', cdo.__version__()) print('Based on CDO version: ', cdo.version()) print('Based on NCO version: ', nco.version()) print('-------------------------') print('') ##--- Input -----------------------------------------## args = sys.argv if len(args) == 2: date_oisst = args[1] if len(date_oisst) == 6: yyyy1 = int(date_oisst[0:4]) mm1 = int(date_oisst[4:6]) yyyy2 = yyyy1 mm2 = mm1 else: print('---How to use this script---') print('python interpolate_oisst.py ') print('or') print('python interpolate_oisst.py ') print(' = start month; = end month') print('') print('Example: python interpolate_oisst.py 202312') print('Example: python interpolate_oisst.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): yyyy1 = int(date_oisst_start[0:4]) mm1 = int(date_oisst_start[4:6]) yyyy2 = int(date_oisst_end[0:4]) mm2 = int(date_oisst_end[4:6]) else: print('---How to use this script---') print('python interpolate_oisst.py ') print('or') print('python interpolate_oisst.py ') print(' = start month; = end month') print('') print('Example: python interpolate_oisst.py 202312') print('Example: python interpolate_oisst.py 199501 202312') print('') sys.exit() else: print('---How to use this script---') print('python interpolate_oisst.py ') print('or') print('python interpolate_oisst.py ') print(' = start month; = end month') print('') print('Example: python interpolate_oisst.py 202312') print('Example: python interpolate_oisst.py 199501 202312') print('') sys.exit() if mm2 == 12: start_month = mm1 start_year = yyyy1 end_month = 1 end_year = yyyy2 + 1 else: start_month = mm1 start_year = yyyy1 end_month = mm2 + 1 end_year = yyyy2 #print(start_month, start_year) #print(end_month, end_year) ##--- Set version ---------------------------------------## ver1.1.0 ver = 'v2.1' #ver = 'v2' ##--- Set path ------------------------------------------## # ver1.2 url = 'https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/' out_month = '/nird/datapeak/NS9873K/norcpm/validation/observation/NOAA/OISST/monthly/' out_org = '/nird/datapeak/NS9873K/norcpm/validation/observation/NOAA/OISST/original/' + ver + '/' #ver1.1.0 out_tmp = '/nird/datapeak/NS9873K/norcpm/validation/observation/NOAA/OISST/tmp/' grid_file = '/nird/datapeak/NS9873K/norcpm/validation/observation/NOAA/OISST/grid/grid_360x180.txt' err_file = '/nird/datapeak/NS9873K/norcpm/validation/observation/NOAA/OISST/tmp/err.nc' if not os.path.exists(out_tmp): os.makedirs(out_tmp) if not os.path.exists(out_month): os.makedirs(out_month) ##--- Set year, month ------------------------------------## start_date = '{:0>4d}-{:0>2d}'.format(start_year,start_month) end_date = '{:0>4d}-{:0>2d}'.format(end_year,end_month) ym_str_list = pd.date_range(start_date, end_date, freq='M').strftime("%Y%m").tolist() #print(ym_str_list) print('start month: ', ym_str_list[0]) print('end month: ', ym_str_list[-1]) ##--- Compute monthly mean from daily OISST ----------------------------------------## print('') print('Start computing monthly mean.....') print('') for ym in ym_str_list: YEAR = ym[0:4] MONTH = ym[4:6] print('') print('month:', ym) print('') out_monmean = out_tmp + 'monthly_mean_' + YEAR + '_' + MONTH + '.nc' out_monmean_reg = out_tmp + 'monthly_mean_reg_' + YEAR + '_' + MONTH + '.nc' file_name_list = sorted(glob.glob(out_org + 'oisst-avhrr-v02r01.' + ym + '*')) #print(file_name_list) ndays = len(file_name_list) #print(ndays) SE = round(1/(math.sqrt(ndays)), 3) print('daily files: ', file_name_list) #print(out_monmean_file) print('') ### Calculate monthly mean ### nco.ncra(input = file_name_list, output = out_monmean) ### Remapping: 0.25 degree -> 1 degree ### cdo.remapcon(grid_file, input = out_monmean, output = out_monmean_reg) ### Calculate standard error and save it as a separate output file (err.nc) ### cdo.mulc(SE, input = '-selvar,err '+out_monmean_reg, output = err_file) # -selvar: select variable ### Append standard error ### opt = ['-v err'] # -v: variable nco.ncks(input = err_file, output = out_monmean_reg, append=True, options=opt) ### Calculate weighted average ### opt2 = ['-O -a zlev'] # -O: Overwrite nco.ncwa(input = out_monmean_reg, output = out_monmean_reg, options=opt2) ### Rename variable: ice -> icec ### rDict = {'ice': 'icec'} nco.ncrename(input = out_monmean_reg, options = [ Rename("variable", rDict) ]) ### Delete FillValue and missing value for icec ### opt3 = ['-a _FillValue,icec,d,, -a missing_value,icec,d,,'] # d: Delete nco.ncatted(input = out_monmean_reg, options = opt3) ### Handle sea ice concentration ### ''' If SST > -2 and icec = -999 -> icec = 0 Sea ice concentration -> % ''' opt4 = ['-O -s \'where(sst>-2. && icec==-999.) icec=0.; where(icec>0.) icec=icec*100.;\''] nco.ncap2(input = out_monmean_reg, output = out_monmean_reg, options = opt4) ### Add FillValue and missing value for icec ### opt5 = ['-a _FillValue,icec,o,f,-999. -a missing_value,icec,o,f,-999.'] # f: Float nco.ncatted(input = out_monmean_reg, options = opt5) ### Add add_offset and scale_factor for sst and icec ### opt6 = ['-a add_offset,sst,o,f,0. -a scale_factor,sst,o,f,1. -a add_offset,icec,o,f,0. -a scale_factor,icec,o,f,1.'] nco.ncatted(input = out_monmean_reg, options = opt6) ### Move NetCDF file to monthly directory ### new_path = shutil.copy(out_monmean_reg, out_month) ### Rename NetCDF file ### old_file_name = out_month + 'monthly_mean_reg_' + YEAR + '_' + MONTH + '.nc' new_file_name = out_month + YEAR + '_' + MONTH + '.nc' os.rename(old_file_name, new_file_name) print('') print('monthly file: ', new_file_name) print('') print('----------') print('') ### Remove temp files ### os.remove(out_monmean) os.remove(out_monmean_reg) print('Completed!')