Sandy Cove Analysis Aug 19, 2017

In [2]:
import arrow
from datetime import datetime, timedelta
import matplotlib
from matplotlib.colors import LogNorm
import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
import pandas as pd
import xarray as xr

from salishsea_tools import data_tools, places
import shared

%matplotlib inline

Statistics Functions

In [3]:
def rmse(predictions, targets):
    return np.sqrt(np.nanmean((predictions - targets) ** 2))
def bias(predictions, targets):
    return np.nanmean(predictions - targets)
def willmott(predictions, targets):
    xbar = np.nanmean(targets)
    return (1 - (np.nansum((predictions - targets)**2))  / 
                             np.nansum((np.abs(predictions - xbar)
                                     + np.abs(targets - xbar))**2))

Set-up Location and Dates

In [4]:
filename = 'SandyCove'
place = 'Sandy Cove'
mean_sea_lvl = places.PLACES[place]['mean sea lvl']
In [5]:
start = datetime(2017, 8, 21)
end = datetime(2018, 5, 8)

Get Model Results

In [6]:
ssh = xr.open_dataset('/results/SalishSea/nowcast-green/20aug17/'+filename+'.nc')
for aday in arrow.Arrow.range('day', start, end):
    day = aday.format('DDMMMYY').lower()
    newdata = xr.open_dataset('/results/SalishSea/nowcast-green/'+day+'/'+filename+'.nc')
    ssh = xr.concat([ssh, newdata], dim='time_counter')
In [7]:
# Save just in case
ssh.to_netcdf('sandycove.nc')
In [9]:
ssh.sossheig.plot();

Get Observations

In [10]:
obs_1min = data_tools.get_chs_tides(
            'obs',
            place,
            arrow.get(str(ssh.time_counter[0].values)) -
            timedelta(seconds=5 * 60),
            arrow.get(str(ssh.time_counter[-1].values))
        )
obs_10min_avg = xr.DataArray(
            obs_1min.resample('10min', loffset='5min').mean()
        )
obs = xr.Dataset({
            'water_level': obs_10min_avg.rename({
                'dim_0': 'time'
            })
        })
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-10-9ef6b0405ed4> in <module>()
      4             arrow.get(str(ssh.time_counter[0].values)) -
      5             timedelta(seconds=5 * 60),
----> 6             arrow.get(str(ssh.time_counter[-1].values))
      7         )
      8 obs_10min_avg = xr.DataArray(

/ocean/sallen/allen/research/Meopar/Tools/SalishSeaTools/salishsea_tools/data_tools.py in get_chs_tides(data_type, stn_id, begin, end)
    486             end.format('YYYY-MM-DD HH:mm:ss'),
    487             first_value, n_values, metadata, 'station_id={}'.format(stn_number),
--> 488             order
    489         )
    490         response = zeep.helpers.serialize_object(response)

/home/sallen/anaconda/envs/py36/lib/python3.6/site-packages/zeep/client.py in __call__(self, *args, **kwargs)
     43         return self._proxy._binding.send(
     44             self._proxy._client, self._proxy._binding_options,
---> 45             self._op_name, args, kwargs)
     46 
     47 

/home/sallen/anaconda/envs/py36/lib/python3.6/site-packages/zeep/wsdl/bindings/soap.py in send(self, client, options, operation, args, kwargs)
    111 
    112         response = client.transport.post_xml(
--> 113             options['address'], envelope, http_headers)
    114 
    115         operation_obj = self.get(operation)

/home/sallen/anaconda/envs/py36/lib/python3.6/site-packages/zeep/transports.py in post_xml(self, address, envelope, headers)
     93         """
     94         message = etree_to_string(envelope)
---> 95         return self.post(address, message, headers)
     96 
     97     def load(self, url):

/home/sallen/anaconda/envs/py36/lib/python3.6/site-packages/zeep/transports.py in post(self, address, message, headers)
     65             data=message,
     66             headers=headers,
---> 67             timeout=self.operation_timeout)
     68 
     69         if self.logger.isEnabledFor(logging.DEBUG):

/home/sallen/anaconda/envs/py36/lib/python3.6/site-packages/requests/sessions.py in post(self, url, data, json, **kwargs)
    533         """
    534 
--> 535         return self.request('POST', url, data=data, json=json, **kwargs)
    536 
    537     def put(self, url, data=None, **kwargs):

/home/sallen/anaconda/envs/py36/lib/python3.6/site-packages/requests/sessions.py in request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)
    486         }
    487         send_kwargs.update(settings)
--> 488         resp = self.send(prep, **send_kwargs)
    489 
    490         return resp

