Search
Stinsage temperature
import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
#import pyaerocom as pya 
import cartopy as ccrs

Loading the data

# 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
 
/opt/conda/lib/python3.7/site-packages/xarray/conventions.py:494: SerializationWarning: variable 'tas' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.
  use_cftime=use_cftime,

Weigthed averages

# 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)
/opt/conda/lib/python3.7/site-packages/xarray/conventions.py:494: SerializationWarning: variable 'areacella' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.
  use_cftime=use_cftime,
# 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)
<xarray.Dataset>
Dimensions:    (lat: 192, lon: 288, nbnd: 2)
Coordinates:
  * lat        (lat) float64 -90.0 -89.06 -88.12 -87.17 ... 88.12 89.06 90.0
  * lon        (lon) float64 0.0 1.25 2.5 3.75 5.0 ... 355.0 356.2 357.5 358.8
Dimensions without coordinates: nbnd
Data variables:
    areacella  (lat, lon) float32 ...
    lat_bnds   (lat, nbnd) float32 ...
    lon_bnds   (lon, nbnd) float32 ...
Attributes:
    Conventions:            CF-1.7 CMIP-6.2
    activity_id:            CMIP
    case_id:                15
    cesm_casename:          b.e21.BHIST.f09_g17.CMIP6-historical.001
    contact:                cesm_cmip6@ucar.edu
    creation_date:          2019-01-16T21:43:39Z
    data_specs_version:     01.00.29
    experiment:             all-forcing simulation of the recent past
    experiment_id:          historical
    forcing_index:          1
    frequency:              fx
    grid:                   native 0.9x1.25 finite volume grid (192x288 latxlon)
    grid_label:             gn
    initialization_index:   1
    institution:            National Center for Atmospheric Research, Climate...
    institution_id:         NCAR
    license:                CMIP6 model data produced by <The National Center...
    mip_era:                CMIP6
    model_doi_url:          https://doi.org/10.5065/D67H1H0V
    nominal_resolution:     100 km
    parent_activity_id:     CMIP
    parent_experiment_id:   piControl
    parent_mip_era:         CMIP6
    parent_source_id:       CESM2
    parent_time_units:      days since 0001-01-01 00:00:00
    parent_variant_label:   r1i1p1f1
    physics_index:          1
    product:                model-output
    realization_index:      1
    realm:                  atmos land
    source:                 CESM2 (2017): atmosphere: CAM6 (0.9x1.25 finite v...
    source_id:              CESM2
    source_type:            AOGCM BGC
    table_id:               fx
    tracking_id:            hdl:21.14100/75033260-7bd2-46b9-9d98-3456f176d588
    variable_id:            areacella
    variant_info:           CMIP6 20th century experiments (1850-2014) with C...
    variant_label:          r1i1p1f1
    sub_experiment:         none
    sub_experiment_id:      none
    branch_time_in_parent:  219000.0
    branch_time_in_child:   674885.0
    branch_method:          standard
    further_info_url:       https://furtherinfo.es-doc.org/CMIP6.NCAR.CESM2.h...
# 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 
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
/opt/conda/lib/python3.7/site-packages/xarray/core/dataarray.py in _getitem_coord(self, key)
    626         try:
--> 627             var = self._coords[key]
    628         except KeyError:

KeyError: 'tas'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
<ipython-input-15-b58f6b6f39cf> in <module>
      1 # Global average
----> 2 tas_cesm_global = tas_cesm_yearly['tas'][:,:,:].mean()

/opt/conda/lib/python3.7/site-packages/xarray/core/dataarray.py in __getitem__(self, key)
    636     def __getitem__(self, key: Any) -> "DataArray":
    637         if isinstance(key, str):
--> 638             return self._getitem_coord(key)
    639         else:
    640             # xarray-style array indexing

/opt/conda/lib/python3.7/site-packages/xarray/core/dataarray.py in _getitem_coord(self, key)
    629             dim_sizes = dict(zip(self.dims, self.shape))
    630             _, key, var = _get_virtual_variable(
--> 631                 self._coords, key, self._level_coords, dim_sizes
    632             )
    633 

/opt/conda/lib/python3.7/site-packages/xarray/core/dataset.py in _get_virtual_variable(variables, key, level_vars, dim_sizes)
    146         ref_var = dim_var.to_index_variable().get_level_variable(ref_name)
    147     else:
--> 148         ref_var = variables[ref_name]
    149 
    150     if var_name is None:

KeyError: 'tas'
tas_cesm_global
<xarray.DataArray 'tas' ()>
array(277.92297, dtype=float32)
# Time series of CESM data 
<xarray.DataArray 'tas' (year: 165, lat: 192, lon: 288)>
array([[[224.4655 , 224.46547, 224.46547, ..., 224.46547, 224.46547,
         224.46547],
        [224.89612, 224.86778, 224.72174, ..., 224.89752, 224.90123,
         224.89964],
        [225.38698, 225.36232, 225.35796, ..., 225.5744 , 225.53694,
         225.48076],
        ...,
        [256.6074 , 256.6331 , 256.657  , ..., 256.54022, 256.56216,
         256.58295],
        [256.29343, 256.30157, 256.30997, ..., 256.26462, 256.27466,
         256.28442],
        [255.96272, 255.96405, 255.96526, ..., 255.95795, 255.95967,
         255.96129]],

       [[225.02513, 225.02513, 225.02513, ..., 225.02513, 225.02513,
         225.02513],
        [225.47853, 225.44695, 225.2979 , ..., 225.48708, 225.48834,
         225.48439],
        [225.96367, 225.9362 , 225.92958, ..., 226.15318, 226.1164 ,
         226.05931],
        ...,
        [256.69946, 256.71533, 256.72977, ..., 256.6564 , 256.67044,
         256.68372],
        [256.5364 , 256.54178, 256.54733, ..., 256.5177 , 256.5243 ,
         256.53058],
        [256.38513, 256.38577, 256.3864 , ..., 256.38254, 256.38348,
         256.3843 ]],

       [[224.43323, 224.43323, 224.43323, ..., 224.43323, 224.43323,
         224.43323],
        [224.79431, 224.76717, 224.6226 , ..., 224.7906 , 224.79597,
         224.7961 ],
        [225.23415, 225.21263, 225.21199, ..., 225.40553, 225.37396,
         225.32281],
        ...,
        [254.9683 , 254.99547, 255.02144, ..., 254.89632, 254.91945,
         254.94197],
        [254.91895, 254.92598, 254.93333, ..., 254.89545, 254.90369,
         254.91154],
        [254.80573, 254.8065 , 254.80719, ..., 254.80293, 254.80396,
         254.80492]],

       ...,

       [[226.8747 , 226.87465, 226.8747 , ..., 226.87465, 226.87463,
         226.87468],
        [227.09892, 227.06682, 226.91841, ..., 227.1061 , 227.10767,
         227.10448],
        [227.57092, 227.54103, 227.53595, ..., 227.75331, 227.71944,
         227.6643 ],
        ...,
        [258.8723 , 258.8749 , 258.8771 , ..., 258.86423, 258.8664 ,
         258.86896],
        [259.0241 , 259.0245 , 259.02524, ..., 259.02356, 259.02414,
         259.02408],
        [259.27103, 259.27042, 259.26987, ..., 259.27322, 259.27243,
         259.2717 ]],

       [[225.6037 , 225.6037 , 225.6037 , ..., 225.6037 , 225.6037 ,
         225.6037 ],
        [225.68701, 225.65607, 225.508  , ..., 225.69873, 225.69786,
         225.693  ],
        [225.88914, 225.86302, 225.85878, ..., 226.0803 , 226.04248,
         225.98372],
        ...,
        [261.38132, 261.3901 , 261.39807, ..., 261.3549 , 261.36337,
         261.37164],
        [261.39246, 261.395  , 261.39743, ..., 261.38327, 261.38666,
         261.38968],
        [261.30823, 261.30875, 261.3092 , ..., 261.30634, 261.307  ,
         261.30765]],

       [[226.04913, 226.04913, 226.04913, ..., 226.04913, 226.04913,
         226.04913],
        [226.4693 , 226.43915, 226.2913 , ..., 226.47685, 226.47803,
         226.47467],
        [226.95909, 226.93245, 226.92613, ..., 227.14818, 227.11066,
         227.05309],
        ...,
        [259.1596 , 259.1677 , 259.17502, ..., 259.13437, 259.14218,
         259.15036],
        [259.24295, 259.24844, 259.2543 , ..., 259.2261 , 259.23203,
         259.23755],
        [259.37125, 259.37112, 259.371  , ..., 259.37183, 259.3716 ,
         259.37143]]], dtype=float32)
Coordinates:
  * lat      (lat) float64 -90.0 -89.06 -88.12 -87.17 ... 87.17 88.12 89.06 90.0
  * lon      (lon) float64 0.0 1.25 2.5 3.75 5.0 ... 355.0 356.2 357.5 358.8
  * year     (year) int64 1850 1851 1852 1853 1854 ... 2010 2011 2012 2013 2014