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

Understanding eddy momentum flux I: Plan views (maps)

from 7km G5NR data, using G5NRutils.py

Suvarchal Cheedela and Brian Mapes, Oct 2017

Part of this nbviewer repo

---

---

Sections

Open your case file (.zidv) 1. First interactive: u,v,w, uw,vw,uv 2. Second: u’,v’,w’, u’w’, v’w’,u’v’ 3. Third: SKEdot as a function of resolution 4. Fourth: Hover to see value

---

Geoviews is an Earth-specific Holoviews. It displays 1-2 dimensions of a multi-dimensional object, and the other dims become interactive sliders.

xarray is the netCDF-friendly front end for reading in the data. Cartopy has the background maps.

import holoviews as hv
import geoviews as gv
import geoviews.feature as gf
from cartopy import crs as ccrs
import xarray as xr
import numpy as np
from datetime import datetime

Bokeh is a plotting package with pan-zoom interactivity

hv.notebook_extension('bokeh') #other backend is matplotlib but can be changed on the fly just before plottting

Open your 3D data set for a case, as an xarray dataset object

# 7km data grids in local dir from unzipping of .zidv bundle
xrdata_3d = xr.open_dataset('data_0_3D7km30minuteInst.nc')

### Or, if you have no data locally,
# load from a zidv bundle
#import G5NR_utils
#xrdata_3d= G5NR_utils.load_from_zidv('http://weather.rsmas.miami.edu/repository/entry/get/skedot_10.5_lat_40.9_lon_-64.0_time_200508120330.zidv?entryid=a5d8d9fb-0eb1-49bc-9378-a1dd70709dc4')
# remote whole dataset
#xrdata_3d = xr.open_dataset('http://weather.rsmas.miami.edu/repository/opendap/synth:1142722f-a386-4c17-a4f6-0f685cd19ae3:L0c1TlIvR0VPUzUtTmF0dXJlLVJ1bi1JbnN0MzBtaW4tN2ttX3ByZXNzdXJlX21ldGhvZDIubmNtbA==/entry.das')

print(xrdata_3d.data_vars)

Data variables:
    u        (time, lev, lat, lon) float64 -40.06 -39.62 -39.19 -38.81 ...
    v        (time, lev, lat, lon) float64 6.734 6.562 6.39 6.25 6.14 6.109 ...
    w        (time, lev, lat, lon) float64 -0.3402 -0.411 -0.4897 -0.5185 ...
    airdens  (time, lev, lat, lon) float64 2.424e-05 2.424e-05 2.424e-05 ...

Subset data in time and space if desired before lazy read

Use .sel( ) to select by coordinate values, or .isel( ) to select by index. Use slice(value1,value2) to select ranges.

#xrdata_3d.sel(time=slice(datetime(2005,6,1,0,30),datetime(2005,6,1,2,30)),lat=slice(0,5),lon=slice(-14,-10))
#xrdata_3d.isel(time=0)

Merge total fields (U,V,W) with total flux (UW, VW, UV) as a Geoview (Holoview) Image object.

1. First, create and name a new xarray dataset of the product terms.

2. Next, merge these into the xarray set.

3. Finally, convert the merged set into a Geoviews object.

# XARRAY computations
uw=xrdata_3d.u*xrdata_3d.w
uw.name='uw'
vw=xrdata_3d.v*xrdata_3d.w
vw.name='vw'
uv=xrdata_3d.u*xrdata_3d.v
uv.name='uv'

# Merge derived with raw fields
xrdata_3d_all=xr.merge([xrdata_3d,uw,vw,uv])

# GEOVIEWS conversion
subgrid_gv=gv.Dataset(xrdata_3d_all)

# GEOVIEWS display (Image) creation
# dynamic = True is faster to execute but slower to interact: images are made on the fly during interaction
# dynamic = False is slower to execute, but faster to interact: all images are precomputed now

U_img   =subgrid_gv.to(hv.Image,kdims=['lon','lat'],vdims=['u'],label='U',dynamic=True).redim.range(u=(-15,15))
V_img   =subgrid_gv.to(hv.Image,kdims=['lon','lat'],vdims=['v'],label='V',dynamic=True).redim.range(v=(-15,15))
W_img   =subgrid_gv.to(hv.Image,kdims=['lon','lat'],vdims=['w'],label='W',dynamic=True).redim.range(w=(-2,2))
UW_img=subgrid_gv.to(hv.Image,kdims=['lon','lat'],vdims=['uw'],label='UW',dynamic=True).redim.range(uw=(-50,50))
VW_img=subgrid_gv.to(hv.Image,kdims=['lon','lat'],vdims=['vw'],label='VW',dynamic=True).redim.range(vw=(-50,50))
UV_img=subgrid_gv.to(hv.Image,kdims=['lon','lat'],vdims=['uv'],label='UV',dynamic=True).redim.range(vw=(-50,50))

