Table of Contents¶

Under Fred Dupont's suggestion we extracted with option -ip1" followed by 12000

In [226]:

import importlib
import cmocean.cm as cm
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
import xarray as xr

from salishsea_tools import viz_tools

%matplotlib inline

Read Files¶

In [2]:

rpndata = xr.open_dataset('./2014110206_007.nc')
rpndata2 = xr.open_dataset('./2014110206_011.nc')
rpndata3 = xr.open_dataset('./2014110206_012.nc')
opsdata = xr.open_dataset('ops_y2014m11d02.nc')

In [3]:

#rpndata = xr.open_dataset('/data/dlatorne/MEOPAR/GEMLAM-netcdf/2014110206_007.nc')
#opsdata = xr.open_dataset('/results/forcing/atmospheric/GEM2.5/operational/ops_y2014m11d02.nc')

In [4]:

IST, IEN = 110, 365+1
JST, JEN = 20, 285+1

In [5]:

rpndata2

Out[5]:

<xarray.Dataset>
Dimensions:       (time_counter: 1, x: 675, y: 476, z: 1)
Coordinates:
  * time_counter  (time_counter) datetime64[ns] 2014-11-02T17:00:00
Dimensions without coordinates: x, y, z
Data variables:
    FB            (time_counter, z, y, x) float32 ...
    depth         (z) float32 ...
    nav_lat       (y, x) float32 ...
    nav_lon       (y, x) float32 ...
    NT            (time_counter, z, y, x) float32 ...
    PN            (time_counter, z, y, x) float32 ...
    RT            (time_counter, z, y, x) float32 ...
    TD            (time_counter, z, y, x) float32 ...
    TT            (time_counter, z, y, x) float32 ...
    UU            (time_counter, z, y, x) float32 ...
    VV            (time_counter, z, y, x) float32 ...
    PR            (time_counter, z, y, x) float32 ...
    RN            (time_counter, z, y, x) float32 ...

In [6]:

opsdata

Out[6]:

<xarray.Dataset>
Dimensions:       (time_counter: 24, x: 256, y: 266)
Coordinates:
  * time_counter  (time_counter) datetime64[ns] 2014-11-02 ...
  * x             (x) float64 0.0 2.5e+03 5e+03 7.5e+03 1e+04 1.25e+04 ...
  * y             (y) float64 0.0 2.5e+03 5e+03 7.5e+03 1e+04 1.25e+04 ...
Data variables:
    atmpres       (time_counter, y, x) float32 ...
    nav_lat       (y, x) float64 ...
    nav_lon       (y, x) float64 ...
    precip        (time_counter, y, x) float32 ...
    qair          (time_counter, y, x) float32 ...
    solar         (time_counter, y, x) float32 ...
    tair          (time_counter, y, x) float32 ...
    therm_rad     (time_counter, y, x) float32 ...
    u_wind        (time_counter, y, x) float32 ...
    v_wind        (time_counter, y, x) float32 ...
Attributes:
    Conventions:          CF-1.0
    GRIB2_grid_template:  20
    NCO:                  4.4.2
    History:              Thu Dec 31 12:21:07 2015: ncatted -O -a time_origin...

Look at Domains, Interpolation¶

In [119]:

plt.plot(rpndata['nav_lon'][JST:JEN, IST:IEN], rpndata['nav_lat'][JST:JEN, IST:IEN], 'or');
plt.plot(opsdata['nav_lon'], opsdata['nav_lat'], 'ob');
print (rpndata['nav_lon'][JST:JEN, IST:IEN].shape, opsdata['nav_lon'].shape)

(266, 256) (266, 256)

In [237]:

TT = rpndata['TT'][0, 0, JST:JEN, IST:IEN]
mymask = TT > 282 - 273.15
plt.pcolormesh(mymask)

Out[237]:

<matplotlib.collections.QuadMesh at 0x1ce266f60>

Interpolation Code¶

In [8]:

rpn_lon = np.array(rpndata['nav_lon'][JST:JEN, IST:IEN]).flatten()
rpn_lat = np.array(rpndata['nav_lat'][JST:JEN, IST:IEN]).flatten()
points = np.array( (rpn_lon, rpn_lat)).T
print (points.shape)

(68096, 2)

In [9]:

def interpolate_rpn_to_ops(field, points, opsdata):
    xnew = opsdata['nav_lon'][:]
    ynew = opsdata['nav_lat'][:]
    flat_field = np.array(field).flatten()
    rpn_on_ops_field = interpolate.griddata(points, flat_field, (xnew, ynew))
    return rpn_on_ops_field

Temperature¶

Raw Comparison¶

In [227]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
(rpndata['TT'][0, 0, JST:JEN, IST:IEN]+273.15).plot(ax=axs[0], vmax=288, vmin=260, cmap=cm.oxy)
axs[0].set_xlabel('Raw RPN')
opsdata['tair'][13].plot(ax=axs[1], vmax=288, vmin=260)
axs[1].set_xlabel('OPS')
(rpndata['TT'][0, 0, JST:JEN, IST:IEN]+273.15-opsdata['tair'][13]).plot(ax=axs[2], vmax=2, vmin=-2, cmap='bwr')
axs[2].set_xlabel('Difference')

Out[227]:

<matplotlib.text.Text at 0x1b91830f0>

Interpolate¶

In [11]:

field = rpndata['TT'][0, 0, JST:JEN, IST:IEN]+273.15
rpn_TT_on_ops = interpolate_rpn_to_ops(field, points, opsdata)

Interpolated Comparison¶

