0001 function [data] = mrg_dfsu_exceed(exceed, chunk_size)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046 if ~exist('exceed', 'var')
0047 exceed = [.10 .20 .30 .50];
0048 end
0049
0050
0051 if ~exist('chunk_size', 'var')
0052 chunk_size = 1000;
0053 end
0054
0055
0056 NET.addAssembly('DHI.Generic.MikeZero.DFS');
0057 NET.addAssembly('DHI.Generic.MikeZero.EUM');
0058 import DHI.Generic.MikeZero.DFS.*;
0059 import DHI.Generic.MikeZero.DFS.dfsu.*;
0060 import DHI.Generic.MikeZero.*
0061
0062
0063 [filename,filepath] = uigetfile('*.dfsu','Select the .dfsu file to analyse');
0064 infile = [filepath,filename];
0065 dfsu_file = DfsFileFactory.DfsuFileOpen(infile);
0066
0067
0068 items = {};
0069 for i = 0:dfsu_file.ItemInfo.Count-1
0070 item = dfsu_file.ItemInfo.Item(i);
0071 items{i+1,1} = char(item.Name);
0072 items{i+1,2} = char(item.Quantity.Unit);
0073 items{i+1,3} = char(item.Quantity.UnitAbbreviation);
0074 end
0075
0076
0077 choice = menu(sprintf('Which of these objects do you want to process?'),items{:,1}, 'Ummm... None of the above!');
0078 if choice == 0 || choice == length(items(:,1))+1
0079 error('Do or do not. There is no try.');
0080 else
0081 dfsu_item = choice;
0082 end
0083
0084 nsteps = dfsu_file.NumberOfTimeSteps;
0085 nelements = dfsu_file.NumberOfElements;
0086
0087
0088 data = struct();
0089
0090
0091 node_x = single(dfsu_file.X);
0092 node_y = single(dfsu_file.Y);
0093 node_z = single(dfsu_file.Z);
0094 ele_table = mzNetFromElmtArray(dfsu_file.ElementTable);
0095 [ele_X,ele_Y,ele_Z] = mzCalcElmtCenterCoords(ele_table,node_x,node_y,node_z);
0096 data.X = ele_X.';
0097 data.Y = ele_Y.';
0098 data.Z = ele_Z.';
0099
0100
0101 for a = 1:length(exceed)
0102 data.(['quant_',regexprep(num2str(exceed(a)),'\.','')]) = NaN(1,nelements, 'single');
0103 end
0104
0105
0106 temp_single_ts = NaN(1,nelements, 'single');
0107 temp_el_chunk = NaN(chunk_size,nsteps, 'single');
0108 seq = int32([1:chunk_size:nelements]);
0109 seq_len = size(seq,2);
0110 seq(end+1) = int32(nelements);
0111
0112
0113 sDelete = single(1e-35);
0114
0115
0116 timings = NaN(1,seq_len);
0117
0118
0119 loc = 0;
0120
0121
0122
0123
0124
0125 for j=1:seq_len
0126 tic
0127 if j == seq_len
0128
0129 last_chunk_size = seq(j+1)-seq(j)+1;
0130 temp_el_chunk = ones(last_chunk_size,nsteps);
0131 for i=0:nsteps-1
0132 temp_single_ts(:) = single(dfsu_file.ReadItemTimeStep(dfsu_item,i).Data)';
0133 temp_el_chunk(:,i+1) = temp_single_ts(seq(j):seq(j+1));
0134 end
0135
0136 nel_in_chunk = size(temp_el_chunk,1);
0137 disp(['Processing elements ', num2str(loc+1), ' to ' num2str(size(temp_el_chunk,1)+loc), ' of ', num2str(nelements)])
0138
0139
0140 temp_el_chunk(temp_el_chunk==sDelete) = NaN;
0141
0142
0143 quant = quantile(temp_el_chunk,exceed,2);
0144
0145 for a = 1:size(quant,2)
0146 data.(['quant_',regexprep(num2str(exceed(a)),'\.','')])(loc+1:size(temp_el_chunk,1)+loc) = quant(:,a);
0147 end
0148 loc = size(temp_el_chunk,1)+loc;
0149 else
0150 for i=0:nsteps-1
0151 temp_single_ts(:) = single(dfsu_file.ReadItemTimeStep(dfsu_item,i).Data)';
0152 temp_el_chunk(:,i+1) = temp_single_ts(seq(j):seq(j+1)-1);
0153 end
0154
0155 nel_in_chunk = size(temp_el_chunk,1);
0156
0157 disp(['Processing elements ', num2str(loc+1), ' to ' num2str(size(temp_el_chunk,1)+loc), ' of ', num2str(nelements)])
0158
0159
0160 temp_el_chunk(temp_el_chunk==sDelete) = NaN;
0161
0162 quant = quantile(temp_el_chunk,exceed,2);
0163 for a = 1:size(quant,2)
0164 data.(['quant_',regexprep(num2str(exceed(a)),'\.','')])(loc+1:size(temp_el_chunk,1)+loc) = quant(:,a);
0165 end
0166 loc = size(temp_el_chunk,1)+loc;
0167 end
0168 timings(j) = toc;
0169 mean_time_sec = nanmean(timings);
0170 remain = (seq_len-j)*mean_time_sec/60;
0171 disp(['Approximatly ', num2str(remain), ' minutes remain.'])
0172 end
0173
0174
0175 if exist('mrg_struct_to_csv.m', 'file') == 2
0176 disp('Writing exceedance values to a CSV file.')
0177 out_file = [filename(1:end-5),'_',regexprep(items{choice,1},' ',''),'_exceed.csv'];
0178 out_file_full = [filepath,out_file];
0179 mrg_struct_to_csv(data, out_file_full);
0180 disp(['Output CSV file (', out_file, ') written successfully.'])
0181 else
0182 disp('mrg_struct_to_csv.m not avilable in your search path. No CSV file will be generated.')
0183 end
0184
0185
0186
0187 disp('Writing exceedance values to a DFSU file.')
0188
0189
0190 NET.addAssembly('DHI.Generic.MikeZero.DFS');
0191 NET.addAssembly('DHI.Generic.MikeZero.EUM');
0192 import DHI.Generic.MikeZero.DFS.*;
0193 import DHI.Generic.MikeZero.DFS.dfsu.*;
0194 import DHI.Generic.MikeZero.*
0195 import DHI.Generic.MikeZero.EUM.*
0196
0197 builder = DfsuBuilder.Create(DfsuFileType.Dfsu2D);
0198 builder.SetTimeInfo(dfsu_file.StartDateTime, 1);
0199 builder.SetNodes(dfsu_file.X,dfsu_file.Y,dfsu_file.Z,dfsu_file.Code);
0200 builder.SetElements(dfsu_file.ElementTable);
0201 builder.SetProjection(dfsu_file.Projection);
0202
0203 data_names = fieldnames(data);
0204
0205
0206 for n = 1:length(data_names)
0207 builder.AddDynamicItem(data_names{n},eumQuantity(eumItem.eumIItemUndefined,eumUnit.eumUUnitUndefined));
0208 end
0209
0210 out_file = [filename(1:end-5),'_',regexprep(items{choice,1},' ',''),'_exceed', filename(end-4:end)];
0211 out_file_full = [filepath,out_file];
0212
0213 dfs_out = builder.CreateFile(out_file_full);
0214
0215 if ~(size(data.(data_names{1}), 2) == dfs_out.NumberOfElements)
0216 error(['Somehow there is a mismatch in the number of elements in the ',...
0217 'DFSU output file and the length of the data. This is an error'])
0218 end
0219
0220
0221 for n = 1:length(data_names)
0222 data_out = single(data.(data_names{n}));
0223 data_out(isnan(data_out)) = sDelete;
0224 dfs_out.WriteItemTimeStepNext(n-1, NET.convertArray(data_out));
0225 end
0226
0227 dfs_out.Close()
0228 disp(['Output DFSU file (', out_file, ') written successfully.'])
0229
0230 end
0231
0232
0233
0234