Click here to launch an interactive online version of this notebook:

Cumulus Friction in the G5NR

regression of SKEdot on (SKE x precip)

Part of this nbviewer repo

---

Sections

1. regression

2. displaying

3. results

---

Imports and setups

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from datetime import datetime
import holoviews as hv
from holoviews import streams
import G5NR_utils

hv.notebook_extension('bokeh')
renderer=hv.renderer('bokeh')

da=xr.open_dataset('SKEdot_merged_90x45_flip.nc').sel(time=slice(datetime(2005,6,1),datetime(2005,9,30))) #for summer

#tricks because stats.linregress takes only 1d arrays
def regression_dataset(varx,vary):
    xname=varx.name
    yname=vary.name

    newda=xr.merge([varx,vary]).stack(lonlat=['lon','lat'])
    gpd=newda.groupby('lonlat')
    regs=list(map(lambda X: stats.linregress(X[1][xname],X[1][yname]),gpd))

    slps_2d=np.array([x[0] for x in regs]).reshape(90,45).T
    intcp_2d=np.array([x[1] for x in regs]).reshape(90,45).T
    corr_2d=np.array([x[2] for x in regs]).reshape(90,45).T
    err_2d=np.array([x[4] for x in regs]).reshape(90,45).T

    slope=xr.DataArray(slps_2d,dims=['lat','lon'],coords={'lon':da.lon,'lat':da.lat})
    slope.name='Slope'

    intercept=xr.DataArray(intcp_2d,dims=['lat','lon'],coords={'lon':da.lon,'lat':da.lat})
    intercept.name='Intercept'

    correlation=xr.DataArray(corr_2d,dims=['lat','lon'],coords={'lon':da.lon,'lat':da.lat})
    correlation.name='Correlation'

    std_err=xr.DataArray(err_2d,dims=['lat','lon'],coords={'lon':da.lon,'lat':da.lat})
    std_err.name='Standard Error'

    return xr.merge([slope,intercept,correlation,std_err])

friction_map=hv.QuadMesh(hv.Image(hv.Dataset(reg_da['Slope']*10),kdims=['lon','lat'],vdims=['Slope']),vdims='Slope',kdims=['lon','lat'],label='Cumulus friction cm^-1')


friction_map=friction_map.redim.range(Slope=(-0.5,0.5))
friction_map=friction_map(plot={'width':800,'height':400,'tools':['hover','tap','box_select'],'colorbar':True,'toolbar':'above'},style={'cmap':'RdBu_r'})
tap_latlon=streams.SingleTap(source=friction_map,x=0,y=0)

def scatter(x,y):

    kedot=da.KEDOT.sel(lon=x,lat=y,method='nearest')
    skeprec=da.SKE2PREC.sel(lon=x,lat=y,method='nearest')*0.1 #to convert Precip part from mm/s to cm/s

    slope, intercep, rval, pval, std = stats.linregress(skeprec, kedot)
    scatter=hv.Points((skeprec.values,kedot.values))

    scatter=scatter.redim.label(x='SKEPREC',y='KEDOT')
    scatter=scatter.redim.unit(SKEPREC='J m-2 cm s-1',KEDOT='W m-2')

    xs = np.linspace(*scatter.range(0)+(2,))
    reg = slope*xs+intercep
    reg_line=hv.Curve((xs, reg),label='')

    scatter_reg=(scatter*reg_line).relabel(label="Lon "+format(x,"0.1f")+" Lat "+format(y,"0.1f"))
    scatter_reg=scatter_reg.relabel(group="slope "+format(slope,"0.2f"))
    return scatter_reg
scatter_plot=hv.DynamicMap(scatter,kdims=[],streams=[tap_latlon])

Results! interactive

hv.opts({'Path':{'style':{'color':'Black'}}}) #to make all coast lines black

friction_map*G5NR_utils.coastline + scatter_plot(plot={'height':400,'width':400,'tools':['hover']},norm={'axiswise':True})