In [12]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpn_TT_on_ops, vmax=288, vmin=260)
fig.colorbar(colors, ax=axs[0])
opsdata['tair'][13].plot(ax=axs[1], vmax=288, vmin=260)
colors = axs[2].pcolormesh(rpn_TT_on_ops-opsdata['tair'][13], cmap='bwr', vmax=2, vmin=-2)
fig.colorbar(colors, ax=axs[2]);
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated - OPS');

In [13]:

jz0, jz1 = 100, 160
iz0, iz1 = 100, 160
fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpn_TT_on_ops[iz0:iz1, jz0:jz1], vmax=288, vmin=260)
fig.colorbar(colors, ax=axs[0])
opsdata['tair'][13, iz0:iz1, jz0:jz1].plot(ax=axs[1], vmax=288, vmin=260)
colors = axs[2].pcolormesh(rpn_TT_on_ops[iz0:iz1, jz0:jz1] 
                           - opsdata['tair'][13, iz0:iz1, jz0:jz1], cmap='bwr', vmax=2, vmin=-2)
fig.colorbar(colors, ax=axs[2]);
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated - OPS');

Different Time¶

Raw¶

In [14]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
(rpndata2['TT'][0, 0, JST:JEN, IST:IEN]+273.15).plot(ax=axs[0], vmax=288, vmin=260)
axs[0].set_xlabel('Raw RPN')
opsdata['tair'][17].plot(ax=axs[1], vmax=288, vmin=260)
axs[1].set_xlabel('OPS')
(rpndata2['TT'][0, 0, JST:JEN, IST:IEN]+273.15-opsdata['tair'][17]).plot(ax=axs[2], vmax=2, vmin=-2, cmap='bwr')
axs[2].set_xlabel('Difference Raw RPN - OPS');

Interpolated Comparison¶

In [15]:

field = rpndata2['TT'][0, 0, JST:JEN, IST:IEN]+273.15
rpn_TT_on_ops2 = interpolate_rpn_to_ops(field, points, opsdata)

In [16]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpn_TT_on_ops2, vmax=288, vmin=260)
fig.colorbar(colors, ax=axs[0])
opsdata['tair'][17].plot(ax=axs[1], vmax=288, vmin=260)
colors = axs[2].pcolormesh(rpn_TT_on_ops2-opsdata['tair'][17], cmap='bwr', vmax=2, vmin=-2)
fig.colorbar(colors, ax=axs[2]);
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated - OPS');

In [17]:

jz0, jz1 = 100, 160
iz0, iz1 = 100, 160
fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpn_TT_on_ops2[iz0:iz1, jz0:jz1], vmax=288, vmin=260)
fig.colorbar(colors, ax=axs[0])
opsdata['tair'][17, iz0:iz1, jz0:jz1].plot(ax=axs[1], vmax=288, vmin=260)
colors = axs[2].pcolormesh(rpn_TT_on_ops2[iz0:iz1, jz0:jz1] 
                           - opsdata['tair'][17, iz0:iz1, jz0:jz1], cmap='bwr', vmax=2, vmin=-2)
fig.colorbar(colors, ax=axs[2]);
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated - OPS');

One to One (both times)¶

In [18]:

fig, ax = plt.subplots(1, 1, figsize=(4, 4))
ax.plot(np.array(opsdata['tair'][13]).flatten(),rpn_TT_on_ops.flatten(), 'x')
ax.plot(np.array(opsdata['tair'][17]).flatten(),rpn_TT_on_ops2.flatten(), '+')
ax.set_xlabel('OPS')
ax.set_ylabel('Interpolated RPN');

Pressure¶

In [19]:

field = rpndata['PN'][0, 0, JST:JEN, IST:IEN]*100
rpn_PN_on_ops = interpolate_rpn_to_ops(field, points, opsdata)

In [20]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpn_PN_on_ops)
fig.colorbar(colors, ax=axs[0])
opsdata['atmpres'][13].plot(ax=axs[1])
colors = axs[2].pcolormesh(rpn_PN_on_ops-opsdata['atmpres'][13], cmap='bwr', vmax=10, vmin=-10)
fig.colorbar(colors, ax=axs[2]);
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated - OPS');

In [21]:

jz0, jz1 = 100, 160
iz0, iz1 = 100, 160
fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpn_PN_on_ops[iz0:iz1, jz0:jz1])
fig.colorbar(colors, ax=axs[0])
opsdata['atmpres'][13, iz0:iz1, jz0:jz1].plot(ax=axs[1])
colors = axs[2].pcolormesh(rpn_PN_on_ops[iz0:iz1, jz0:jz1]
                           -opsdata['atmpres'][13, iz0:iz1, jz0:jz1], cmap='bwr', vmax=10, vmin=-10)
fig.colorbar(colors, ax=axs[2]);
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated - OPS');

In [22]:

fig, ax = plt.subplots(1, 1, figsize=(4, 4))
ax.plot(np.array(opsdata['atmpres'][13]).flatten(),rpn_PN_on_ops.flatten(), 'x')
ax.set_xlabel('OPS')
ax.set_ylabel('Interpolated RPN');

Solar Data¶

Raw¶

In [23]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
(rpndata2['FB'][0, 0, JST:JEN, IST:IEN]).plot(ax=axs[0], vmax=300)
axs[0].set_xlabel('Raw RPN')
opsdata['solar'][17].plot(ax=axs[1], vmax=300)
axs[1].set_xlabel('OPS')
(rpndata2['FB'][0, 0, JST:JEN, IST:IEN]-opsdata['solar'][17]).plot(ax=axs[2], cmap='bwr')
axs[2].set_xlabel('Difference Raw RPN - OPS');

Interpolated¶

UTC hour 17¶

In [24]:

field = rpndata2['FB'][0, 0, JST:JEN, IST:IEN]
rpn_FB_on_ops = interpolate_rpn_to_ops(field, points, opsdata)

In [25]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpn_FB_on_ops, vmax=300)
fig.colorbar(colors, ax=axs[0])
opsdata['solar'][17].plot(ax=axs[1], vmax=300)
colors = axs[2].pcolormesh(rpn_FB_on_ops-opsdata['solar'][17], cmap='bwr', vmax=150, vmin=-150)
fig.colorbar(colors, ax=axs[2]);
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated - OPS');

UTC hour 18¶

In [26]:

field = rpndata3['FB'][0, 0, JST:JEN, IST:IEN]
rpn_FB_on_ops3 = interpolate_rpn_to_ops(field, points, opsdata)

In [27]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpn_FB_on_ops3, vmax=400)
fig.colorbar(colors, ax=axs[0])
opsdata['solar'][18].plot(ax=axs[1], vmax=400)
colors = axs[2].pcolormesh(rpn_FB_on_ops3 - opsdata['solar'][18], cmap='bwr', vmin=-200, vmax=200)
fig.colorbar(colors, ax=axs[2]);
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated - OPS');

In [28]:

fig, axs = plt.subplots(1, 4, figsize=(16, 4))
axs[0].hist2d(np.array(opsdata['solar'][17]).flatten(),rpn_FB_on_ops.flatten(), bins=24, range=[[0, 400], [0, 400]] )
axs[1].hist2d(np.array(opsdata['solar'][18]).flatten(),rpn_FB_on_ops3.flatten(), bins=24, range=[[0, 400], [0, 400]])
axs[2].hist2d(np.array(0.5*(opsdata['solar'][18] + opsdata['solar'][17])).flatten(),rpn_FB_on_ops.flatten(), bins=24, range=[[0, 400], [0, 400]])
axs[3].hist2d(np.array(opsdata['solar'][18]).flatten(),(0.5*(rpn_FB_on_ops+rpn_FB_on_ops3)).flatten(), bins=24, range=[[0, 400], [0, 400]])
for ax in axs:
    ax.set_xlabel('OPS')
    ax.set_ylabel('Interpolated RPN')
    ax.plot([0, 400], [0, 400], 'w-');

Between Time¶

In [29]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(0.5*(rpn_FB_on_ops3+rpn_FB_on_ops), vmax=400)
fig.colorbar(colors, ax=axs[0])
opsdata['solar'][18].plot(ax=axs[1], vmax=400)
colors = axs[2].pcolormesh(0.5*(rpn_FB_on_ops3+rpn_FB_on_ops) - opsdata['solar'][18], cmap='bwr', vmin=-100, vmax=100)
fig.colorbar(colors, ax=axs[2]);
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated - OPS');

In [30]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpn_FB_on_ops, vmax=400)
fig.colorbar(colors, ax=axs[0])
opsdata['solar'][18].plot(ax=axs[1], vmax=400)
colors = axs[2].pcolormesh(rpn_FB_on_ops - 0.5*(opsdata['solar'][18]+opsdata['solar'][17]), cmap='bwr', vmin=-100, vmax=100)
fig.colorbar(colors, ax=axs[2]);
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated - OPS');

This result makes sense because the OPS numbers are calculated from accumulated whereas RPN must be instantantaneous at end of hour.

Precipitation¶

Units for CORE are kg/m2 = 1000 kg/m3 * ? m/s =

In [31]:

field = (rpndata3['PR'][0, 0, JST:JEN, IST:IEN] - rpndata2['PR'][0, 0, JST:JEN, IST:IEN])/3.6  # volume to mass and per hour
rpn_RT_on_ops = interpolate_rpn_to_ops(field, points, opsdata)

In [32]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpn_RT_on_ops, cmap='plasma')
fig.colorbar(colors, ax=axs[0])
opsdata['precip'][18].plot(ax=axs[1], cmap='plasma', vmax=0.003, vmin=0)
colors = axs[2].pcolormesh(rpn_RT_on_ops-opsdata['precip'][18], cmap='bwr', vmax=0.0002, vmin=-0.0002)
fig.colorbar(colors, ax=axs[2]);
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated - OPS');

In [33]:

fig, ax = plt.subplots(1, 1, figsize=(4, 4))
ax.plot(np.array(opsdata['precip'][18]).flatten(),rpn_RT_on_ops.flatten(), 'x')
ax.plot([0, 0.0035], [0, 0.0035], '-m')
ax.set_xlabel('OPS')
ax.set_ylabel('Interpolated RPN');

Humidity¶

Td of moist air at pressure p and with mixing ratio r is the temperature at which moist air, saturated with respect to water at the given pressure, has a saturation mixing ratio rw equal to the given mixing ratio r.

Atmospheric pressure, p, in hectopascals (hPa); Temperature, t, in degrees Celsius (°C) or T in kelvin (K);

In [34]:

# from WMO 2012 https://www.oceanbestpractices.net/handle/11329/83?show=full
def calculate_q(P, Td):
    P = P / 100 # convert to hectopascals
    # saturation water vapour at the dew point in the pure phase
    ewTd = 6.112 * np.exp(17.62 * Td/(243.12 + Td))
    # epwTd = xvw * p
    xvw = ewTd / P                     
    rwTd = 0.62198 * xvw / (1 - xvw)
    # But at Td r = rw!
    r = rwTd
    q = r/(1+r)
    return q, ewTd                    

In [35]:

# Buck formula  --- Gives identical results.
def calculate_q_b(P, Td):
    P = P / 100 # convert to mb
    # saturation water vapour at the dew point in the pure phase
    ewTd = 6.1121 * np.exp((17.368 - Td/234.5) * (Td /( 238.88 + Td)))
    # epwTd = xvw * p
    xvw = ewTd / P                     
    rwTd = 0.62198 * xvw / (1 - xvw)
    # But at Td r = rw!
    r = rwTd
    q = r/(1+r)
    return q

