<jupyter_start><jupyter_text>Notebook to summarize IOS Fall 2014 comparisons. Profile comaprisons condensed to a scatter plot.<jupyter_code>import os
import glob

import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import datetime

import comparisons

from salishsea_tools import tidetools, viz_tools, teos_tools
from nowcast import analyze

import seaborn as sns

import ACTDR
import pandas as pd
%matplotlib inline
import math
sns.set_color_codes()
sns.set_context('talk')
nowcast_dir = '/results/SalishSea/nowcast/'
grid_B = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')    
mesh = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/mesh_mask_SalishSea2.nc')
FIRST_NOWCAST = datetime.datetime(2014, 10,27)
def get_model_cast(year, month, day, lon, lat, field):
    """Get the model daily mean, max and min corresponding to given year, month, day and lat/lon"""
    
    if field =='Salinity':
        var = 'vosaline'
    elif field == 'Temperature':
        var = 'votemper'
        
    date = datetime.datetime(year, month, day)
    early_days = False
    if date < FIRST_NOWCAST:
        early_days = True
    #Look up grid point
    bathy, X, Y= tidetools.get_bathy_data(grid_B)
    j, i = tidetools.find_closest_model_point(lon, lat, X, Y, bathy)
    # Grab data
    depths = mesh.variables['gdept'][0,:,j,i]
    if early_days:
        var_model = comparisons.early_model_data(date, j, i, '1h', 'grid_T', var, nowcast_dir)
    else:
        results_dir = os.path.join(nowcast_dir, date.strftime('%d%b%y').lower())
        grid_T_h = comparisons.results_dataset('1h', 'grid_T', results_dir)
        var_model = grid_T_h.variables[var][:, :, j, i]
    var_model, tmask = mask_model(j, i, var_model)
    # Daily mean/min/max
    mean_daily = np.mean(var_model, axis=0)
    max_daily = np.max(var_model, axis=0)
    min_daily = np.min(var_model, axis=0)
 
    return mean_daily, max_daily, min_daily, depths, tmask
def mask_model(j, i, var_model):
    """Mask a model variable at a grid point"""
    
    tmask = np.abs(1-mesh.variables['tmask'][:, :, j, i])
    tmask = tmask + np.zeros(var_model.shape)
    var_masked = np.ma.array(var_model, mask=tmask)
    return var_masked, tmask[0,:]
def prepare_cast_comparison(dep_obs, var_obs, year, month, day, lon, lat, field):
    """Gather model and observed cast data for comparison"""
    
    # model data 
    mean_mod, max_mod, min_mod, dep_mod, tmask = get_model_cast(year, month, day, lon, lat, field)
    
    # interp observations
    interp_obs = comparisons.interpolate_depth(var_obs, dep_obs, dep_mod)
    interp_obs = np.ma.array(interp_obs, mask=tmask)
    
    return interp_obs, mean_mod, max_mod, min_mod, dep_mod
ACTDR.load_dat('SOG_2000.dat')
data = pd.DataFrame(ACTDR.CTD_DAT)
data_fall = data[(data.Year>=2014) & (data.Month >=10)  ]
sns.set_context('talk')
def scatter_compare(data, field, cmin, cmax, ax, units):
    sal_flags=[]
    for dep_obs, var_obs, lon, lat, year, month, day in zip(data['Depth'],
                                                        data[field],
                                                        data['Longitude'],
                                                        data['Latitude'],
                                                        data['Year'],
                                                        data['Month'],
                                                        data['Day']):
        interp_obs, mean_mod, max_mod, min_mod, dep_mod = prepare_cast_comparison(dep_obs, var_obs,
                                                                                  year, month, day,
                                                                                  lon, lat, field)
        thres=1; date=datetime.datetime(year, month, day)
        if field == 'Salinity':
            flag, k=classify_surface_salinity(interp_obs, mean_mod, max_mod, min_mod, thres)
            if flag == 0:
                print('Salinity within {} at depth {}, {}, {}, {}'.format(thres, dep_mod[k], lon, lat, date))
            if flag == 1:
                print('Too salty at depth {}, {}, {}, {}'.format(dep_mod[k], lon, lat, date))
            if flag == -1:
                print('Too fresh at depth {}, {}, {}, {}'.format(dep_mod[k], lon, lat, date))
            sal_flags.append(flag)
        ax.errorbar(interp_obs, mean_mod, yerr=[mean_mod-min_mod,max_mod-mean_mod ],fmt='k:', marker='',
                     zorder=0, ecolor='gray',lw=1,)
        sc = ax.scatter(interp_obs, mean_mod, edgecolor='',c= dep_mod, 
                         norm=mcolors.LogNorm(), cmap='Spectral', vmin=0.5, vmax=300, s=50)
        ax.set_xlim([cmin, cmax])
        ax.set_ylim([cmin, cmax])
        ax.plot([cmin, cmax],[cmin, cmax], 'r-',zorder=0)
    cb = plt.colorbar(sc, ax=ax)
    ticks = [1, 10, 25, 50, 100, 200, 400]
    cb.set_ticks(ticks)
    cb.set_ticklabels(ticks)
    cb.set_label('Depth [m]')
    ax.set_xlabel('Observed {} [{}]'.format(field, units))
    ax.set_ylabel('Modelled {} [{}]'.format(field, units))
    return sal_flags
