Home > mrg > MRG_MIKE > mrg_dfsu_exceed.m

mrg_dfsu_exceed

PURPOSE ^

Performs an exceedance analysis on a DFSU files.

SYNOPSIS ^

function [data] = mrg_dfsu_exceed(exceed, chunk_size)

DESCRIPTION ^

 Performs an exceedance analysis on a DFSU files.

 INPUT
   exceed          A vector representing the quantiles to extract.  
                   Defualts to [.10 .20 .30 .50]
   chunk_size      An integer specifying the size of the 'chunks' to break 
                   the file up into.  Larger chunks may speed up runtime, 
                   but very large chunks may lead to out of memory errors.  
                   Defaults to 1000.

 OUTPUT
   A MATLAB structure with n+3 items, where n is the number of 
   quantiles requested, plus X, Y and Z data from the DFSU file.  
   A CSV file, derived from the MATLAB structure.
   A DFSU file with the same spatial extent as the analysed file, but 
   with a single timestep will be written as 'inputfilename_exceed.dfsu'.

 REQUIREMENTS
   The DHI/MIKE Matlab toolbox 2011 (developed with v. 20110304)
   mrg_struct_to_csv.m function (assuming you want csv output, else it will be skipped)

 NOTES
     Extracts quantiles from each cell in a DFSU file (i.e. an analysis
   through the time domain).

 LICENCE
   Created by Daniel Pritchard (www.pritchard.co)
   Distributed under a creative commons CC BY-SA licence. See here:
   http://creativecommons.org/licenses/by-sa/3.0/

 DEVELOPMENT
   v 1.0   6/7/2012
           DP. Initial attempt and distribution. 
   v 1.1   12/7/2012
           DP. Added better filenaming and removed error-causing 'factory' call. 
   v 1.2   19/7/2012
           DP. Now deals with delete values by using single() throughout (use of
           double() causes issues with rounding).
           DP. Much better estimation of time remaining.
   v 1.3   14/02/2013
           DP. Documentation

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [data] = mrg_dfsu_exceed(exceed, chunk_size)
0002 % Performs an exceedance analysis on a DFSU files.
0003 %
0004 % INPUT
0005 %   exceed          A vector representing the quantiles to extract.
0006 %                   Defualts to [.10 .20 .30 .50]
0007 %   chunk_size      An integer specifying the size of the 'chunks' to break
0008 %                   the file up into.  Larger chunks may speed up runtime,
0009 %                   but very large chunks may lead to out of memory errors.
0010 %                   Defaults to 1000.
0011 %
0012 % OUTPUT
0013 %   A MATLAB structure with n+3 items, where n is the number of
0014 %   quantiles requested, plus X, Y and Z data from the DFSU file.
0015 %   A CSV file, derived from the MATLAB structure.
0016 %   A DFSU file with the same spatial extent as the analysed file, but
0017 %   with a single timestep will be written as 'inputfilename_exceed.dfsu'.
0018 %
0019 % REQUIREMENTS
0020 %   The DHI/MIKE Matlab toolbox 2011 (developed with v. 20110304)
0021 %   mrg_struct_to_csv.m function (assuming you want csv output, else it will be skipped)
0022 %
0023 % NOTES
0024 %     Extracts quantiles from each cell in a DFSU file (i.e. an analysis
0025 %   through the time domain).
0026 %
0027 % LICENCE
0028 %   Created by Daniel Pritchard (www.pritchard.co)
0029 %   Distributed under a creative commons CC BY-SA licence. See here:
0030 %   http://creativecommons.org/licenses/by-sa/3.0/
0031 %
0032 % DEVELOPMENT
0033 %   v 1.0   6/7/2012
0034 %           DP. Initial attempt and distribution.
0035 %   v 1.1   12/7/2012
0036 %           DP. Added better filenaming and removed error-causing 'factory' call.
0037 %   v 1.2   19/7/2012
0038 %           DP. Now deals with delete values by using single() throughout (use of
0039 %           double() causes issues with rounding).
0040 %           DP. Much better estimation of time remaining.
0041 %   v 1.3   14/02/2013
0042 %           DP. Documentation
0043 
0044 %% Go!
0045 % Exceedance values (a vector)
0046 if ~exist('exceed', 'var')
0047     exceed = [.10 .20 .30 .50];
0048 end
0049 
0050 % Chunk size...
0051 if ~exist('chunk_size', 'var')
0052     chunk_size = 1000;
0053 end
0054 
0055 % Load libraries
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 % Get file
0063 [filename,filepath] = uigetfile('*.dfsu','Select the .dfsu file to analyse');
0064 infile = [filepath,filename];
0065 dfsu_file = DfsFileFactory.DfsuFileOpen(infile);
0066 
0067 % Read some item information
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 % Ask the user to select the item to process
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; % A number delimiting the item position in the DFSU file.  The same as mat_item (i.e. no zero indexing).  i.e. this is completly redundant!
0082 end
0083 
0084 nsteps = dfsu_file.NumberOfTimeSteps;
0085 nelements = dfsu_file.NumberOfElements;
0086 
0087 % Setup the output object
0088 data = struct();
0089 
0090 % Add some nicities, like X, Y, Z
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 % Add (empty) objects for the exceedance values
0101 for a = 1:length(exceed)
0102     data.(['quant_',regexprep(num2str(exceed(a)),'\.','')]) = NaN(1,nelements, 'single');
0103 end
0104 
0105 % Pre-allocate some memory, and setup our sequence...
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 % Single delete value
0113 sDelete = single(1e-35);
0114 
0115 % Timings are for debugging and display only...
0116 timings = NaN(1,seq_len);
0117 
0118 % loc = Keeping track of the loctaion in the data file where we are at
0119 loc = 0;
0120 % OK this is ugly, but I think I implimented this after 4 beers at about
0121 % 10 pm sitting at a coffee table in a rented apartment in Denmark.  If you
0122 % wanted professional code you wouldn't rely on drunk, tired, algal
0123 % physiologists!
0124 
0125 for j=1:seq_len
0126     tic
0127     if j == seq_len
0128         %This is the last chunk.  Act accordingly.
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         % OK now we have a chunk.  We need to process it...
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         % Convert delete values to NANs
0140         temp_el_chunk(temp_el_chunk==sDelete) = NaN;
0141         
0142         % Do the heavy lifting
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         % OK now we have a chunk.  We need to process it...
0155         nel_in_chunk = size(temp_el_chunk,1);
0156         %disp(['Processing element ', num2str(j), ' of ', num2str(nelements), ' in chunk ', num2str(j), ' of ', num2str(seq_len), '.'])
0157         disp(['Processing elements ', num2str(loc+1), ' to ' num2str(size(temp_el_chunk,1)+loc), ' of ', num2str(nelements)])
0158         
0159         % Convert delete values to NANs
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 %% Here we write out a CSV file
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 %% Here we  stick data into a DFSU file
0186 % First DRAFT of the DFSU - based output...  Seems to work OK...
0187 disp('Writing exceedance values to a DFSU file.')
0188 
0189 % Apparently we need to load these again...  Hmmm... TODO: Figure this out
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 % Being lazy here - no need to write out X, Y and Z - TODO: Remove this...
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 % Follow up from above - Re: Layzness
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

Generated on Thu 29-May-2014 21:29:53 by m2html © 2005