In [36]:

field = rpndata2['TD'][0, 0, JST:JEN, IST:IEN]  # volume to mass and per hour
rpn_TD_on_ops = interpolate_rpn_to_ops(field, points, opsdata)
field = rpndata2['PN'][0, 0, JST:JEN, IST:IEN]*100
rpn_PN_on_ops = interpolate_rpn_to_ops(field, points, opsdata)

In [37]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpn_TD_on_ops, cmap='plasma')
fig.colorbar(colors, ax=axs[0])
opsdata['qair'][17].plot(ax=axs[1])
axs[2].plot(np.array(opsdata['qair'][17]).flatten(),rpn_TD_on_ops.flatten(), 'x')
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated vs OPS');

In [38]:

rpn_Q_on_ops, rpn_vp_on_ops = calculate_q(rpn_PN_on_ops, rpn_TD_on_ops)
rpn_Qb_on_ops = calculate_q_b(rpn_PN_on_ops, rpn_TD_on_ops)

In [39]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpn_Q_on_ops)
fig.colorbar(colors, ax=axs[0])
opsdata['qair'][17].plot(ax=axs[1])
colors = axs[2].pcolormesh(rpn_Q_on_ops-opsdata['qair'][17], vmax=0.002, vmin=-0.002, cmap='bwr')
fig.colorbar(colors, ax=axs[2]);
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated vs OPS');

In [40]:

fig, axs = plt.subplots(1, 2, figsize=(9, 4))
axs[0].plot(np.array(opsdata['qair'][17]).flatten(),rpn_Q_on_ops.flatten(), 'x')
axs[0].plot(np.array(opsdata['qair'][17]).flatten(),rpn_Qb_on_ops.flatten(), '+')
axs[0].plot([0, 0.01], [0, 0.01], '-m')
axs[0].set_xlabel('OPS')
axs[0].set_ylabel('Interpolated RPN');

axs[1].hist2d(np.array(opsdata['qair'][17]).flatten(),
                        rpn_Q_on_ops.flatten() - np.array(opsdata['qair'][17]).flatten(), vmax=5e3)
axs[1].set_xlabel('OPS')
axs[1].set_ylabel('Interpolated RPN - OPS');

In [41]:

fig, ax = plt.subplots(1, 1, figsize=(4, 4))
ax.hist2d(rpn_TT_on_ops2.flatten(), -np.array(opsdata['qair'][17]).flatten()+rpn_Q_on_ops.flatten(), vmax=5e3)
ax.set_xlabel('OPS')
ax.set_ylabel('Interpolated RPN');

Looks fine over the water.

Long Wave¶

In [179]:

sigma = 5.6697e-8

In [180]:

field = rpndata2['NT'][0, 0, JST:JEN, IST:IEN]  # volume to mass and per hour
rpn_NT_on_ops2 = interpolate_rpn_to_ops(field, points, opsdata)

In [181]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpn_NT_on_ops, cmap='plasma')
fig.colorbar(colors, ax=axs[0])
opsdata['therm_rad'][17].plot(ax=axs[1])
axs[2].plot(np.array(opsdata['therm_rad'][17]).flatten(),rpn_NT_on_ops.flatten(), 'x')
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated vs OPS');

In [239]:

def calculate_stats(a, b, mymask):
    bias = np.ma.array(b-a, mask=mymask).mean()
    rmse = np.sqrt(np.ma.array((b - a)**2, mask=mymask).mean())
    return bias, rmse

Longwave from SOG¶

In [189]:

# from sog
def sog_longwave(rpn_NT_on_ops, rpn_TT_on_ops, sigma):
    # cf  cloud fraction 
    # atemp_value  current air temperature [K]

    atemp_value = rpn_TT_on_ops
    cf = rpn_NT_on_ops
    
    r = 0.03
    Ce = 9.37e-6

    # Downward radiation from atmosphere
    # *** There are several different ways to calculate this
    lw_in = ((1 - r) * (1 + 0.170 * cf**2) * Ce * atemp_value**2 
         * sigma * atemp_value**4)

    return lw_in

Swinbank Unsworth¶

In [176]:

def swinbank_unsworth(rpn_NT_on_ops, rpn_TT_on_ops, sigma):
    To = rpn_TT_on_ops2
    Lclr = 5.31e-13*To**6   # Swinbank
    eclr = Lclr/(sigma * To**4)
    cf = rpn_NT_on_ops
    ewc = (1 - 0.84*cf)*eclr + 0.84*cf # Unsworth
    Lwc = ewc * sigma * To**4
    return Lwc

Prata Unsworth¶

In [183]:

def prata_unsworth(rpn_vp_on_ops, rpn_NT_on_ops, rpn_TT_on_ops, sigma):
    eo = rpn_vp_on_ops/1000.
    To = rpn_TT_on_ops2
    w = 465*eo/To  
    eclr = 1 - (1 + w) * np.exp(-(1.2 + 3*w)**0.5)
    cf = rpn_NT_on_ops
    ewc = (1 - 0.84*cf)*eclr + 0.84*cf # Unsworth
    Lwc = ewc * sigma * To**4
    return Lwc

Dilley Unsworth¶

In [243]:

def dilley_unsworth(rpn_vp_on_ops, rpn_NT_on_ops, rpn_TT_on_ops, sigma):
# https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2008WR007394
# Dilley clear sky algorithm
    eo = rpn_vp_on_ops/1000.
    To = rpn_TT_on_ops2
    w = 465*eo/To
    Lclr = 59.38 + 113.7*(To/273.16)**6 + 96.96*np.sqrt(w/2.5) # Dilley
    eclr = Lclr/(sigma * To**4)
