Calculate yearly anomalies for several years on netcdf files in python

I need to calculate monthly, seasonal and annualy anomalies of air temperature of netcdf monthly files of 44 years with some function that allows me to obtain anomalies in the period on a monthly, seasonal and annualy automatically and save the results in a folder. I only know how to do It for one year and not for several years with a function.

from netCDF4 import Dataset, num2date
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeat
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from mpl_toolkits.basemap import Basemap

ds = Dataset('./interim_t2m_19792017.nc')
lats = ds.variables['latitude'][:]  # extract/copy the data
lons = ds.variables['longitude'][:]
time = ds.variables['time']
var = ds.variables['t2m'] 

lon, lat = np.meshgrid(lons, lats)
dates = num2date(time[:], time.units)
dates_pd = pd.to_datetime(dates)
periods = dates_pd.to_period(freq='M')

def plt_map(data):
    m = Basemap(projection='mill',llcrnrlat=-80,urcrnrlat=80,
            llcrnrlon=0,urcrnrlon=360,lat_ts=20,resolution='c')
    x, y = m(lon, lat)
    plt.figure(figsize=(10,7))
    m.drawcoastlines()
    m.drawparallels(np.arange(-80.,81.,20.))
    m.drawmeridians(np.arange(-180.,181.,20.))
    m.drawmapboundary(fill_color='white')
    m.contourf(x,y,data, extend="both",cmap="jet");
    plt.colorbar(orientation='horizontal', pad=0.05)
plt_map(var[0,:,:])

mask_2016 = periods.year==2016
data = var[mask_2016,:,:].mean(axis=0)
plt_map(data)

Answer

I know you are looking for a python answer, but this is bread and butter of CDO (climate data operators), that allows you to do these sort of calculations in one or two commands from the terminal window.

For example, to get the annual means of your Era Interim data you can do

cdo yearmean interim_t2m_19792017.nc erai_yearmean.nc

and then to calculate the annual anomaly you need to do the long term average and subtract it

cdo timmean interim_t2m_19792017.nc erai_timemean.nc
cdo sub erai_yearmean.nc erai_timemean.nc yearanom.nc

You can combine all these above 3 commands using “piping”, but I keep them separate here as it is easier to see what is going on.

you can get the mean monthly seasonal cycle with:

cdo ymonmean interim_t2m_19792017.nc erai_ymonmean.nc

this gives you a file with the average of all the januarys, Feb etc (12 time slices). And then you can calculate the monthly anomaly, each with respect to its own monthly mean with

cdo monmean interim_t2m_19792017.nc erai_monmean.nc
cdo sub erai_monmean.nc erai_ymonmean.nc erai_monanom.nc

There are also functions for seasonal averages.

See the online documentation for further details: https://code.mpimet.mpg.de/projects/cdo/

Lastly, msi_gerva is correct in the comment, it is not clear in the question what the anomalies are with respect to, as you could also calculate monthly anomalies with respect to the annual or long term mean. Moreover, you ask for annual anomalies and say you only know how to do it for one year, but I don’t think that makes much sense, as the anomalies would be zero. It may be helpful to clarify the question more precisely.