Climatology Functions

Climatology Functions#

Warning

This is not meant to be a standalone notebook. This notebook is part of the process we have for adding entries to the NCL Index and is not meant to be used as tutorial or example code.

Functions covered#

NCL code#

; Note the use of the calendar attribute to assign calendar information

; calcDayAnomTLL
; Adapted from www.ncl.ucar.edu/Document/Functions/Contributed/calcDayanomTLL.shtml
f = addfile("CMIP6_sea_ice_daily_subset.nc", "r")
time = f->time
time@calendar = "noleap"
TIME = cd_calendar(time,0)
year = toint(TIME(:,0))
year@calendar = "noleap"
month = toint(TIME(:,1))
day = toint(TIME(:,2))
ddd = day_of_year(year, month, day)
yyyyddd = year*1000+ddd
yyyyddd@calendar = "noleap"
aice_d = f->aice_d
aiceClmDay = clmDayTLL(aice_d, yyyyddd)
aiceDayAnom = calcDayAnomTLL(aice_d, yyyyddd, aiceClmDay)
print(aiceDayAnom(0,10,10))

; calcMonAnomTLL
; Adapted from www.ncl.ucar.edu/Document/Functions/Contributed/calcMonAnomTLL.shtml
f = addfile("CMIP6_sea_ice_monthly_subset.nc", "r")
aice = f->aice
aiceClmMon = clmMonTLL(aice)
aiceMonAnom = calcMonAnomTLL(aice,aiceClmMon)
print(aiceMonAnom(0,10,10))

; clmDayTLL
; Adapted from www.ncl.ucar.edu/Document/Functions/Contributed/clmDayTLL.shtml
f = addfile("CMIP6_sea_ice_daily_subset.nc", "r")
time = f->time
time@calendar = "noleap"
TIME = cd_calendar(time,0)
year = toint(TIME(:,0))
year@calendar = "noleap"
month = toint(TIME(:,1))
day = toint(TIME(:,2))
ddd = day_of_year(year, month, day)
yyyyddd = year*1000+ddd
yyyyddd@calendar = "noleap"
aice_d = f->aice_d
aiceClmDay = clmDayTLL(aice_d, yyyyddd)
print(aiceClmDay(0,10,10))

; clmMonTLL
; Adapted from www.ncl.ucar.edu/Document/Functions/Contributed/clmMonTLL.shtml
f = addfile("CMIP6_sea_ice_monthly_subset.nc", "r")
aice = f->aice
aiceClmMon = clmMonTLL(aice)
print(aiceClmMon(0,10,10))

; month_to_season
; Adapted from www.ncl.ucar.edu/Document/Functions/Contributed/calcDayAnomTLL.shtml
f = addfile("CMIP6_sea_ice_monthly_subset.nc", "r")
aice = f->aice
aiceSeason = month_to_season(aice, "ASO")
print(aiceSeason(0,10,10))

; rmMonAnnCycTLL
; Adapted from www.ncl.ucar.edu/Document/Functions/Contributed/month_to_season.shtml
f = addfile("CMIP6_sea_ice_monthly_subset.nc", "r")
aice = f->aice
aiceRmMonAnnCyc = rmMonAnnCycTLL(aice)
print(aiceRmMonAnnCyc(0,10,10))

; stdMonTLL
; Adapted from www.ncl.ucar.edu/Document/Functions/Contributed/stdMonTLL.shtml
f = addfile("CMIP6_sea_ice_monthly_subset.nc", "r")
aice = f->aice
aiceStdMon = stdMonTLL(aice)
print(aiceStdMon(0,10,10))

Python Functionality#

import math
import xarray as xr
import geocat.datafiles as gdf
from geocat.comp import month_to_season

results = {}

ds_m = xr.open_dataset(
    gdf.get("applications_files/inputs/CMIP6_sea_ice_monthly_subset.nc")
)
ds_d = xr.open_dataset(
    gdf.get("applications_files/inputs/CMIP6_sea_ice_daily_subset.nc")
)