# Unsworth
    cf = rpn_NT_on_ops
    ewc = (1 - 0.84*cf)*eclr + 0.84*cf
#
    Lwc = ewc * sigma * To**4
    return Lwc

Berliand from NEMO CLIO¶

In [218]:

def Berliand(rpn_lat, rpn_TT_on_ops, rpn_NT_on_ops, rpn_vp_on_ops, sigma):
#DATA budyko / 1.00, 0.98, 0.95, 0.92, 0.89, 0.86, 0.83, 0.80, 0.78, 0.75,   &
#      &          0.72, 0.69, 0.67, 0.64, 0.61, 0.58, 0.56, 0.53, 0.50 /
    budyko = [1.00, 0.98, 0.95, 0.92, 0.89, 0.86, 0.83, 0.80, 0.78, 0.75,   
               0.72, 0.69, 0.67, 0.64, 0.61, 0.58, 0.56, 0.53, 0.50]
    sbudyko = np.zeros_like(rpn_lat)
# and the correction factor for taking into account  the effect of clouds 
    for jj in range(rpn_lat.shape[1]):
        for ji in range(rpn_lat.shape[0]):
#         DO jj = 1, jpj
#            DO ji = 1 , jpi
#    zalat          = ( 90.e0 - ABS( gphit(ji,jj) ) ) /  5.e0
            zalat = (90 - np.abs(np.array(rpn_lat[ji, jj]))) / 5
#               indxb          = 1 + INT( zalat )
            indxb = zalat.astype(int)
# correction factor to account for the effect of clouds
            sbudyko[ji,jj] = budyko[indxb]
#--------------------------------------!
#  long-wave radiation over the ocean  !  ( Berliand 1952 ; all latitudes )
#--------------------------------------!
#ztatm = sf(jp_tair)%fnow(ji,jj,1)
#ztatm3  = ztatm * ztatm * ztatm
    ztatm4 = rpn_TT_on_ops**4
    zcldeff = 1.0 - sbudyko * rpn_NT_on_ops**2
    ztaevbk = ztatm4 * zcldeff * ( 0.39 - 0.05 * np.sqrt(rpn_vp_on_ops) ) 
    zqlw = 0.97 * sigma * ( ztaevbk + 4. * ztatm4 ) / 4.
    return zqlw

Parkison and Washington¶

In [222]:

def parkinsonwashington(cf, TT, sigma):
    Lw = sigma * TT**4 * (1 
                          - 0.261 * np.exp(-7.77e-4 * (273 - TT)**2)) * (1 + 0.275*cf) 
    return Lw

In [185]:

field = rpndata2['NT'][0, 0, JST:JEN, IST:IEN]  # volume to mass and per hour
rpn_NT_on_ops2 = interpolate_rpn_to_ops(field, points, opsdata)
field = rpndata2['PN'][0, 0, JST:JEN, IST:IEN]*100
rpn_PN_on_ops2 = interpolate_rpn_to_ops(field, points, opsdata)
field = rpndata2['TD'][0, 0, JST:JEN, IST:IEN]  # volume to mass and per hour
rpn_TD_on_ops2 = interpolate_rpn_to_ops(field, points, opsdata)
rpn_Q_on_ops2, rpn_vp_on_ops2 = calculate_q(rpn_PN_on_ops2, rpn_TD_on_ops2)

In [244]:

clio = Berliand(rpndata2['nav_lat'][JST:JEN, IST:IEN], rpn_TT_on_ops2, rpn_NT_on_ops2, rpn_vp_on_ops2, sigma)
sog = sog_longwave(rpn_NT_on_ops2, rpn_TT_on_ops2, sigma)
swinbank = swinbank_unsworth(rpn_NT_on_ops2, rpn_TT_on_ops2, sigma)
prata = prata_unsworth(rpn_vp_on_ops2, rpn_NT_on_ops2, rpn_TT_on_ops2, sigma)
dilley = dilley_unsworth(rpn_vp_on_ops2, rpn_NT_on_ops2, rpn_TT_on_ops2, sigma)
pw = parkinsonwashington(rpn_NT_on_ops2, rpn_TT_on_ops2, sigma)

In [245]:

fig, axs = plt.subplots(6, 3, figsize=(18, 24))

colors = axs[0, 0].pcolormesh(clio, vmax=375, vmin=175)
fig.colorbar(colors, ax=axs[0, 0])
colors = axs[0, 1].pcolormesh(clio - np.array(opsdata['therm_rad'][17]), vmax=100, vmin=-100, cmap='bwr');
fig.colorbar(colors, ax=axs[0, 1])
axs[0, 2].hist2d(np.array(opsdata['therm_rad'][17]).flatten(), clio.flatten(), bins=12, range=[[175, 375], [175, 375]])
axs[0, 2].plot([175, 375], [175, 375], 'm-')
axs[0, 2].set_title('CLIO')

colors = axs[1, 0].pcolormesh(sog, vmax=375, vmin=175)
fig.colorbar(colors, ax=axs[1, 0])
colors = axs[1, 1].pcolormesh(sog - np.array(opsdata['therm_rad'][17]), vmax=100, vmin=-100, cmap='bwr');
fig.colorbar(colors, ax=axs[1, 1])
axs[1, 2].hist2d(np.array(opsdata['therm_rad'][17]).flatten(), sog.flatten(), bins=12, range=[[175, 375], [175, 375]])
axs[1, 2].plot([175, 375], [175, 375], 'm-')
axs[1, 2].set_title('SOG')