def compare_region(data, lon_min, lon_max, lat_min, lat_max, title, smin=18, smax=32, tmin=8, tmax=16):
    cmap = sns.diverging_palette(220, 20,as_cmap=True)
    data_region  = comparisons.isolate_region(data,lon_min, lon_max, lat_min, lat_max)
    #sal
    fig2,ax2 = plt.subplots(1,1)
    field = 'Salinity'; 
    sal_flags=scatter_compare(data_region, field, smin,smax,ax2,'g/kg')
    ax2.set_title('{} - {}'.format(title, field))
    #temp
    fig3,ax3 = plt.subplots(1,1)
    field = 'Temperature';
    scatter_compare(data_region, field, tmin,tmax,ax3,'deg C')
    ax3.set_title('{} - {}'.format(title, field))
    #map
    fig1, ax1 = plt.subplots(1,1)
    ax1.scatter(data_region['Longitude'], data_region['Latitude'], s=50,c=sal_flags, cmap=cmap,vmin=-1,vmax=1)
    viz_tools.plot_coastline(ax1,grid_B, coords='map')
    ax1.set_xlim([lon_min, lon_max])
    ax1.set_ylim([lat_min, lat_max])
    ax1.set_title(title)
   

    return fig1, ax1
def classify_surface_salinity(interp_obs, mean_mod, max_mod, min_mod, thres):

    interp_obs = np.ma.masked_invalid(interp_obs)
    nonmasked = np.ma.flatnotmasked_edges(interp_obs)
    k=nonmasked[0]
    if (min_mod[k] -interp_obs[k] >thres) :
        flag = 1 #salty surface
    elif (max_mod[k] - interp_obs[k] )<-thres:
        flag = -1 #fresh surface
    else:
        flag =0  #ok
    
    return flag, k<jupyter_output><empty_output><jupyter_text>Strait of Georgia<jupyter_code>lon_min=-124
lon_max=-123
lat_min=48.8
lat_max=49.35
fig1,ax1 = compare_region(data_fall, lon_min, lon_max, lat_min, lat_max, 'Southern Strait of Georgia')
ax1.plot(49.2611, -123.2531, 'r*')<jupyter_output>Too fresh at depth 3.50003051758, -123.527, 49.1026666667, 2014-10-01 00:00:00
Too fresh at depth 3.50003051758, -123.310333333, 49.1673333333, 2014-10-01 00:00:00
Salinity within 1 at depth 3.50003051758, -123.429833333, 49.2795, 2014-10-02 00:00:00
Salinity within 1 at depth 5.50015068054, -123.7185, 49.2735, 2014-10-02 00:00:00
Salinity within 1 at depth 3.50003051758, -123.082333333, 48.8425, 2014-10-03 00:00:00
Salinity within 1 at depth 3.50003051758, -123.62, 49.0896666667, 2014-10-03 00:00:00
Too fresh at depth 1.50000309944, -123.568166667, 49.201, 2014-10-02 00:00:00
Too fresh at depth 2.50001144409, -123.303833333, 48.9031666667, 2014-10-03 00:00:00
Too salty at depth 2.50001144409, -123.180333333, 48.8565, 2014-10-28 00:00:00
Too salty at depth 3.50003051758, -123.4385, 49.03, 2014-10-28 00:00:00
Too salty at depth 2.50001144409, -123.372666667, 49.0546666667, 2014-10-28 00:00:00
Too salty at depth 1.50000309944, -123.551166667, 49.1643333333, 2014-10-28 00:00:00
Too salty [...]<jupyter_text>Northern SoG<jupyter_code>lon_min=-125.5
lon_max=-123
lat_min=49.35
lat_max=51
compare_region(data_fall, lon_min, lon_max, lat_min, lat_max, 'Northern Strait of Georgia',
               smin=18,smax=31.5, tmin=7, tmax=16,
               )<jupyter_output>Salinity within 1 at depth 5.50015068054, -123.368666667, 49.422, 2014-10-01 00:00:00
