Stinsage temperature
import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
#import pyaerocom as pya
import cartopy as ccrs
# Observational record (CRU)
path_crutem = 'shared-ns1000k/inputs/CRU/CRUTEM.4.6.0.0.anomalies.nc'
path_hadcrut = 'shared-ns1000k/inputs/CRU/HadCRUT.4.6.0.0.median.nc'
crutem_data = xr.open_dataset(path_crutem)
hadcrut_data = xr.open_dataset(path_hadcrut)
# Historical model data
# CESM2
path_cesm2 = 'shared-cmip6-for-ns1000k/historical/CESM2/r1i1p1f1/tas_Amon_CESM2_historical_r1i1p1f1_gn_185001-201412.nc'
tas_cesm2_data = xr.open_dataset(path_cesm2)
# NorEsm2-LM
# Get the area weight
path_area_weight_cesm2 = 'shared-cmip6-for-ns1000k/historical/CESM2/r1i1p1f1/areacella_fx_CESM2_historical_r1i1p1f1_gn.nc'
areacella_cesm2 = xr.open_dataset(path_area_weight_cesm2)
# Function to calculate the weighted average
def masked_average(xa, dim=None, weights=None, mask=None):
"""
This function will average
:param xa: dataArray
:param dim: dimension or list of dimensions. e.g. 'lat' or ['lat','lon','time']
:param weights: weights (as xarray)
:param mask: mask (as xarray), True where values to be masked.
:return: masked average xarray
"""
_ = xa.copy()
if mask is not None:
dum, mask_alld = xr.broadcast(xa, mask) # broadcast to all dims
xa_ = xa_.where(np.logical_not(mask))
if weights is not None:
dum, weights_alld = xr.broadcast(xa, weights) # broadcast to all dims
weights_alld = weights_alld.where(np.logical_not(mask_alld))
return (xa_*weights_alld).sum(dim=dim)/weights_alld.sum(dim=dim)
else:
return xa_.mean(dim)
elif weights is not None:
dum, weights_alld = xr.broadcast(xa, weights) # broadcast to all dims
return (xa_*weights_alld).sum(dim)/weights_alld.where(xa_.notnull()).sum(dim=dim)
else:
return xa.mean(dim)
# Yearly average of the data
#daily_t2m_1999 = t2m_1999.resample(time = 'D').mean()
tas_cesm_yearly = tas_cesm2_data['tas'].resample(time = '1Y').mean()
# Global average
tas_cesm_global
# Time series of CESM data