Geoviews multiple-displays creation. So simple!

• (U_img+V_img+W_img+UW_img+VW_img+UV_img).cols(3)

• Options include

  • %%output backend=‘bokeh’ or ‘matplotlib’

  • %%opts Image (cmap=‘RdBu_r’) [width=300 height=200 colorbar=True toolbar=‘above’] Feature.Lines (line_color=‘gray’ line_width=0.5)

  • %%output size=100

%%output backend='bokeh' #backend changed to bokeh, gives more options for tooltip
%%opts Image (cmap='RdBu_r') [width=300 height=200 colorbar=True toolbar='above'] Feature.Lines (line_color='gray' line_width=0.5)
%%output size=90

# Here is the Magic Command to Create the Displays!
(U_img+V_img+W_img+UW_img+VW_img+UV_img).cols(3)     # cols specifies how many columns of plots

lev:

time:

---

G5NR_utils.py has computations with 7km gridded u,v,w inputs

(90,45) is the 4 degree filter scale (divides the world into 90x45 arrays)

• In particular, subgrid returns the filter scale mean to create the prime terms up, vp, wp

from G5NR_utils import subgrid

subgrid_xr = subgrid(xrdata_3d,90,45)

Merge the prime terms, and their products, back into the Subgrid xarray dataset

# XARRAY computations
upwp =subgrid_xr.u*subgrid_xr.w
upwp.name ='upwp'
vpwp =subgrid_xr.v*subgrid_xr.w
vpwp.name ='vpwp'
upvp =subgrid_xr.u*subgrid_xr.v
upvp.name ='upvp'

# Merge these derived fields back with the raw fields
subgrid_xr =xr.merge([subgrid_xr,upwp,vpwp,upvp])

subgrid_xr

<xarray.Dataset>
Dimensions:   (lat: 172, lev: 72, lon: 241, time: 5)
Coordinates:
  * time      (time) datetime64[ns] 2005-05-16T01:30:00 2005-05-16T02:00:00 ...
  * lev       (lev) float64 0.01 0.02 0.0327 0.0476 0.066 0.0893 0.1197 ...
  * lon       (lon) float64 144.5 144.6 144.6 144.7 144.8 144.8 144.9 144.9 ...
  * lat       (lat) float64 23.25 23.31 23.38 23.44 23.5 23.56 23.62 23.69 ...
Data variables:
    u         (time, lev, lat, lon) float64 -1.735 -1.298 -0.8601 -0.4851 ...
    v         (time, lev, lat, lon) float64 3.659 3.487 3.315 3.175 3.065 ...
    w         (time, lev, lat, lon) float64 -0.328 -0.3988 -0.4774 -0.5062 ...
    airdens   (time, lev, lat, lon) float64 2.936e-07 2.918e-07 2.89e-07 ...
    lon_bins  (lat, lon) object (142.0, 146.0] (142.0, 146.0] (142.0, 146.0] ...
    lat_bins  (lat) object (22.5, 26.591] (22.5, 26.591] (22.5, 26.591] ...
    upwp      (time, lev, lat, lon) float64 0.569 0.5174 0.4106 0.2456 ...
    vpwp      (time, lev, lat, lon) float64 -1.2 -1.391 -1.583 -1.607 -1.426 ...
    upvp      (time, lev, lat, lon) float64 -6.349 -4.525 -2.852 -1.54 ...

Make a Geoviews Dataset, and 6 Images, for subgrid (prime) terms

# GEOVIEWS conversion
subgrid_gv = gv.Dataset(subgrid_xr)
subgrid_gv.vdims

[Dimension('u'),
 Dimension('v'),
 Dimension('w'),
 Dimension('airdens'),
 Dimension('lon_bins'),
 Dimension('lat_bins'),
 Dimension('upwp'),
 Dimension('vpwp'),
 Dimension('upvp')]

# GEOVIEWS display (Image) creation
# dynamic = True is faster to execute but slower to interact: images are made on the fly during interaction
# dynamic = False is slower to execute, but faster to interact: all images are precomputed now

