# Xhistogram Tutorial

Histograms are the foundation of many forms of data analysis.
The goal of xhistogram is to make it easy to calculate weighted histograms in multiple dimensions over n-dimensional arrays, with control over the axes.
Xhistogram builds on top of xarray, for automatic coordiantes and labels, and dask, for parallel scalability.

## Toy Data

We start by showing an example with toy data. First we use xarray to create some random, normally distributed data.

### 1D Histogram


```python
import xarray as xr
import numpy as np
%matplotlib inline

nt, nx = 100, 30
da = xr.DataArray(np.random.randn(nt, nx), dims=['time', 'x'],
                  name='foo') # all inputs need a name
display(da)
da.plot()
```

By default xhistogram operates on all dimensions of an array, just like numpy. However, it operates on xarray DataArrays, taking labels into account.


```python
from xhistogram.xarray import histogram

bins = np.linspace(-4, 4, 20)
h = histogram(da, bins=[bins])
display(h)
h.plot()
```

**TODO:** 
- Bins needs to be a list; this is annoying, would be good to accept single items
- The `foo_bin` coordinate is the estimated bin center, not the bounds. We need to add the bounds to the coordinates, but we can as long as we are returning DataArray and not Dataset.

Both of the above need GitHub Issues

### Histogram over a single axis


```python
h_x = histogram(da, bins=[bins], dim=['time'])
h_x.plot()
```

**TODO:**
  - Relax / explain requirement that dims is always a list


```python
h_x.mean(dim='x').plot()
```

### Weighted Histogram

Weights can be the same shape as the input:


```python
weights = 0.4 * xr.ones_like(da)
histogram(da, bins=[bins], weights=weights)
```

Or can use Xarray broadcasting:


```python
weights = 0.2 * xr.ones_like(da.x)
histogram(da, bins=[bins], weights=weights)
```

## 2D Histogram

Now let's say we have multiple input arrays. We can calculate their joint distribution:


```python
db = xr.DataArray(np.random.randn(nt, nx), dims=['time', 'x'],
                  name='bar') - 2

histogram(da, db, bins=[bins, bins]).plot()
```

## Real Data Example

### Ocean Volume Census in TS Space

Here we show how to use xhistogram to do a volume census of the ocean in Temperature-Salinity Space

First we open the World Ocean Atlas dataset from the opendap dataset (http://apdrc.soest.hawaii.edu/dods/public_data/WOA/WOA13/1_deg/annual). 

Here we read the annual mean Temparature, Salinity and Oxygen on a 5 degree grid.


```python
# Read WOA using opendap 
Temp_url = 'http://apdrc.soest.hawaii.edu:80/dods/public_data/WOA/WOA13/5_deg/annual/temp'
Salt_url = 'http://apdrc.soest.hawaii.edu:80/dods/public_data/WOA/WOA13/5_deg/annual/salt'
Oxy_url = 'http://apdrc.soest.hawaii.edu:80/dods/public_data/WOA/WOA13/5_deg/annual/doxy'

ds = xr.merge([
    xr.open_dataset(Temp_url).tmn.load(),
    xr.open_dataset(Salt_url).smn.load(),
    xr.open_dataset(Oxy_url).omn.load()])
ds
```

Use histogram to bin data points. Use canonical ocean T/S ranges, and bin size of $0.1^0C$, and $0.025psu$. Similar ranges and bin size as this review paper on Mode Waters: https://doi.org/10.1016/B978-0-12-391851-2.00009-X .


```python
sbins = np.arange(31,38, 0.025)
tbins = np.arange(-2, 32, 0.1)
```


```python
# histogram of number of data points
# histogram of number of data points
hTS = histogram(ds.smn, ds.tmn, bins=[sbins, tbins])
np.log10(hTS.T).plot(levels=31)
```

However, we would like to do a volume census, which requires the data points to be weighted by volume of the grid box. 

\begin{equation}
dV = dz*dx*dy
\end{equation}


```python
# histogram of number of data points weighted by volume resolution
# Note that depth is a non-uniform axis

# Create a dz variable
dz = np.diff(ds.lev)
dz =np.insert(dz, 0, dz[0])
dz = xr.DataArray(dz, coords= {'lev':ds.lev}, dims='lev')

# weight by volume of grid cell (resolution = 5degree, 1degree=110km)
dVol = dz * (5*110e3) * (5*110e3*np.cos(ds.lat*np.pi/180)) 

# Note: The weights are automatically broadcast to the right size
hTSw = histogram(ds.smn, ds.tmn, bins=[sbins, tbins], weights=dVol)
np.log10(hTSw.T).plot(levels=31, vmin=11.5, vmax=16, cmap='brg')
```

The ridges of this above plot are indicative of T/S classes with a lot of volume, and some of them are indicative of Mode Waters (example Eighteen Degree water with T$\sim18^oC$, and S$\sim36.5psu$. 

#### Averaging a variable 

Next we calculate the mean oxygen value in each TS bin. 

\begin{equation}
\overline{A} (m,n) = \frac{\sum_{T(x,y,z)=m, S(x,y,z)=n} (A(x,y,z) dV)}{\sum_{T(x,y,z)=m, S(x,y,z)=n}dV}.
\end{equation}


```python
hTSO2 = (histogram(ds.smn.where(~np.isnan(ds.omn)), 
                   ds.tmn.where(~np.isnan(ds.omn)), 
                   bins=[sbins, tbins], 
                   weights=ds.omn.where(~np.isnan(ds.omn))*dVol)/
                histogram(ds.smn.where(~np.isnan(ds.omn)), 
                          ds.tmn.where(~np.isnan(ds.omn)), 
                          bins=[sbins, tbins], 
                          weights=dVol))

(hTSO2.T).plot(vmin=1, vmax=8)
```

Some interesting patterns in average oxygen emerge. Convectively ventilated cold water have the highest oxygen and mode waters have relatively high oxygen. Oxygen minimum zones are interspersed in the middle of volumetic ridges (high volume waters). 

TODO: NaNs in weights are not handled properly, and a few nanvalues make the weighted sum as nan. Maybe use nansum.

## Dask Integration

Should just work, but need examples.

fd1188c3ef83b8c816097be8dc70352918a00601
ipynb
Jupyter Notebook
doc/tutorial.ipynb
raybellwaves/xhistogram
9992eba11eff9ee704e0f6755d521ab56cab9de3
MIT
doc/tutorial.ipynb
raybellwaves/xhistogram
9992eba11eff9ee704e0f6755d521ab56cab9de3
MIT
doc/tutorial.ipynb
raybellwaves/xhistogram
9992eba11eff9ee704e0f6755d521ab56cab9de3
MIT
Qwen/Qwen-72B
1. YES
2. YES