Ebbflood

examples/river_example.html

@@ -92,10 +92,10 @@
 % they are most often exceeded.
 
 % Calculating the exceedance probability structure (F) from the Discharge data structure
-F = exceedance_probability(data); 
+F = exceedance_probability(data.Discharge) 
 
 % plotting the flow duration curve from Discharge and Exceedance Probability vectors
-plot_flow_duration_curve(data.Discharge,F.F); 
+plot_flow_duration_curve(data.Discharge,F); 
 %% Velocity Duration Curve
 % At each river energy converter location, we must provide a curve that relates 
 % the velocity at the turbine location to the river discharge levels. IEC 301 
@@ -129,10 +129,10 @@
 
 % Use polynomial fit from DV curve to calculate velocity('V') from 
 %  discharge at turbine location
-V = discharge_to_velocity(data, p.coef); % Input of discharge structure and polynomial coeficients 
+V = discharge_to_velocity(data, p.coef) % Input of discharge structure and polynomial coeficients 
 
 % Plot the velocity duration curve (VDC) 
-p2 = plot_velocity_duration_curve(V.V, F.F); % Input velocity vector and exceedance probability vector
+p2 = plot_velocity_duration_curve(V.V, F); % Input velocity vector and exceedance probability vector
 %% Power Duration Curve
 % The power duration curve is created in a nearly identical manner to the VDC. 
 % Here a velocity to power curve will be used to create a polynomial fit. The 
@@ -155,7 +155,7 @@
 P = velocity_to_power(V,p_VP.coef,min(VP_curve.V), max(VP_curve.V)); % Input is velocity structure, poly-coef, min and max of VP curve
 
 % Plot the power duration curve
-p4 = plot_power_duration_curve(P.P, F.F); % Input of power vector and exceedance probability vector
+p4 = plot_power_duration_curve(P.P, F); % Input of power vector and exceedance probability vector
 %% Annual Energy Produced 
 % Finally, annual energy produced is computed, as shown below.
 

mhkit/river/resource/exceedance_probability.m

@@ -5,81 +5,29 @@
 %    
 % Parameters
 % ----------
-%     Q : Discharge data [m3/s]
-%
-%         Pandas dataframe indexed by time [datetime or s]  
-%
-%           To make a pandas data frame from user supplied frequency and spectra
-%           use py.mhkit_python_utils.pandas_dataframe.timeseries_to_pandas(timeseries,time,x)
-%
-%         OR
-%
-%         structure of form:
-%
-%           Q.Discharge
-%
-%           Q.time
+%     Q : Array
+%       Discharge data [m3/s]
 %         
 % Returns   
 % -------
-%     F : Structure 
-%
-%
-%         F.F: Exceedance probability [unitless] 
-%
-%         F.time: time [epoch time (s)]
-%
+%     F : Array
+%       Exceedance probability [unitless] 
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+% First we make a sorted array of discharge values
+sort_q = sort(Q,'descend');
+% Next we make a ranking array
+rank = [1:numel(sort_q)];
+% Now each value in Q is assigned a rank where rank is the same as pandas
+% ranking algorithm with the 'max' option. For example in the data set 
+% [1, 2, 2, 3] the ranks would be [ 1, 2, 3, 4] and the rank of 2 would be
+% 3 because 3 is the maximum rank that a 2 has in the total set. 
+assigned_rank = arrayfun(@(x) rank(find(sort_q==x,1,'last')), Q);
 
-py.importlib.import_module('mhkit_python_utils');
-py.importlib.import_module('mhkit');
-
-if (isa(Q,'py.pandas.core.frame.DataFrame')~=1)
-    x=size(Q.Discharge);
-    li=py.list();
-    if x(2)>1 
-        for i = 1:x(2)
-            app=py.list(Q.Discharge(:,i));
-            li=py.mhkit_python_utils.pandas_dataframe.lis(li,app);
-
-        end
-    elseif x(2) ==1 
-        li=Q.Discharge;
-    end
-
-    if any(isdatetime(Q.time(1)))
-        si=size(Q.time);
-        for i=1:si(2)
-        Q.time(i)=posixtime(Q.time(i));
-        end
-    end
-    Q=py.mhkit_python_utils.pandas_dataframe.timeseries_to_pandas(li,Q.time,int32(x(2)));
-end
-
-EPpd=py.mhkit.river.resource.exceedance_probability(Q);
-
-xx=cell(EPpd.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(EPpd.values)));
-sha=cell(EPpd.values.shape);
-x=int64(sha{1,1});
-y=int64(sha{1,2});
+F = 100 * (assigned_rank / (numel(Q)+1));
 
-vals=reshape(vals,[x,y]);
 
-si=size(vals);
- for i=1:si(2)
-    test=string(py.str(vv{i}));
-    newname=split(test,",");
-
-    F.(newname(1))=vals(:,i);
-
- end
- F.time=double(py.array.array('d',py.numpy.nditer(EPpd.index)));
 
 
 