colors = axs[2, 0].pcolormesh(swinbank, vmax=375, vmin=175)
fig.colorbar(colors, ax=axs[2, 0])
colors = axs[2, 1].pcolormesh(swinbank - np.array(opsdata['therm_rad'][17]), vmax=100, vmin=-100, cmap='bwr');
fig.colorbar(colors, ax=axs[2, 1])
axs[2, 2].hist2d(np.array(opsdata['therm_rad'][17]).flatten(), swinbank.flatten(), bins=12, range=[[175, 375], [175, 375]])
axs[2, 2].plot([175, 375], [175, 375], 'm-')
axs[2, 2].set_title('Swinbank')

colors = axs[3, 0].pcolormesh(prata, vmax=375, vmin=175)
fig.colorbar(colors, ax=axs[3, 0])
colors = axs[3, 1].pcolormesh(prata - np.array(opsdata['therm_rad'][17]), vmax=100, vmin=-100, cmap='bwr');
fig.colorbar(colors, ax=axs[3, 1])
axs[3, 2].hist2d(np.array(opsdata['therm_rad'][17]).flatten(), prata.flatten(), bins=12, range=[[175, 375], [175, 375]])
axs[3, 2].plot([175, 375], [175, 375], 'm-')
axs[3, 2].set_title('Prata')

colors = axs[4, 0].pcolormesh(dilley, vmax=375, vmin=175)
fig.colorbar(colors, ax=axs[4, 0])
colors = axs[4, 1].pcolormesh(dilley - np.array(opsdata['therm_rad'][17]), vmax=100, vmin=-100, cmap='bwr');
fig.colorbar(colors, ax=axs[4, 1])
axs[4, 2].hist2d(np.array(opsdata['therm_rad'][17]).flatten(), dilley.flatten(), bins=12, range=[[175, 375], [175, 375]])
axs[4, 2].plot([175, 375], [175, 375], 'm-')
axs[4, 2].set_title('Dilley');

colors = axs[5, 0].pcolormesh(pw, vmax=375, vmin=175)
fig.colorbar(colors, ax=axs[5, 0])
colors = axs[5, 1].pcolormesh(pw - np.array(opsdata['therm_rad'][17]), vmax=100, vmin=-100, cmap='bwr');
fig.colorbar(colors, ax=axs[5, 1])
axs[5, 2].hist2d(np.array(opsdata['therm_rad'][17]).flatten(), pw.flatten(), bins=12, range=[[175, 375], [175, 375]])
axs[5, 2].plot([175, 375], [175, 375], 'm-')
axs[5, 2].set_title('Parkinson and Washington');

In [246]:

print ('clio', calculate_stats(clio, opsdata['therm_rad'][17], mymask))
print ('sog', calculate_stats(sog, opsdata['therm_rad'][17], mymask))
print ('swinbank',calculate_stats(swinbank, opsdata['therm_rad'][17], mymask))
print ('prata', calculate_stats(prata, opsdata['therm_rad'][17], mymask))
print ('dilley', calculate_stats(dilley, opsdata['therm_rad'][17], mymask))
print ('parkinson washington', calculate_stats(pw, opsdata['therm_rad'][17], mymask))

clio (-52.192977814342356, 70.151511703530915)
sog (30.670266881764679, 38.951972590125393)
swinbank (0.95632445237785202, 15.47978490074234)
prata (8.2796261449058495, 18.075544155914653)
dilley (25.553302493046417, 34.237285177081581)
parkinson washington (-3.6351386709523257, 17.016615342179801)

In [54]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
axs[0].hist2d(w.flatten(), 
                     np.array(opsdata['therm_rad'][17]).flatten()-Lwc.flatten(), vmax=5e3);
axs[1].hist2d(rpn_NT_on_ops.flatten(), 
                     np.array(opsdata['therm_rad'][17]).flatten()-Lwc.flatten(), vmax=5e3);
axs[2].hist2d(rpn_TT_on_ops.flatten(), 
                     np.array(opsdata['therm_rad'][17]).flatten()-Lwc.flatten(), vmax=5e3);

In [55]:

# https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2008WR007394
# Dilley clear sky algorithm
eo = rpn_vp_on_ops/1000.
To = rpn_TT_on_ops2
w = 465*eo/To
Lclr = 59.38 + 113.7*(To/273.16)**6 + 96.96*np.sqrt(w/2.5)
Lclr = 5.31e-13*To**6   # Swinbank
sigma = 5.6697e-8
eclr = Lclr/(sigma * To**4)
# Unsworth
cf = rpn_NT_on_ops
#ewc = (1 - 0.84*cf)*eclr + 0.84*cf
ewc = eclr * (1 + 0.17*cf**2)
Lwc = ewc * sigma * To**4

In [56]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(Lwc, vmax=375, vmin=175)
fig.colorbar(colors, ax=axs[0])
opsdata['therm_rad'][17].plot(ax=axs[1], vmax=375, vmin=175);
axs[2].plot(np.array(opsdata['therm_rad'][17]).flatten(), Lwc.flatten(), 'x')
axs[2].plot([175, 375], [175, 375], 'm-')
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated vs OPS');

In [57]:

import scipy.stats as stats
slope, intercept, r_value, p_value, std_err = stats.linregress(np.array(opsdata['therm_rad'][17]).flatten(), 
                       Lwc.flatten())

In [58]:

print (slope, intercept, r_value, p_value, std_err)

0.840113760758 28.6746921964 0.916138366052 0.0 0.00140869734966

In [59]:

slope, intercept, r_value, p_value, std_err = stats.linregress(cf.flatten(), 
                    Lwc.flatten() - np.array(opsdata['therm_rad'][17]).flatten())

In [60]:

print (slope, intercept, r_value, p_value, std_err)

-20.7416111826 -6.18941675955 -0.406587875756 0.0 0.178605687356

In [61]:

Lwc_c = Lwc - (intercept + slope*cf)

In [62]:

print (calculate_stats(np.array(opsdata['therm_rad'][17]).flatten(), Lwc.flatten()))

(-20.321096193682706, 24.105825224854367)

In [63]:

print (calculate_stats(np.array(opsdata['therm_rad'][17]).flatten(), Lwc_c.flatten()))

(-1.2459163372113835e-05, 15.763028368267756)

In [64]:

fig, axs = plt.subplots(1, 3, figsize=(18, 4))
axs[0].hist2d(w.flatten(), 
                     np.array(opsdata['therm_rad'][17]).flatten()-Lwc_c.flatten(), vmax=5e3);
axs[1].hist2d(rpn_NT_on_ops.flatten(), 
                     np.array(opsdata['therm_rad'][17]).flatten()-Lwc_c.flatten(), vmax=5e3);
axs[2].hist2d(rpn_TT_on_ops.flatten(), 
                     np.array(opsdata['therm_rad'][17]).flatten()-Lwc_c.flatten(), vmax=5e3);

Wind¶

In [65]:

field = rpndata3['UU'][0, 0, JST:JEN, IST:IEN]
# convert to m/s from knots
rpn_UU_on_ops = interpolate_rpn_to_ops(field, points, opsdata)*0.514444
field = rpndata3['VV'][0, 0, JST:JEN, IST:IEN]
rpn_VV_on_ops = interpolate_rpn_to_ops(field, points, opsdata)*0.514444

In [66]:

fig, axs = plt.subplots(3, 3, figsize=(18, 13.5))
colors = axs[0, 0].pcolormesh(rpn_UU_on_ops, cmap='bwr')
fig.colorbar(colors, ax=axs[0, 0])
opsdata['u_wind'][18].plot(ax=axs[0, 1], cmap='bwr')
colors = axs[0, 2].pcolormesh(rpn_UU_on_ops-opsdata['u_wind'][18], cmap='bwr')
fig.colorbar(colors, ax=axs[0, 2]);
axs[0, 0].set_title('RPN Interpolated')
axs[0, 1].set_title('OPS')
axs[0, 2].set_title('RPN Interpolated - OPS');

colors = axs[1, 0].pcolormesh(rpn_VV_on_ops, cmap='bwr')
fig.colorbar(colors, ax=axs[1, 0])
opsdata['v_wind'][18].plot(ax=axs[1, 1], cmap='bwr')
colors = axs[1, 2].pcolormesh(rpn_VV_on_ops-opsdata['v_wind'][18], cmap='bwr')
fig.colorbar(colors, ax=axs[1, 2]);
axs[1, 0].set_title('RPN Interpolated')
axs[1, 1].set_title('OPS')
axs[1, 2].set_title('RPN Interpolated - OPS');

rpnspeed = np.sqrt(rpn_UU_on_ops**2 + rpn_VV_on_ops**2)
opsspeed = np.sqrt(opsdata['u_wind'][18]**2 + opsdata['v_wind'][18]**2)
colors = axs[2, 0].pcolormesh(rpnspeed , cmap='plasma', vmax=12)
fig.colorbar(colors, ax=axs[2, 0])
opsspeed.plot(ax=axs[2, 1], cmap='plasma', vmax=12)
colors = axs[2, 2].pcolormesh(rpnspeed - opsspeed, cmap='bwr', vmax=1, vmin=-1)
fig.colorbar(colors, ax=axs[2, 2]);
axs[1, 0].set_title('RPN Interpolated')
axs[1, 1].set_title('OPS')
axs[1, 2].set_title('RPN Interpolated - OPS');

In [67]:

jz0, jz1 = 100, 160
iz0, iz1 = 100, 160
fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpnspeed[iz0:iz1, jz0:jz1], cmap='plasma', vmax=12)
fig.colorbar(colors, ax=axs[0])
opsspeed[iz0:iz1, jz0:jz1].plot(ax=axs[1], vmax=12, cmap='plasma')
colors = axs[2].pcolormesh(rpnspeed[iz0:iz1, jz0:jz1] 
                           - opsspeed[iz0:iz1, jz0:jz1], cmap='bwr', 
                           vmax=1, vmin=-1)
fig.colorbar(colors, ax=axs[2]);
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated - OPS');

In [68]:

jz0, jz1 = 40, 100
iz0, iz1 = 40, 100
fig, axs = plt.subplots(1, 3, figsize=(18, 4))
colors = axs[0].pcolormesh(rpnspeed[iz0:iz1, jz0:jz1], cmap='plasma', 
                           vmin=0, vmax=12)
fig.colorbar(colors, ax=axs[0])
opsspeed[iz0:iz1, jz0:jz1].plot(ax=axs[1], vmin=0, vmax=12, cmap='plasma')
colors = axs[2].pcolormesh(rpnspeed[iz0:iz1, jz0:jz1] 
                           - opsspeed[iz0:iz1, jz0:jz1], cmap='bwr', 
                           vmax=1, vmin=-1)
fig.colorbar(colors, ax=axs[2]);
axs[0].set_title('RPN Interpolated')
axs[1].set_title('OPS')
axs[2].set_title('RPN Interpolated - OPS');

In [69]:

plt.hist2d(np.array(opsspeed).flatten(), rpnspeed.flatten(), vmax=1e3);
plt.plot([0, 12], [0, 12], 'm')
plt.colorbar()

Out[69]:

<matplotlib.colorbar.Colorbar at 0x11bc0ea58>

In [115]:

expand_lat = np.empty((1, rpndata3['nav_lat'].shape[0], rpndata3['nav_lat'].shape[1]))
expand_lat[0] = rpndata3['nav_lat']
expand_lon = np.empty((1, rpndata3['nav_lon'].shape[0], rpndata3['nav_lon'].shape[1]))
expand_lon[0] = rpndata3['nav_lon']
coords = {'glamu': expand_lon,
         'gphiu': expand_lat}

