Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WPTO Hindcast #42

Closed
wants to merge 13 commits into from
179 changes: 179 additions & 0 deletions examples/WPTO_hindcast_example.html

Large diffs are not rendered by default.

Binary file added examples/WPTO_hindcast_example.mlx
Binary file not shown.
Binary file modified mhkit.mltbx
Binary file not shown.
80 changes: 80 additions & 0 deletions mhkit/tests/Wave_TestIO.m
Original file line number Diff line number Diff line change
Expand Up @@ -76,5 +76,85 @@ function test_ndbc_request_data(testCase)


end

function test_WPTO_point_multiyear(testCase)

hindcast_data = request_wpto_point_data('3-hour',...
{'significant_wave_height'},{44.624076,-124.280097},...
{1990,1991});
file = '../../examples/data/wave/hindcast/multi_year_hindcast.csv';
meta = '../../examples/data/wave/hindcast/multi_year_meta.csv';
expected_data = readtable(file,'delimiter',',');
expected_meta = readtable(meta);
expected_data.time_index = datetime(expected_data.time_index,'InputFormat','yyyy-MM-dd HH:mm:ssXXX',...
'TimeZone','UTC');

%
assertEqual(testCase,expected_data.time_index,hindcast_data.time.');
assertEqual(testCase,expected_data.significant_wave_height_0,hindcast_data.location0.significant_wave_height,'RelTol',0.000001);
assertEqual(testCase,expected_meta.latitude,hindcast_data.location0.latitude,'RelTol',0.000001);
assertEqual(testCase,expected_meta.longitude,hindcast_data.location0.longitude,'RelTol',0.000001);
assertEqual(testCase,expected_meta.water_depth,hindcast_data.location0.water_depth,'RelTol',0.000001);
assertEqual(testCase,expected_meta.timezone,hindcast_data.location0.timezone);
assertEqual(testCase,string(expected_meta.jurisdiction),hindcast_data.location0.jurisdiction);
assertEqual(testCase,expected_meta.distance_to_shore,hindcast_data.location0.distance_to_shore,'RelTol',0.000001);

end

function test_WPTO_point_multiloc(testCase)

hindcast_data = request_wpto_point_data('3-hour',...
{'mean_absolute_period'},{{44.624076,-124.280097},{43.489171,-125.152137}},...
{1995});
file = '../../examples/data/wave/hindcast/single_year_hindcast_multiloc.csv';
meta = '../../examples/data/wave/hindcast/multiloc_meta.csv';
expected_data = readtable(file,'delimiter',',');
expected_meta = readtable(meta);
expected_data.time_index = datetime(expected_data.time_index,'InputFormat','yyyy-MM-dd HH:mm:ssXXX',...
'TimeZone','UTC');

%
assertEqual(testCase,expected_data.time_index,hindcast_data.time.');
assertEqual(testCase,expected_data.mean_absolute_period_0,hindcast_data.location0.mean_absolute_period,'RelTol',0.000001);
assertEqual(testCase,expected_meta.latitude(1),hindcast_data.location0.latitude,'RelTol',0.000001);
assertEqual(testCase,expected_meta.longitude(1),hindcast_data.location0.longitude,'RelTol',0.000001);
assertEqual(testCase,expected_meta.water_depth(1),hindcast_data.location0.water_depth,'RelTol',0.000001);
assertEqual(testCase,expected_meta.timezone(1),hindcast_data.location0.timezone);
assertEqual(testCase,string(expected_meta.jurisdiction(1)),hindcast_data.location0.jurisdiction);
assertEqual(testCase,expected_meta.distance_to_shore(1),hindcast_data.location0.distance_to_shore,'RelTol',0.000001);

assertEqual(testCase,expected_data.mean_absolute_period_1,hindcast_data.location1.mean_absolute_period,'RelTol',0.000001);
assertEqual(testCase,expected_meta.latitude(2),hindcast_data.location1.latitude,'RelTol',0.000001);
assertEqual(testCase,expected_meta.longitude(2),hindcast_data.location1.longitude,'RelTol',0.000001);
assertEqual(testCase,expected_meta.water_depth(2),hindcast_data.location1.water_depth,'RelTol',0.000001);
assertEqual(testCase,expected_meta.timezone(2),hindcast_data.location1.timezone);
assertEqual(testCase,string(expected_meta.jurisdiction(2)),hindcast_data.location1.jurisdiction);
assertEqual(testCase,expected_meta.distance_to_shore(2),hindcast_data.location1.distance_to_shore,'RelTol',0.000001);

end

function test_WPTO_point_multiparm(testCase)

hindcast_data = request_wpto_point_data('1-hour',...
{'energy_period','mean_zero-crossing_period'},{44.624076,-124.280097},...
{1996});
file = '../../examples/data/wave/hindcast/multiparm.csv';
meta = '../../examples/data/wave/hindcast/multiparm_meta.csv';
expected_data = readtable(file,'delimiter',',');
expected_meta = readtable(meta);
expected_data.time_index = datetime(expected_data.time_index,'InputFormat','yyyy-MM-dd HH:mm:ssXXX',...
'TimeZone','UTC');

%
assertEqual(testCase,expected_data.time_index,hindcast_data.time.');
assertEqual(testCase,expected_data.energy_period_0,hindcast_data.location0.energy_period,'RelTol',0.000001);
assertEqual(testCase,expected_data.mean_zero_crossing_period_0,hindcast_data.location0.mean_zero_crossing_period,'RelTol',0.000001);
assertEqual(testCase,expected_meta.latitude(1),hindcast_data.location0.latitude,'RelTol',0.000001);
assertEqual(testCase,expected_meta.longitude(1),hindcast_data.location0.longitude,'RelTol',0.000001);
assertEqual(testCase,expected_meta.water_depth(1),hindcast_data.location0.water_depth,'RelTol',0.000001);
assertEqual(testCase,expected_meta.timezone(1),hindcast_data.location0.timezone);
assertEqual(testCase,string(expected_meta.jurisdiction(1)),hindcast_data.location0.jurisdiction);
assertEqual(testCase,expected_meta.distance_to_shore(1),hindcast_data.location0.distance_to_shore,'RelTol',0.000001);
end
end
end
178 changes: 178 additions & 0 deletions mhkit/wave/IO/hindcast/request_wpto_point_data.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,178 @@
function datast=request_wpto_point_data(data_type,parameter,lat_lon, years,options)

%%%%%%%%%%%%%%%%%%%%
% Returns data from the WPTO wave hindcast hosted on AWS at the specified
% latitude and longitude points. Visit https://registry.opendata.aws/wpto-pds-us-wave/
% for more information about the dataset and available
% locations and years.
%
% Note: To access the WPTO hindcast data, you will need to configure
% h5pyd for data access on HSDS.
% See the WPTO_hindcast example notebook for more details.
%
% Parameters
% ------------
% data_type : string
% data set type of interst
% Options: '3-hour', '1-hour'
%
% parameter : string or cell array of strings
% dataset parameter to be downloaded
% 3-hour dataset options: 'directionality_coefficient', 'energy_period', 'maximum_energy_direction'
% 'mean_absolute_period', 'mean_zero-crossing_period', 'omni-directional_wave_power', 'peak_period'
% 'significant_wave_height', 'spectral_width', 'water_depth'
% 1-hour dataset options: 'directionality_coefficient', 'energy_period', 'maximum_energy_direction'
% 'mean_absolute_period', 'mean_zero-crossing_period', 'omni-directional_wave_power', 'peak_period'
% 'significant_wave_height', 'spectral_width', 'water_depth', 'maximim_energy_direction',
% 'mean_wave_direction', 'frequency_bin_edges'
%
% lat_lon: cell array or cell array of cell arrays
% latitude longitude pairs at which to extract data
%
% years : cell array
% Vector of years to access. The years 1979-2010 available.
% Examples: {1996} or {2004,2006,2007}
%
% tree : str | cKDTree (optional)
% cKDTree or path to .pkl file containing pre-computed tree
% of lat, lon coordinates, default = py.None
%
% unscale : bool (optional)
% Boolean flag to automatically unscale variables on extraction
% Default = py.True
%
% str_decode : bool (optional)
% Boolean flag to decode the bytestring meta data into normal
% strings. Setting this to False will speed up the meta data read.
% Default = py.True
%
% hsds : bool (optional)
% Boolean flag to use h5pyd to handle .h5 'files' hosted on AWS
% behind HSDS. Setting to False will indicate to look for files on
% local machine, not AWS. Default = py.True
%
% Returns
% ---------
% data: Structure
%
%
% data.Data: named according to header row
%
% data.time: given in datetime
%
% data.units: the units for each data entry
%
% OR if a spectra data NDBC file
%
% data.spectra: spectra data
%
% data.time: given in datetime
%
% data.frequency: spectral frequency
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

arguments
data_type
parameter
lat_lon
years
options.tree = py.None;
options.unscale = py.True;
options.str_decode = py.True;
options.hsds = py.True;
end

py.importlib.import_module('mhkit');
py.importlib.import_module('numpy');

datap = py.mhkit.wave.io.hindcast.request_wpto_point_data(data_type,py.list(parameter),...
lat_lon,py.list(years),pyargs("tree",options.tree,"unscale",options.unscale,...
"str_decode",options.str_decode,"hsds",options.hsds));

datac=cell(datap);
datapd=datac{1};

datamat=datac{2};
matstr=datamat;
meta_ind = cell(py.list(py.numpy.nditer(matstr.index,pyargs("flags",{"refs_ok"}))));
meta_col = cell(py.list(py.numpy.nditer(matstr.columns,pyargs("flags",{"refs_ok"}))));
data_col = cell(py.list(py.numpy.nditer(datapd.columns,pyargs("flags",{"refs_ok"}))));
data_ind = cell(py.list(py.numpy.nditer(datapd.index,pyargs("flags",{"refs_ok"}))));
meta_val = cell(py.list(py.numpy.nditer(matstr.values,pyargs("flags",{"refs_ok"}))));
%disp(meta_ind{1})
%disp(class(meta_val{1}))
%disp(cell(py.list(py.numpy.nditer(datamat.columns.values,pyargs("flags",{"refs_ok"})))))

xx=cell(datapd.axes);
v=xx{2};

vv=cell(py.list(py.numpy.nditer(v.values,pyargs("flags",{"refs_ok"}))));

vals=double(py.array.array('d',py.numpy.nditer(datapd.values,pyargs("flags",{"refs_ok"}))));
sha=cell(datapd.values.shape);
x=int64(sha{1,1});
y=int64(sha{1,2});

vals=reshape(vals,[x,y]);
ti=cell(py.list(py.numpy.nditer(datapd.index,pyargs("flags",{"refs_ok"}))));
siti=size(ti);
si=size(vals);
temp = [];

for k = 1:max(size(data_col))
dat = char(py.str(data_col{k}));
num = dat(end);
num2 = num+1;
dat = dat(1:end-2);

numst = string(num);
dat = string(strrep(dat,'-','_'));
test='location'+string(numst);
num = str2num(num);

datast.(test).(dat)=vals(:,k);
%datast.(test).time = data_ind;

for j = 1:max(size(meta_ind))
if int64(meta_ind{j}) == num
m = num+1;
c = 1;
while m <= max(size(meta_val))
met = string(py.str(meta_col{c}));

if met == 'jurisdiction'
datast.(test).(met)=string(py.str(meta_val{m}));
else
datast.(test).(met)=double(py.array.array...
('d',py.numpy.nditer(meta_val{m},pyargs("flags",{"refs_ok"}))));
end
m = m+max(size(meta_ind));

c = c + 1;
end
end

end
%unit=string(matstr.(test));
%datast.units.(test)=unit;
end

si = size(data_ind);
for i=1:si(2)
datast.time(i)=datetime(string(py.str(data_ind{i})),'InputFormat','yyyy-MM-dd HH:mm:ssXXX',...
'TimeZone','UTC')';
end
end
%disp(datast)
% else
% datast.spectrum = vals';
% for i=1:si(2)
% temp = [temp, double(py.array.array('d',py.numpy.nditer(vv{i})))];
% end
% datast.frequency = temp';
% end
%
% for i=1:siti(2)
% datast.time(i)=datetime(string(py.str(ti{i})),'InputFormat','yyyy-MM-dd''T''HH:mm:ss.SSSSSSSSS')';
% end