mhkit/tests/River_TestResource.m

@@ -11,13 +11,13 @@ function test_Froude_number(testCase)
 
         function test_exceedance_probability(testCase)
             % Create arbitrary discharge between 0 and 8(N=9)
-            Q = struct('Discharge',[0;1;2;3;4;5;6;7;8],'time',[0 1 2 3 4 5 6 7 8]);
+            Q = [0;1;2;3;4;5;6;7;8];
             % Rank order for non-repeating elements simply adds 1 to each element
             %if N=9, max F = 100((max(Q)+1)/10) =  90%
             %if N=9, min F = 100((min(Q)+1)/10) =  10%
             f = exceedance_probability(Q);
-            assertEqual(testCase,min(f.F), 10);
-            assertEqual(testCase,max(f.F), 90);
+            assertEqual(testCase,min(f), 10);
+            assertEqual(testCase,max(f), 90);
         end
 
         function test_polynomial_fit(testCase)
@@ -103,9 +103,9 @@ function test_plot_flow_duration_curve(testCase)
             end
 
             Q = struct('Discharge',[0;1;2;3;4;5;6;7;8],'time',[0 1 2 3 4 5 6 7 8]);
-            f = exceedance_probability(Q);
+            f = exceedance_probability(Q.Discharge);
 
-            plot_flow_duration_curve(Q.Discharge, f.F,"savepath",filename);
+            plot_flow_duration_curve(Q.Discharge, f,"savepath",filename);
 
             assertTrue(testCase,isfile(filename));
             delete(filename);
@@ -118,7 +118,7 @@ function test_plot_power_duration_curve(testCase)
             end
 
             Q = struct('Discharge',[0;1;2;3;4;5;6;7;8],'time',[0 1 2 3 4 5 6 7 8]);
-            f = exceedance_probability(Q);
+            f = exceedance_probability(Q.Discharge);
             poly = polynomial_fit(0:1:8, 10*(0:1:8),1);
             V = discharge_to_velocity(Q, poly.coef);
             VV = V.V;
@@ -130,7 +130,7 @@ function test_plot_power_duration_curve(testCase)
             % Power should be 10x greater and exclude the ends of V
             P = velocity_to_power(V, poly.coef, cut_in, cut_out);
 
-            plot_power_duration_curve(P.P, f.F,"savepath",filename);
+            plot_power_duration_curve(P.P, f,"savepath",filename);
 
             assertTrue(testCase,isfile(filename));
             delete(filename);
@@ -145,9 +145,9 @@ function test_plot_velocity_duration_curve(testCase)
             Q = struct('Discharge',[0;1;2;3;4;5;6;7;8],'time',[0 1 2 3 4 5 6 7 8]);
             poly = polynomial_fit(0:1:8, 10*(0:1:8),1);
             V = discharge_to_velocity(Q, poly.coef);
-            f = exceedance_probability(Q);
+            f = exceedance_probability(Q.Discharge);
 
-            plot_velocity_duration_curve(V.V, f.F,"savepath",filename);
+            plot_velocity_duration_curve(V.V, f,"savepath",filename);
 
             assertTrue(testCase,isfile(filename));
             delete(filename);
@@ -190,7 +190,6 @@ function test_plot_velocity_vs_power(testCase)
             end
 
             Q = struct('Discharge',[0;1;2;3;4;5;6;7;8],'time',[0 1 2 3 4 5 6 7 8]);
-            f = exceedance_probability(Q);
             poly = polynomial_fit(0:1:8, 10*(0:1:8),1);
             V = discharge_to_velocity(Q, poly.coef);
             VV = V.V;

mhkit/tests/Tidal_TestResource.m

@@ -3,10 +3,10 @@
     methods (Test) 
 
         function test_exceedance_probability(testCase)
-            Q = struct('Discharge',[0;1;2;3;4;5;6;7;8],'time',[1 2 3 4 5 6 7 8 9]);
+            Q = [0;1;2;3;4;5;6;7;8];
             f = exceedance_probability(Q);
-            assertEqual(testCase,min(f.F), 10);
-            assertEqual(testCase,max(f.F), 90);
+            assertEqual(testCase,min(f), 10);
+            assertEqual(testCase,max(f), 90);
         end
 
         function test_principal_flow_directions(testCase)
@@ -88,8 +88,9 @@ function test_plot_phase_probability(testCase)
                 delete(filename);
             end
 
-            file_name = '../../examples/data/tidal/s08010.json';
-            data = read_noaa_json(file_name);
+            relative_file_name = '../../examples/data/tidal/s08010.json';
+            full_file_name = fullfile(fileparts(mfilename('fullpath')), relative_file_name);
+            data = read_noaa_json(full_file_name);
             data.s = data.s/100;
             width_direction = 1;
             [direction1, direction2] = ...
