Sandy Cove Analysis Aug 19, 2017
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
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))
filename = 'SandyCove'
place = 'Sandy Cove'
mean_sea_lvl = places.PLACES[place]['mean sea lvl']
start = datetime(2017, 8, 21)
end = datetime(2018, 5, 8)
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')
# Save just in case
ssh.to_netcdf('sandycove.nc')
ssh.sossheig.plot();
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:
# Save just in case
obs.to_netcdf('SCobc.nc')
# Read it in for speed
obs = xr.open_dataset('SCobc.nc')
ssh.sossheig.plot()
(obs.water_level - mean_sea_lvl).plot();
#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
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')
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)
# Save just in case
ttide_ds.to_netcdf('extract_pred.nc')
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))
(0, 200)
ssh_correction = ttide_ds.pred_noshallow - ttide_ds.pred_8
ssh_corrected = np.array(ssh.sossheig[:, 0, 0].values) + np.array(ssh_correction.values)
#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
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')
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)
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
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')