Too fresh at depth 3.50003051758, -124.739833333, 49.7013333333, 2014-10-04 00:00:00
Salinity within 1 at depth 3.50003051758, -124.675, 49.6345, 2014-10-04 00:00:00
Salinity within 1 at depth 3.50003051758, -124.472333333, 49.5313333333, 2014-10-05 00:00:00
Salinity within 1 at depth 3.50003051758, -124.545666667, 49.4931666667, 2014-10-05 00:00:00
Salinity within 1 at depth 3.50003051758, -124.446166667, 49.381, 2014-10-05 00:00:00
Salinity within 1 at depth 3.50003051758, -124.143333333, 49.3916666667, 2014-10-06 00:00:00
Too fresh at depth 3.50003051758, -123.238666667, 49.6, 2014-10-07 00:00:00
Salinity within 1 at depth 4.5000705719, -123.622833333, 49.3938333333, 2014-10-07 00:00:00
Too fresh at depth 3.50003051758, -124.021666667, 49.799, 2014-10-08 00:00:00
Salinity within 1 at depth 4.5000705719, -124.233166667, 49.7033333333, 2014-10-08 00:00:00
Too salty at depth 4.5000705719, -124.2733333[...]<jupyter_text>Strait of Juan de Fuca<jupyter_code>lon_min=-125
lon_max=-123
lat_min=48
lat_max=48.48
compare_region(data_fall, lon_min, lon_max, lat_min, lat_max, 'Strait of Juan de Fuca',
              smin=30,smax=34, tmin=7, tmax=15)<jupyter_output>Too salty at depth 2.50001144409, -124.067166667, 48.3101666667, 2014-10-27 00:00:00
Too salty at depth 2.50001144409, -123.716666667, 48.255, 2014-10-27 00:00:00
Salinity within 1 at depth 2.50001144409, -123.300666667, 48.2331666667, 2014-10-27 00:00:00
Salinity within 1 at depth 2.50001144409, -123.1655, 48.2651666667, 2014-10-27 00:00:00
Salinity within 1 at depth 3.50003051758, -123.043166667, 48.3806666667, 2014-10-28 00:00:00
Salinity within 1 at depth 1.50000309944, -123.664, 48.29, 2014-10-16 00:00:00
Salinity within 1 at depth 4.5000705719, -124.081166667, 48.3896666667, 2014-10-16 00:00:00
Salinity within 1 at depth 4.5000705719, -124.341, 48.4721666667, 2014-10-16 00:00:00<jupyter_text>Haro<jupyter_code>lon_min=-123.49
lon_max=-122
lat_min=48.48
lat_max=48.8
fig1,ax1=compare_region(data_fall, lon_min, lon_max, lat_min, lat_max, 'Haro Strait',
              smin=27,smax=32, tmin=9, tmax=13)<jupyter_output>Salinity within 1 at depth 1.50000309944, -123.243666667, 48.6295, 2014-10-27 00:00:00
Salinity within 1 at depth 2.50001144409, -123.213166667, 48.5645, 2014-10-27 00:00:00
Salinity within 1 at depth 2.50001144409, -123.1545, 48.4863333333, 2014-10-27 00:00:00
Salinity within 1 at depth 3.50003051758, -123.027, 48.772, 2014-10-28 00:00:00
{'fraction_non_alphanumeric': 0.13188474256022673, 'fraction_numerical': 0.1721303731695796, 'mean_word_length': 3.2158502588610114, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 7, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '36803', 'license': 'Apache-2.0', 'n_tokens_mistral': 5255, 'n_tokens_neox': 4014, 'n_words': 894, 'path': 'Nancy/strat/comparisons/Summarizing IOS comparisons.ipynb', 'repo_name': 'SalishSeaCast/analysis', 'size': '10585'}