In [138]:

importlib.reload(viz_tools)
rpn_UU = rpndata3['UU'][0, 0, 1:, 1:]
rpn_VV = rpndata3['VV'][0, 0, 1:, 1:]
print (rpn_UU.shape, rpn_VV.shape)
print (expand_lon.shape, expand_lat.shape)
u_out, v_out, a, b, c, cosA, A = viz_tools.rotate_vel2(rpn_UU, rpn_VV, coords, origin='grid', coord_dict=True)

(475, 674) (475, 674)
(1, 476, 675) (1, 476, 675)
shapes (475, 674) (475, 674) (475, 674)

/Users/sallen/Documents/MEOPAR/tools/SalishSeaTools/salishsea_tools/viz_tools.py:539: RuntimeWarning: invalid value encountered in arccos
  A = np.arccos(cosA)

In [139]:

fig, axs = plt.subplots(3, 3, figsize=(12, 12))

colours = axs[0, 0].pcolormesh(coords['glamu'][0])
fig.colorbar(colours, ax=axs[0, 0])
axs[0, 0].set_title('Longitude')

colours = axs[0, 1].pcolormesh(coords['gphiu'][0])
fig.colorbar(colours, ax=axs[0, 1])
axs[0, 1].set_title('Latitude')

colours = axs[0, 2].pcolormesh(a)
fig.colorbar(colours, ax=axs[0, 2])
axs[0, 2].set_title('a')

colours = axs[1, 0].pcolormesh(b)
fig.colorbar(colours, ax=axs[1, 0])
axs[1, 0].set_title('b')

colours = axs[1, 1].pcolormesh(c)
fig.colorbar(colours, ax=axs[1, 1])
axs[1, 1].set_title('c')

colours = axs[1, 2].pcolormesh(cosA, vmax=1, vmin=0.9)
fig.colorbar(colours, ax=axs[1, 2])
axs[1, 2].set_title('cosA')

colours = axs[2, 0].pcolormesh(A, vmin=-0.2, vmax=0.2, cmap='bwr')
fig.colorbar(colours, ax=axs[2, 0])
axs[2, 0].set_title('A')

colours = axs[2, 1].pcolormesh(u_out, cmap='bwr')
fig.colorbar(colours, ax=axs[2, 1])
axs[2, 1].set_title('u_out')

colours = axs[2, 2].pcolormesh(v_out, cmap='bwr')
fig.colorbar(colours, ax=axs[2, 2])
axs[2, 2].set_title('v_out')

Out[139]:

<matplotlib.text.Text at 0x11c9be080>

In [140]:

rpn_ueast_on_ops = interpolate_rpn_to_ops(u_out[JST-1:JEN-1, IST-1:IEN-1], points, opsdata)*0.514444
rpn_vnorth_on_ops = interpolate_rpn_to_ops(v_out[JST-1:JEN-1, IST-1:IEN-1], points, opsdata)*0.514444

In [142]:

fig, axs = plt.subplots(3, 3, figsize=(18, 13.5))
colors = axs[0, 0].pcolormesh(rpn_ueast_on_ops, cmap='bwr', vmax=10, vmin=-10)
fig.colorbar(colors, ax=axs[0, 0])
opsdata['u_wind'][18].plot(ax=axs[0, 1], cmap='bwr')
colors = axs[0, 2].pcolormesh(rpn_ueast_on_ops-opsdata['u_wind'][18], cmap='bwr', vmax=1, vmin=-1)
fig.colorbar(colors, ax=axs[0, 2]);
axs[0, 0].set_title('RPN Interpolated')
axs[0, 1].set_title('OPS')
axs[0, 2].set_title('RPN Interpolated - OPS');

colors = axs[1, 0].pcolormesh(rpn_vnorth_on_ops, cmap='bwr', vmax=10, vmin=-10)
fig.colorbar(colors, ax=axs[1, 0])
opsdata['v_wind'][18].plot(ax=axs[1, 1], cmap='bwr')
colors = axs[1, 2].pcolormesh(rpn_vnorth_on_ops-opsdata['v_wind'][18], cmap='bwr', vmax=1, vmin=-1)
fig.colorbar(colors, ax=axs[1, 2]);
axs[1, 0].set_title('RPN Interpolated')
axs[1, 1].set_title('OPS')
axs[1, 2].set_title('RPN Interpolated - OPS');

rpnspeed = np.sqrt(rpn_ueast_on_ops**2 + rpn_vnorth_on_ops**2)
opsspeed = np.sqrt(opsdata['u_wind'][18]**2 + opsdata['v_wind'][18]**2)
colors = axs[2, 0].pcolormesh(rpnspeed , cmap='plasma', vmax=12)
fig.colorbar(colors, ax=axs[2, 0])
opsspeed.plot(ax=axs[2, 1], cmap='plasma', vmax=12)
colors = axs[2, 2].pcolormesh(rpnspeed - opsspeed, cmap='bwr', vmax=1, vmin=-1)
fig.colorbar(colors, ax=axs[2, 2]);
axs[1, 0].set_title('RPN Interpolated')
axs[1, 1].set_title('OPS')
axs[1, 2].set_title('RPN Interpolated - OPS');

In [97]:

plt.hist2d(np.array(opsdata['u_wind'][17]).flatten(), rpn_ueast_on_ops.flatten(), vmax=1e3);
plt.plot([-10, 10], [-10, 10], 'm')
plt.colorbar()

Out[97]:

<matplotlib.colorbar.Colorbar at 0x11ba99748>

In [ ]: