<jupyter_start><jupyter_text>Temperature DiagnosticsOnce you copy this repository, feel free to delete this notebook!  Imports<jupyter_code># Display output of plots directly in Notebook
%matplotlib inline
import warnings
warnings.filterwarnings("ignore")

import intake
import numpy as np
import pandas as pd
import xarray as xr
from ncar_jobqueue import NCARCluster
from distributed import Client<jupyter_output><empty_output><jupyter_text>Spin up a Cluster<jupyter_code>cluster = NCARCluster(memory='10 GB')
cluster.scale(20)
client = Client(cluster)
client<jupyter_output><empty_output><jupyter_text>Data Ingest<jupyter_code># Read in the data using xarray or some other package
data_catalog = intake.open_esm_datastore('data/silver-linings-test.json')<jupyter_output><empty_output><jupyter_text>Subset for 2m Temperature<jupyter_code>data_subset = data_catalog.search(frequency='month_1', variable='TREFHT')
data_subset.df<jupyter_output><empty_output><jupyter_text>Read in the dictionary of datasets<jupyter_code>dsets = data_subset.to_dataset_dict(cdf_kwargs={'chunks':{'time':-1}})
ds = dsets['atm.cam.h0.b.e21.BW.f09_g17.SSP245-TSMLT-GAUSS-DEFAULT.003']
ds<jupyter_output><empty_output><jupyter_text>Data OperationFirst, we set up a few helper functions...<jupyter_code>def area_grid(lat, lon):
    """
    Calculate the area of each grid cell
    Area is in square meters
    
    Input
    -----------
    lat: vector of latitude in degrees
    lon: vector of longitude in degrees
    
    Output
    -----------
    area: grid-cell area in square-meters with dimensions, [lat,lon]
    
    Notes
    -----------
    Based on the function in
    https://github.com/chadagreene/CDT/blob/master/cdt/cdtarea.m
    """
    from numpy import meshgrid, deg2rad, gradient, cos
    from xarray import DataArray

    xlon, ylat = meshgrid(lon, lat)
    R = earth_radius(ylat)

    dlat = deg2rad(gradient(ylat, axis=0))
    dlon = deg2rad(gradient(xlon, axis=1))

    dy = dlat * R
    dx = dlon * R * cos(deg2rad(ylat))

    area = dy * dx

    xda = DataArray(
        area,
        dims=["lat", "lon"],
        coords={"lat": lat, "lon": lon},
        attrs={
            "long_name": "area_per_pixel",
            "description": "area per pixel",
            "units": "m^2",
        },
    )
    return xda
def earth_radius(lat):
    '''
    calculate radius of Earth assuming oblate spheroid
    defined by WGS84
    
    Input
    ---------
    lat: vector or latitudes in degrees  
    
    Output
    ----------
    r: vector of radius in meters
    
    Notes
    -----------
    WGS84: https://earth-info.nga.mil/GandG/publications/tr8350.2/tr8350.2-a/Chapter%203.pdf
    '''
    from numpy import deg2rad, sin, cos

    # define oblate spheroid from WGS84
    a = 6378137
    b = 6356752.3142
    e2 = 1 - (b**2/a**2)
    
    # convert from geodecic to geocentric
    # see equation 3-110 in WGS84
    lat = deg2rad(lat)
    lat_gc = np.arctan( (1-e2)*np.tan(lat) )

    # radius equation
    # see equation 3-107 in WGS84
    r = (
        (a * (1 - e2)**0.5) 
         / (1 - (e2 * np.cos(lat_gc)**2))**0.5 
        )

    return r
def center_time(ds):
    """make time the center of the time bounds"""
    ds = ds.copy()
    attrs = ds.time.attrs
    encoding = ds.time.encoding
    tb_name, tb_dim = _get_tb_name_and_tb_dim(ds)
    
    ds['time'] = ds[tb_name].compute().mean(tb_dim).squeeze()
    attrs['note'] = f'time recomputed as {tb_name}.mean({tb_dim})'
    ds.time.attrs = attrs
    ds.time.encoding = encoding
    return ds

def _get_tb_name_and_tb_dim(ds):
    """return the name of the time 'bounds' variable and its second dimension"""
    assert 'bounds' in ds.time.attrs, 'missing "bounds" attr on time'
    tb_name = ds.time.attrs['bounds']        
    assert tb_name in ds, f'missing "{tb_name}"'    
    tb_dim = ds[tb_name].dims[-1]
    return tb_name, tb_dim
def calc_area_weighted_mean(ds):
    ds = center_time(ds.sortby('time'))
    # Do some sort of calculation on the data
    ds_out = (
        (ds.resample(time="AS").mean("time") * da_area).sum(dim=("lat", "lon"))
    ) / total_area
    return ds_out

def convert_to_df(ds):
    return ds.TREFHT.to_series().unstack().T<jupyter_output><empty_output><jupyter_text>Compute the area for the weights<jupyter_code># area dataArray
da_area = area_grid(ds['lat'], ds['lon'])

# total area
total_area = da_area.sum(['lat','lon'])<jupyter_output><empty_output><jupyter_text>Setup which variables to average<jupyter_code>variables = ['TREFHT']<jupyter_output><empty_output><jupyter_text>Run the computation on each dataset<jupyter_code>xr.set_options(keep_attrs=True)
ds_list = []
for key in dsets.keys():
    ds = dsets[key]
    mean = calc_area_weighted_mean(ds)
    out = mean[variables]
    out.attrs['intake_esm_varname'] = variables
    out.attrs['case'] = ds.case
    ds_list.append(out)<jupyter_output><empty_output><jupyter_text>Here we add additional case information<jupyter_code>cases = []
for ds in ds_list:
    cases.append(ds.case)
merged_ds = xr.concat(ds_list, dim='case')
merged_ds['case'] = cases
merged_ds.persist()<jupyter_output><empty_output><jupyter_text>Use datetime index instead of cftime<jupyter_code>datetimeindex = merged_ds.indexes['time'].to_datetimeindex()
merged_ds['time'] = datetimeindex<jupyter_output><empty_output><jupyter_text>Save the dataset<jupyter_code>merged_ds.to_netcdf('data/global_average_temperature.nc')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12940961433010417, 'fraction_numerical': 0.017181502467556207, 'mean_word_length': 3.617721518987342, 'pattern_counts': {'":': 5, '<': 49, '<?xml version=': 0, '>': 49, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '412607', 'license': 'MIT', 'n_tokens_mistral': 1996, 'n_tokens_neox': 1923, 'n_words': 423, 'path': 'notebooks/01_Compute_Averages.ipynb', 'repo_name': 'NCAR/silverlinings-diagnostics', 'size': '5471'}