Calculate daily mean SSH at Neah Bay
import arrow
import datetime
import matplotlib.pyplot as plt
import os
import pandas as pd
import xarray as xr
%matplotlib inline
basedir = '/results/forcing/sshNeahBay/obs'
xb, yb = 50, 0
xr.open_dataset('/results/forcing/sshNeahBay/obs/ssh_y2020m12d02.nc', decode_cf=False)
<xarray.Dataset> Dimensions: (time_counter: 24, xbT: 100, yb: 1) Coordinates: * time_counter (time_counter) float32 1.0 2.0 3.0 4.0 ... 21.0 22.0 23.0 24.0 Dimensions without coordinates: xbT, yb Data variables: nav_lat (yb, xbT) float32 48.31 48.32 48.32 48.33 ... 48.7 48.7 48.71 nav_lon (yb, xbT) float32 -124.6 -124.6 -124.6 ... -125.0 -125.0 sossheig (time_counter, yb, xbT) float32 -0.2549 -0.2549 ... -0.07914 vobtcrtx (time_counter, yb, xbT) float32 0.0 0.0 0.0 ... 0.0 0.0 0.0 vobtcrty (time_counter, yb, xbT) float32 0.0 0.0 0.0 ... 0.0 0.0 0.0 nbidta (yb, xbT) int32 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0 nbjdta (yb, xbT) int32 370 371 372 373 374 ... 465 466 467 468 469 nbrdta (yb, xbT) int32 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0 Attributes: Conventions: CF-1.6 title: Neah Bay SSH hourly values institution: Dept of Earth, Ocean & Atmospheric Sciences, University of ... source: /results/forcing/sshNeahBay/txt/sshNB_2020-12-03_19.txt references: https://github.com/SalishSeaCast/SalishSeaNowcast/blob/mast... history: [2020-12-03 11:40:03] Created netCDF4 zlib=True dataset. comment: Observation from Neah Bay storm surge website generated by ...
array([ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24.], dtype=float32)
array([[48.314762, 48.318752, 48.322746, 48.326736, 48.330727, 48.334717, 48.338707, 48.342697, 48.346687, 48.350677, 48.354668, 48.358658, 48.362648, 48.36664 , 48.370625, 48.374615, 48.3786 , 48.38259 , 48.386578, 48.390568, 48.394554, 48.39854 , 48.40253 , 48.406517, 48.410503, 48.41449 , 48.418476, 48.422462, 48.42645 , 48.430435, 48.434418, 48.438404, 48.44239 , 48.446373, 48.45036 , 48.45434 , 48.45833 , 48.46231 , 48.466293, 48.47028 , 48.474262, 48.478245, 48.482227, 48.48621 , 48.490192, 48.494175, 48.498158, 48.50214 , 48.50612 , 48.5101 , 48.51408 , 48.518063, 48.52204 , 48.526024, 48.530003, 48.53398 , 48.537964, 48.541943, 48.54592 , 48.5499 , 48.55388 , 48.557858, 48.561832, 48.56581 , 48.56979 , 48.57377 , 48.577744, 48.581722, 48.585697, 48.589672, 48.59365 , 48.597626, 48.6016 , 48.605576, 48.60955 , 48.613525, 48.6175 , 48.621475, 48.62545 , 48.629425, 48.633396, 48.63737 , 48.641346, 48.645317, 48.64929 , 48.653263, 48.657234, 48.661205, 48.665176, 48.669147, 48.67312 , 48.67709 , 48.68106 , 48.68503 , 48.689003, 48.692974, 48.69694 , 48.700912, 48.70488 , 48.70885 ]], dtype=float32)
array([[-124.628914, -124.63225 , -124.63558 , -124.63892 , -124.64226 , -124.6456 , -124.64893 , -124.65227 , -124.65561 , -124.65894 , -124.662285, -124.66562 , -124.66896 , -124.672295, -124.67564 , -124.67897 , -124.68231 , -124.685646, -124.68899 , -124.69233 , -124.69566 , -124.699005, -124.70234 , -124.70568 , -124.70902 , -124.71236 , -124.7157 , -124.71904 , -124.722374, -124.725716, -124.72906 , -124.73239 , -124.73573 , -124.739075, -124.74242 , -124.74575 , -124.74909 , -124.75243 , -124.755775, -124.75911 , -124.76245 , -124.76579 , -124.769135, -124.77248 , -124.77581 , -124.77915 , -124.78249 , -124.785835, -124.78918 , -124.79252 , -124.79586 , -124.799194, -124.802536, -124.80588 , -124.80922 , -124.81256 , -124.8159 , -124.819244, -124.822586, -124.82593 , -124.82927 , -124.83261 , -124.83595 , -124.839294, -124.842636, -124.84598 , -124.84932 , -124.85266 , -124.856 , -124.859344, -124.862686, -124.86603 , -124.86937 , -124.87271 , -124.87605 , -124.879395, -124.88274 , -124.88608 , -124.88942 , -124.89277 , -124.89611 , -124.89945 , -124.902794, -124.906136, -124.90948 , -124.91282 , -124.91616 , -124.91951 , -124.92285 , -124.92619 , -124.929535, -124.93288 , -124.936226, -124.93957 , -124.94291 , -124.94625 , -124.94959 , -124.95294 , -124.95628 , -124.959625]], dtype=float32)
array([[[-0.25489 , -0.25489 , ..., -0.25489 , -0.25489 ]], [[-0.26432 , -0.26432 , ..., -0.26432 , -0.26432 ]], ..., [[-0.090255, -0.090255, ..., -0.090255, -0.090255]], [[-0.079137, -0.079137, ..., -0.079137, -0.079137]]], dtype=float32)
array([[[0., 0., ..., 0., 0.]], [[0., 0., ..., 0., 0.]], ..., [[0., 0., ..., 0., 0.]], [[0., 0., ..., 0., 0.]]], dtype=float32)
array([[[0., 0., ..., 0., 0.]], [[0., 0., ..., 0., 0.]], ..., [[0., 0., ..., 0., 0.]], [[0., 0., ..., 0., 0.]]], dtype=float32)
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=int32)
array([[370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469]], dtype=int32)
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=int32)
year = 2018
start = datetime.datetime(year, 1, 1)
endtime = datetime.datetime(year, 12, 31)
timerange = arrow.Arrow.range('day', start, endtime)
for i, day in enumerate(timerange):
ymd = day.format('YYYYMMDD')
filename = f'ssh_y{day.year}m{day.month:02d}d{day.day:02d}.nc'
fullfile = os.path.join(basedir, filename)
neah = xr.open_dataset(fullfile, decode_cf=False)
ssh = neah['sossheig'].isel(yb=yb, xbT=xb).mean()
neah.close()
if i == 0:
ssh_year = ssh.copy(deep=True)
ssh.close()
else:
ssh_year = xr.concat([ssh_year, ssh], dim='time_counter')
ssh.close()
if i % 100 == 0:
print (i)
0 100 200 300
ssh_year.plot();
ssh_year.to_netcdf(f'ssh_{year}.nc')
ssh2015 = xr.open_dataset('ssh_2015.nc')
ssh2016 = xr.open_dataset('ssh_2016.nc')
ssh2017 = xr.open_dataset('ssh_2017.nc')
ssh2018 = xr.open_dataset('ssh_2018.nc')
ssh = xr.concat([ssh2015, ssh2016, ssh2017, ssh2018], dim='time_counter', coords='minimal')
ssh.sossheig.plot();
ssh_pd = ssh.to_dataframe()
ssh_pd = ssh_pd.drop('xbT', 1)
ssh_pd = ssh_pd.drop('yb', 1)
ssh_pd.to_csv('day_avg_ssh.csv')
ssh_pd
sossheig | |
---|---|
time_counter | |
0 | -0.125347 |
1 | -0.038566 |
2 | -0.063083 |
3 | 0.001806 |
4 | 0.104299 |
... | ... |
1456 | -0.048620 |
1457 | -0.078037 |
1458 | 0.019949 |
1459 | -0.058769 |
1460 | -0.205427 |
1461 rows × 1 columns
low_pass_ssh = ssh_pd.rolling(4, center=True).mean()
low_pass_ssh.to_csv('low_pass_ssh.csv')
low_pass_ssh.plot()
<AxesSubplot:xlabel='time_counter'>
ssh2015.close()
ssh2016.close()
ssh2017.close()
ssh2018.close()