@@ -109,10 +110,10 @@ function test_plot_phase_exceedance(testCase)
                 delete(filename);
             end
 
-            file_name = '../../examples/data/tidal/s08010.json';
-            data = read_noaa_json(file_name);
+            relative_file_name = '../../examples/data/tidal/s08010.json';
+            full_file_name = fullfile(fileparts(mfilename('fullpath')), relative_file_name);
+            data = read_noaa_json(full_file_name);
             data.s = data.s/100;
-            data.Discharge = data.s;
             width_direction = 1;
             [direction1, direction2] = ...
                 principal_flow_directions(data.d,width_direction);

mhkit/tests/Wave_TestIO.m

@@ -110,7 +110,7 @@ function test_swan_read_block(testCase)
 
         % WPTO multiple locations
         function test_WPTO_point_multiloc(testCase)
-
+            if ispc
             api_key = '3K3JQbjZmWctY0xmIfSYvYgtIcM3CN0cb1Y2w9bf';
             hindcast_data = request_wpto('1-hour',...
                 ["energy_period"],[44.624076,-124.280097;43.489171,-125.152137],...
@@ -140,11 +140,13 @@ function test_WPTO_point_multiloc(testCase)
             assertEqual(testCase,expected_meta.timezone(2),hindcast_data(2).metadata.timezone);
             assertEqual(testCase,expected_meta.jurisdiction{2},hindcast_data(2).metadata.jurisdiction);
             assertEqual(testCase,expected_meta.distance_to_shore(2),hindcast_data(2).metadata.distance_to_shore,'RelTol',0.000001);
-
+            else
+
+            end
         end
 
         function test_WPTO_point_multiparm(testCase)
-
+            if ismac
             api_key = '3K3JQbjZmWctY0xmIfSYvYgtIcM3CN0cb1Y2w9bf';
             hindcast_data = request_wpto('3-hour',...
                 ["mean_absolute_period","significant_wave_height"],[44.624076,-124.280097],...
@@ -166,6 +168,10 @@ function test_WPTO_point_multiparm(testCase)
             assertEqual(testCase,expected_meta.timezone(1),hindcast_data.metadata.timezone);
             assertEqual(testCase,expected_meta.jurisdiction{1},hindcast_data.metadata.jurisdiction);
             assertEqual(testCase,expected_meta.distance_to_shore(1),hindcast_data.metadata.distance_to_shore,'RelTol',0.000001);
+
+            else
+
+            end
         end
     end
 end

mhkit/tidal/graphics/plot_tidal_phase_exceedance.m

@@ -80,34 +80,27 @@
 	isEbb = ~((data.d < upper_split) & (data.d >= lower_split));
 end
 
-s_ebb = struct('time', {data.time(isEbb)},...
-    'Discharge', {data.Discharge(isEbb)});
-s_flood = struct('time', {data.time(~isEbb)},...
-    'Discharge', {data.Discharge(~isEbb)});
+s_ebb = data.s(isEbb);
+s_flood = data.s(~isEbb);
 
-F = exceedance_probability(data).F;
-F_ebb = exceedance_probability(s_ebb).F;
-F_flood = exceedance_probability(s_flood).F;
+F_ebb = exceedance_probability(s_ebb);
+F_flood = exceedance_probability(s_flood);
 
 decimals = round(options.bin_size/0.1);
 s_new = [round(min(data.s),decimals):options.bin_size:...
     round(max(data.s),decimals)+options.bin_size];
 s_new = s_new(1:end-1); % this accounts for the difference between matlab 
                         % and numpy.arange which uses a half open interval
 
-py.importlib.import_module('scipy');
-kwa = pyargs('bounds_error',false);
+[~,i_ebb] = unique(s_ebb);
+[~,i_flood] = unique(s_flood);
 
-f_ebb = py.scipy.interpolate.interp1d(s_ebb.Discharge,...
+f_ebb = py.scipy.interpolate.interpolate.interp1d(s_ebb.Discharge,...
     F_ebb, kwa);
-f_flood = py.scipy.interpolate.interp1d(s_flood.Discharge,...
+f_flood = py.scipy.interpolate.interpolate.interp1d(s_flood.Discharge,...
     F_flood, kwa);
 
-F_ebb = double(f_ebb(s_new));
-F_flood = double(f_flood(s_new));
-
-py.importlib.import_module('numpy');
-F_max_total = py.numpy.nanmax(F_ebb) + py.numpy.nanmax(F_flood);
+F_max_total = max(F_ebb(~isnan(F_ebb))) + max(F_flood(~isnan(F_flood)));
 
 s_plot = repmat(s_new,2,1);
 figure = area(s_plot.', [F_ebb/F_max_total*100; F_flood/F_max_total*100].');
@@ -125,4 +118,4 @@
     exportgraphics(gca, options.savepath);
 end 
 
-hold off
+hold off