u_img    =subgrid_gv.to(hv.Image,kdims=['lon','lat'],vdims=['u'],label='up',dynamic=True).redim.range(u=(-15,15))
v_img    =subgrid_gv.to(hv.Image,kdims=['lon','lat'],vdims=['v'],label='vp',dynamic=True).redim.range(v=(-15,15))
w_img    =subgrid_gv.to(hv.Image,kdims=['lon','lat'],vdims=['w'],label='wp',dynamic=True).redim.range(w=(-2,2))
upwp_img =subgrid_gv.to(hv.Image,kdims=['lon','lat'],vdims=['upwp'],label='upwp',dynamic=True).redim.range(upwp=(-20,20))
vpwp_img =subgrid_gv.to(hv.Image,kdims=['lon','lat'],vdims=['vpwp'],label='vpwp',dynamic=True).redim.range(vpwp=(-20,20))
upvp_img =subgrid_gv.to(hv.Image,kdims=['lon','lat'],vdims=['upvp'],label='upvp',dynamic=True).redim.range(upvp=(-20,20))

Geoviews display for the subgrid scale filtered (“eddy”) products:

%%output backend='bokeh' #backend changed to bokeh, gives more options for tooltip
%%opts Image (cmap='RdBu_r') [width=300 height=200 colorbar=True toolbar='above'] Feature.Lines (line_color='gray' line_width=0.5)
%%output size=90

(u_img + v_img + w_img + upwp_img + vpwp_img + upvp_img).cols(3)

lev:

time:

xrdata_3d.time

<xarray.DataArray 'time' (time: 5)>
array(['2005-05-16T01:30:00.000000000', '2005-05-16T02:00:00.000000000',
       '2005-05-16T02:30:00.000000000', '2005-05-16T03:00:00.000000000',
       '2005-05-16T03:30:00.000000000'], dtype='datetime64[ns]')
Coordinates:
  * time     (time) datetime64[ns] 2005-05-16T01:30:00 2005-05-16T02:00:00 ...
Attributes:
    grads_dim:            t
    grads_mapping:        linear
    grads_size:           36576
    grads_min:            21:30z15may2005
    grads_step:           30mn
    long_name:            time
    minimum:              21:30z15may2005
    maximum:              21z16jun2007
    resolution:           0.0208333
    _CoordinateAxisType:  Time

---

from G5NR_utils import SKEDot, regrid

TIMELEVEL = 2 # hardwired for now

def plot_skedot(nlon,nlat):
    skedot = SKEDot(xrdata_3d.airdens.isel(time=TIMELEVEL),xrdata_3d.u.isel(time=TIMELEVEL),xrdata_3d.v.isel(time=TIMELEVEL),xrdata_3d.w.isel(time=TIMELEVEL),nlon,nlat)
    gv_dataset=gv.Dataset(skedot.SKEDOT).redim.range(SKEDOT=(-1,1)) # W m-2 units

    return hv.Image(gv_dataset)

Resolution as a slider (or menu)

mesh_latlon={'4 deg':(90,45),'2 deg':(180,91),'1 deg':(360,181),'.5 deg':(720,361)}

# A dictionary of images as a function of resolution
img_dict = {res:plot_skedot(mesh_latlon[res][0],mesh_latlon[res][1]) for res in mesh_latlon}

# HoloMap is a static map, very slow to display for heavy dataset
# DynamicMap can be constructed directly from hv.DynamicMap(plot_function) which means calculations are
# done on the fly, so it can be more responsive
# here holomap is used because gives this dropdown widget, havent figured how to do this with adynamic map
# but before displaying holomap it can be converted to a dynamicmap

hmap = hv.HoloMap(img_dict,kdims=['Resolution'])

The plot: SKEdot as a function of resolution

Should match the case selection criterion (e.g. positive SKEdot at center)

%%output backend='bokeh'
%%opts Image (cmap='RdBu_r') [width=600 height=400 colorbar=True toolbar='above']
#hv.util.Dynamic(hmap) #dynamic map
hmap #static map

Resolution:

Bleeding edge in 10/17, needed pip install --upgrade git+https://github.com/ioam/holoviews.git

%%opts QuadMesh (cmap='RdBu_r') [width=600 height=400 colorbar=True xaxis= None yaxis= None tools=['hover'] toolbar='above']
hv.QuadMesh(hmap['4 deg'],kdims=['lon','lat'],vdims=['SKEDOT']).redim.range(SKEDOT=(-1,1))