# calcDayAnomTLL
aice = ds_d.aice_d[:, 10, 10]
DayTLL = aice.groupby(aice.time.dt.dayofyear)
clmDayTLL = DayTLL.mean(dim="time")
calcDayAnomTLL = DayTLL - clmDayTLL
results["calcDayAnomTLL"] = calcDayAnomTLL.item(0)

# calcMonAnomTLL
aice = ds_m.aice[:, 10, 10]
MonTLL = aice.groupby(aice.time.dt.month)
clmMonTLL = MonTLL.mean(dim="time")
calcMonAnomTLL = MonTLL - clmMonTLL
results["calcMonAnomTLL"] = calcMonAnomTLL.item(0)

# clmDayTLL
aice = ds_d.aice_d[:, 10, 10]
DayTLL = aice.groupby(aice.time.dt.dayofyear)
clmDayTLL = DayTLL.mean(dim="time")
results["clmDayTLL"] = clmDayTLL.item(0)

# clmMonTLL
aice = ds_m.aice[:, 10, 10]
MonTLL = aice.groupby(aice.time.dt.month)
clmMonTLL = MonTLL.mean(dim="time")
results["clmMonTLL"] = clmMonTLL.item(0)

# month_to_season
aice = ds_m.aice[:, 10, 10]
mon_to_season = month_to_season(aice, "ASO")
results["month_to_season"] = mon_to_season.item(0)

# rmMonAnnCycTLL
aice = ds_m.aice[:, 10, 10]
MonTLL = aice.groupby(aice.time.dt.month)
clmMonTLL = MonTLL.mean(dim="time")
rmMonAnnCycTLL = MonTLL - clmMonTLL
results["rmMonAnnCycTLL"] = rmMonAnnCycTLL.item(0)

# stdMonTLL
aice = ds_m.aice[:, 10, 10]
MonTLL = aice.groupby(aice.time.dt.month)
stdMonTLL = MonTLL.std(ddof=1)
results["stdMonTLL"] = stdMonTLL.item(0)

Comparison#

# dictionary of climatology functions results
ncl_results = {
    "calcDayAnomTLL": -0.02774364,
    "calcMonAnomTLL": 0.003319263,
    "clmDayTLL": 0.7408184,
    "clmMonTLL": 0.8063018,
    "month_to_season": 0.9931215,
    "rmMonAnnCycTLL": 0.003319263,
    "stdMonTLL": 0.07392025,
}

for c in ncl_results.keys() & results.keys():
    print(f"{c}: \n\tpython:\t\t{results[c]}\n\tncl:\t\t{ncl_results[c]}\n")
    assert math.isclose(ncl_results[c], results[c], rel_tol=1e-06, abs_tol=1e-06)

calcMonAnomTLL: 
	python:		0.0033192038536071777
	ncl:		0.003319263

stdMonTLL: 
	python:		0.07392024993896484
	ncl:		0.07392025

calcDayAnomTLL: 
	python:		-0.02774369716644287
	ncl:		-0.02774364

clmDayTLL: 
	python:		0.7408185005187988
	ncl:		0.7408184

clmMonTLL: 
	python:		0.8063018918037415
	ncl:		0.8063018

month_to_season: 
	python:		0.9931214451789856
	ncl:		0.9931215

rmMonAnnCycTLL: 
	python:		0.0033192038536071777
	ncl:		0.003319263

Differences#

for c in ncl_results.keys() & results.keys():
    print(f"{c}:")
    print(f"\t{results[c] - ncl_results[c]}")

calcMonAnomTLL:
	-5.9146392822079924e-08
stdMonTLL:
	-6.103516303479495e-11
calcDayAnomTLL:
	-5.7166442871126044e-08
clmDayTLL:
	1.0051879884009907e-07
clmMonTLL:
	9.180374149764248e-08
month_to_season:
	-5.482101439469034e-08
rmMonAnnCycTLL:
	-5.9146392822079924e-08
Contents