<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Utilities for large scale climate data analysis
The xarrayutils.utils module contains several functions that have proven useful in several of my day to day projects with both observational and model data.
Linear regression
One of the operations many scientists do is calculating a linear trend along a specified dimension (e.g. time) on each grid point of a dataset. linear_trend makes this very easy. For demonstration purposes lets load some monthly gridded Argo data from APDRC
Step1:  Lets find out how much the salinity in each grid point changed over the full period (20 years)
Step2:  Now we can plot the slope as a map
Step3:  linear_trend converts the dimension over which to integrate into logical indicies, so the units of the plot above are (salinity/timestep of the original product), so here PSS/month.
Correlation maps
But what about a bit more complex task? Lets find out how mixedlayer salinity and temperature correlate. For this we use xr_linregress (for which linear_trend is just a thin wrapper)
Step4:  This works in any dimension the dataset has
Step5:  This map shows that in lower latitudes spatial patterns of salinity are generally anticorrlated with temperature and vice versa in the high latitudes.
Hatch sign agreement
It can often be useful to indicate if the sign along an averaged (or otherwise aggregated) dimension. For instance to show if forced changes have a consistent sign in reference to a multi-member mean. sign_agreement makes this easy. Consider this synthetic example
Step6:  Taking the mean of these fields, suggests that values increase in the upper-left, upper-right, lower-right quadrant, and the missing values in the third layer distort the mean.
Step7:  Lets produce a mask to see where all elements along the member dimension have the same sign
Step8:  You could use this information to indicate the areas of the average, where the members do not agree by hatching
Step9:  Masking values in the mixed layer
Sometimes it is helpful to analyze data by excluding the values in the mixed layer. This can be easily done with mask_mixedlayer. Let's see how
Step10:  We can remove the values in the mixed layer
Step11:  Or to have the mixed layer values only
Step12:  In this case the cell bounds are not available for the model output, but mask_mixedlayer has the option to take those into account and e.g. remove cells only if the lower boundary is within the mixed layer. All you need to do is correctly specify z_bounds with the variable/coordinate name of the cell bounds.
Removing bottom values
For similar reasons you might sometimes want to mask out the values closest to the ocean bottom. remove_bottom_values will mask out the bottom-most grid cell at every location.
<ASSISTANT_TASK:>
Python Code:
import xarray as xr
import numpy as np
%matplotlib inline
path = 'http://apdrc.soest.hawaii.edu:80/dods/public_data/Argo_Products/monthly_mean/monthly_mixed_layer'
ds = xr.open_dataset(path, use_cftime=True)
ds
ds.load()
Explanation: Utilities for large scale climate data analysis
The xarrayutils.utils module contains several functions that have proven useful in several of my day to day projects with both observational and model data.
Linear regression
One of the operations many scientists do is calculating a linear trend along a specified dimension (e.g. time) on each grid point of a dataset. linear_trend makes this very easy. For demonstration purposes lets load some monthly gridded Argo data from APDRC
End of explanation
from xarrayutils.utils import linear_trend
# create an array 
salinity_regressed = linear_trend(ds.mls, 'time')
salinity_regressed
Explanation: Lets find out how much the salinity in each grid point changed over the full period (20 years)
End of explanation
salinity_regressed.slope.plot(robust=True)
Explanation: Now we can plot the slope as a map
End of explanation
from xarrayutils.utils import linear_trend, xr_linregress
tempxsalt = xr_linregress(ds.mlt, ds.mls, dim='time')
tempxsalt.r_value.plot()
Explanation: linear_trend converts the dimension over which to integrate into logical indicies, so the units of the plot above are (salinity/timestep of the original product), so here PSS/month.
Correlation maps
But what about a bit more complex task? Lets find out how mixedlayer salinity and temperature correlate. For this we use xr_linregress (for which linear_trend is just a thin wrapper):
End of explanation
tempxsalt = xr_linregress(ds.mlt, ds.mls, dim='lon')
tempxsalt.r_value.plot(x='time')
Explanation: This works in any dimension the dataset has:
End of explanation
x = np.linspace(-np.pi, np.pi, 25)
y = np.linspace(-np.pi, np.pi, 25)
xx, yy = np.meshgrid(x,y)
data1 = np.sin(xx)
data2 = np.sin(yy)
data3 = np.ones_like(xx)
np.fill_diagonal(data3,np.nan)
np.fill_diagonal(data3[1:],np.nan)
np.fill_diagonal(data3[:,1:],np.nan)
np.fill_diagonal(data3[:,2:],np.nan)
da = xr.DataArray(np.array([data1, data2, data3]), dims=['member','x', 'y'])
da.plot(col='member')
Explanation: This map shows that in lower latitudes spatial patterns of salinity are generally anticorrlated with temperature and vice versa in the high latitudes.
Hatch sign agreement
It can often be useful to indicate if the sign along an averaged (or otherwise aggregated) dimension. For instance to show if forced changes have a consistent sign in reference to a multi-member mean. sign_agreement makes this easy. Consider this synthetic example:
End of explanation
da.mean('member').plot()
Explanation: Taking the mean of these fields, suggests that values increase in the upper-left, upper-right, lower-right quadrant, and the missing values in the third layer distort the mean.
End of explanation
from xarrayutils.utils import sign_agreement
sign_agreement(da, da.mean('member'), 'member', threshold=1.0).plot()
Explanation: Lets produce a mask to see where all elements along the member dimension have the same sign:
End of explanation
da.mean('member').plot()
sign_agreement(
    da, da.mean('member'), 'member'
).plot.contourf(
    colors='none',
    hatches=['..', None],
    levels=[0,0.5],
    add_colorbar=False
)
Explanation: You could use this information to indicate the areas of the average, where the members do not agree by hatching:
End of explanation
import intake
url = "https://raw.githubusercontent.com/NCAR/intake-esm-datastore/master/catalogs/pangeo-cmip6.json"
col = intake.open_esm_datastore(url)
cat = col.search(
    table_id='Omon',
    grid_label='gn',
    experiment_id='historical',
    member_id='r1i1p1f1',
    variable_id=['thetao','mlotst'],#, 
    source_id=["ACCESS-ESM1-5"]
)
ddict = cat.to_dataset_dict(
    zarr_kwargs={'consolidated':True, 'decode_times':True},
)
ds = ddict['CMIP.CSIRO.ACCESS-ESM1-5.historical.Omon.gn']
ds
Explanation: Masking values in the mixed layer
Sometimes it is helpful to analyze data by excluding the values in the mixed layer. This can be easily done with mask_mixedlayer. Let's see how:
First load a CMIP6 dataset from the cloud
End of explanation
from xarrayutils.utils import mask_mixedlayer
ds_wo_ml = mask_mixedlayer(ds, ds.mlotst)
ds_wo_ml
Explanation: We can remove the values in the mixed layer
End of explanation
from xarrayutils.utils import mask_mixedlayer
ds_ml_only = mask_mixedlayer(ds, ds.mlotst, mask='inside')
ds_ml_only
import matplotlib.pyplot as plt
roi = dict(i=150, time=0, j=slice(50, 200))
plt.figure(figsize=[10, 15])
for di, (data, label) in enumerate(zip([ds, ds_wo_ml, ds_ml_only], ['full data', 'mixed layer removed', 'mixed layer only'])):
    plt.subplot(3,1,di+1)
    data.thetao.isel(lev=slice(0,20),**roi).plot(yincrease=False)
    ds.mlotst.isel(**roi).plot(x='j')
    plt.title(label)
Explanation: Or to have the mixed layer values only
End of explanation
from xarrayutils.utils import remove_bottom_values
ds_bot_removed = remove_bottom_values(ds)
ds_bot_removed
roi = dict(i=150, time=0, j=slice(50, 100))
plt.figure(figsize=[10, 15])
for di, (data, label) in enumerate(zip([ds, ds_bot_removed], ['full data', 'bottom layer removed'])):
    plt.subplot(2,1,di+1)
    data.thetao.isel(lev=slice(38, -2),**roi).plot(yincrease=False, robust=True)
    plt.title(label)
    plt.grid()
Explanation: In this case the cell bounds are not available for the model output, but mask_mixedlayer has the option to take those into account and e.g. remove cells only if the lower boundary is within the mixed layer. All you need to do is correctly specify z_bounds with the variable/coordinate name of the cell bounds.
Removing bottom values
For similar reasons you might sometimes want to mask out the values closest to the ocean bottom. remove_bottom_values will mask out the bottom-most grid cell at every location.
End of explanation
<END_TASK>