/home/sallen/anaconda/envs/py36/lib/python3.6/site-packages/requests/sessions.py in send(self, request, **kwargs)
    639 
    640         if not stream:
--> 641             r.content
    642 
    643         return r

/home/sallen/anaconda/envs/py36/lib/python3.6/site-packages/requests/models.py in content(self)
    779                 self._content = None
    780             else:
--> 781                 self._content = bytes().join(self.iter_content(CONTENT_CHUNK_SIZE)) or bytes()
    782 
    783         self._content_consumed = True

/home/sallen/anaconda/envs/py36/lib/python3.6/site-packages/requests/models.py in generate()
    701             if hasattr(self.raw, 'stream'):
    702                 try:
--> 703                     for chunk in self.raw.stream(chunk_size, decode_content=True):
    704                         yield chunk
    705                 except ProtocolError as e:

/home/sallen/anaconda/envs/py36/lib/python3.6/site-packages/requests/packages/urllib3/response.py in stream(self, amt, decode_content)
    426         """
    427         if self.chunked and self.supports_chunked_reads():
--> 428             for line in self.read_chunked(amt, decode_content=decode_content):
    429                 yield line
    430         else:

/home/sallen/anaconda/envs/py36/lib/python3.6/site-packages/requests/packages/urllib3/response.py in read_chunked(self, amt, decode_content)
    591                 if self.chunk_left == 0:
    592                     break
--> 593                 chunk = self._handle_chunk(amt)
    594                 decoded = self._decode(chunk, decode_content=decode_content,
    595                                        flush_decoder=False)

/home/sallen/anaconda/envs/py36/lib/python3.6/site-packages/requests/packages/urllib3/response.py in _handle_chunk(self, amt)
    556             returned_chunk = value
    557         else:  # amt > self.chunk_left
--> 558             returned_chunk = self._fp._safe_read(self.chunk_left)
    559             self._fp._safe_read(2)  # Toss the CRLF at the end of the chunk.
    560             self.chunk_left = None

/home/sallen/anaconda/envs/py36/lib/python3.6/http/client.py in _safe_read(self, amt)
    610         s = []
    611         while amt > 0:
--> 612             chunk = self.fp.read(min(amt, MAXAMOUNT))
    613             if not chunk:
    614                 raise IncompleteRead(b''.join(s), amt)

/home/sallen/anaconda/envs/py36/lib/python3.6/socket.py in readinto(self, b)
    584         while True:
    585             try:
--> 586                 return self._sock.recv_into(b)
    587             except timeout:
    588                 self._timeout_occurred = True

/home/sallen/anaconda/envs/py36/lib/python3.6/site-packages/requests/packages/urllib3/contrib/pyopenssl.py in recv_into(self, *args, **kwargs)
    254     def recv_into(self, *args, **kwargs):
    255         try:
--> 256             return self.connection.recv_into(*args, **kwargs)
    257         except OpenSSL.SSL.SysCallError as e:
    258             if self.suppress_ragged_eofs and e.args == (-1, 'Unexpected EOF'):

/home/sallen/anaconda/envs/py36/lib/python3.6/site-packages/OpenSSL/SSL.py in recv_into(self, buffer, nbytes, flags)
   1332             result = _lib.SSL_peek(self._ssl, buf, nbytes)
   1333         else:
-> 1334             result = _lib.SSL_read(self._ssl, buf, nbytes)
   1335         self._raise_ssl_error(self._ssl, result)
   1336 

KeyboardInterrupt: 
In [ ]:
# Save just in case
obs.to_netcdf('SCobc.nc')
In [11]:
# Read it in for speed
obs = xr.open_dataset('SCobc.nc')
In [15]:
ssh.sossheig.plot()
(obs.water_level - mean_sea_lvl).plot();

Base Comparison

In [46]:
#plt.plot(obs.water_level, ssh.sossheig.data[:, 0, 0], 'o')
print('RMSE', rmse(ssh.sossheig.data[:, 0, 0], obs.water_level - mean_sea_lvl))
print('bias', bias(ssh.sossheig.data[:, 0, 0], obs.water_level - mean_sea_lvl))
print('Willmott', willmott(ssh.sossheig.data[:, 0, 0], obs.water_level - mean_sea_lvl))
RMSE 0.105347847994
bias 0.0307859186572
Willmott 0.996905645365
In [47]:
ssh_obs = np.array(obs.water_level) - mean_sea_lvl
ssh_model = np.array(ssh.sossheig.data[:, 0, 0])
matplotlib.rcParams.update({'font.size': 16})
fig, ax = plt.subplots(1, 1, figsize = (8.5, 7))
c, xedge, yedge, im = ax.hist2d(ssh_obs[~np.isnan(ssh_obs)],
                                ssh_model[~np.isnan(ssh_obs)], bins = 100, norm=LogNorm(), cmap='copper')
im
ax.plot([-3.2, 2.2], [-3.2, 2.2], 'w')
ax.set_xlim((-3.2, 2.2))
ax.set_ylim((-3.2, 2.2))
ax.set_xlabel('Observations (m)')
ax.set_ylabel('Raw Model (m)')
ax.set_title('Sea Surface Height at Sandy Cove,\n Aug 20, 2017- May 8, 2018')
cb = fig.colorbar(im, ax=ax)
cb.set_label('Point Density')

Get tidal Predictions

In [21]:
tidal_predictions = '/results/nowcast-sys/SalishSeaNowcast/tidal_predictions/'
model_ssh_period = slice(
        str(ssh.time_centered[0].values), 
        str(ssh.time_centered[-1].values))

# Predicted tide water levels dataset from ttide
ttide = shared.get_tides(place, tidal_predictions)
ttide.rename(columns={' pred_noshallow ': 'pred_noshallow'}, inplace=True)
ttide.index = ttide.time
ttide_full = xr.Dataset.from_dataframe(ttide)
ttide_ds = ttide_full.sel(time=model_ssh_period)
In [ ]:
# Save just in case
ttide_ds.to_netcdf('extract_pred.nc')
In [22]:
fig, ax = plt.subplots(1, 1, figsize=(20, 4))
ax.plot(ssh.sossheig[:, 0, 0])
ax.plot(obs.water_level - mean_sea_lvl)
ax.plot(ttide_ds.pred_all)
ax.set_xlim((0, 200))
Out[22]:
(0, 200)

Correct the Model

In [23]:
ssh_correction = ttide_ds.pred_noshallow - ttide_ds.pred_8
ssh_corrected = np.array(ssh.sossheig[:, 0, 0].values) + np.array(ssh_correction.values)
In [45]:
#plt.plot(obs.water_level, ssh_corrected, 'o')
print('RMSE', rmse(ssh_corrected, obs.water_level - mean_sea_lvl))
print('bias', bias(ssh_corrected, obs.water_level - mean_sea_lvl))
print('Willmott', willmott(ssh_corrected, obs.water_level - mean_sea_lvl))
RMSE 0.0834871724631
bias 0.0309168653373
Willmott 0.998085656952
In [44]:
ssh_obs = np.array(obs.water_level) - mean_sea_lvl
matplotlib.rcParams.update({'font.size': 16})
fig, ax = plt.subplots(1, 1, figsize = (8.5, 7))
c, xedge, yedge, im = ax.hist2d(ssh_obs[~np.isnan(ssh_obs)],
                                ssh_corrected[~np.isnan(ssh_obs)], bins = 100, norm=LogNorm(), cmap='copper')
im
ax.plot([-3.2, 2.2], [-3.2, 2.2], 'w')
ax.set_xlim((-3.2, 2.2))
ax.set_ylim((-3.2, 2.2))
ax.set_xlabel('Observations (m)')
ax.set_ylabel('Corrected Model (m)')
ax.set_title('Sea Surface Height at Sandy Cove,\n Aug 20, 2017- May 8, 2018')
cb = fig.colorbar(im, ax=ax)
cb.set_label('Point Density')

Look at Residuals

In [25]:
residual_model = ssh_corrected - np.array(ttide_ds.pred_all.values)
residual_obs = np.array(obs.water_level) - mean_sea_lvl - np.array(ttide_ds.pred_all.values)
In [26]:
print('RMSE', rmse(residual_model, residual_obs))
print('bias', bias(residual_model, residual_obs))
print('Willmott', willmott(residual_model, residual_obs))
RMSE 0.0834871724631
bias 0.0309168653373
Willmott 0.947058353941
In [27]:
matplotlib.rcParams.update({'font.size': 16})
fig, ax = plt.subplots(1, 1, figsize = (8.5, 7))
c, xedge, yedge, im = ax.hist2d(residual_obs[~np.isnan(residual_obs)], 
                                residual_model[~np.isnan(residual_obs)], bins = 100, norm=LogNorm(), cmap='copper')
im
ax.plot([-0.48, 0.72], [-0.48, 0.72], 'w')
ax.set_xlim((-0.48, 0.72))
ax.set_ylim((-0.48, 0.72))
ax.set_xlabel('Observations (m)')
ax.set_ylabel('Corrected Model (m)')
ax.set_title('Tidal Residuals at Sandy Cove,\n Aug 20, 2017- May 8, 2018')
cb = fig.colorbar(im, ax=ax)
cb.set_label('Point Density')
In [ ]: