source: issm/trunk/src/m/contrib/musselman/read_netCDF.m@ 27889

%{
Given a NetCDF4 file, this set of functions will perform the following:
    1. Enter each group of the file.
    2. For each variable in each group, update an empty model with the variable's data
    3. Enter nested groups and repeat


If the model you saved has subclass instances that are not in the standard model() class
you can:
    1. Copy lines 30-35, set the "results" string to the name of the subclass instance,
    2. Copy and modify the make_results_subclasses() function to create the new subclass 
        instances you need. 
From there, the rest of this script will automatically create the new subclass 
instance in the model you're writing to and store the data from the netcdf file there.
%}


function model_copy = read_netCDF(filename, varargin)
    if nargin > 1
        verbose = true;
    else
        verbose = false;
    end
    
    if verbose
        fprintf('NetCDF42C v1.1.14\n');
    end
    % make a model framework to fill that is in the scope of this file
    model_copy = model();

    % Check if path exists
    if exist(filename, 'file')
        if verbose
            fprintf('Opening %s for reading\n', filename);
        end

        % Open the given netCDF4 file
        NCData = netcdf.open(filename, 'NOWRITE');
        % Remove masks from netCDF data for easy conversion: NOT WORKING
        %netcdf.setMask(NCData, 'NC_NOFILL');

        % see if results is in there, if it is we have to instantiate some classes
        try
            results_group_id = netcdf.inqNcid(NCData, "results");
            model_copy = make_results_subclasses(model_copy, NCData, verbose);
        catch
        end % 'results' group doesn't exist 

        % see if inversion is in there, if it is we may have to instantiate some classes
        try
            inversion_group_id = netcdf.inqNcid(NCData, "inversion");
            model_copy = check_inversion_class(model_copy, NCData, verbose);
        catch
        end % 'inversion' group doesn't exist 
        
        % loop over first layer of groups in netcdf file
        for group = netcdf.inqGrps(NCData)
            group_id = netcdf.inqNcid(NCData, netcdf.inqGrpName(group));
            %disp(netcdf.inqGrpNameFull(group_id))
            % hand off first level to recursive search
            model_copy = walk_nested_groups(group_id, model_copy, NCData, verbose);
        end
        
        % Close the netCDF file
        netcdf.close(NCData);
        if verbose
            disp('Model Successfully Copied')
        end
    else
        fprintf('File %s does not exist.\n', filename);
    end
end


function model_copy = make_results_subclasses(model_copy, NCData, verbose)
    resultsGroup = netcdf.inqNcid(NCData, "results");
    variables = netcdf.inqVarIDs(resultsGroup);
    for name = variables
        class_instance = netcdf.inqVar(resultsGroup, name);
        class_instance_names_raw = netcdf.getVar(resultsGroup, name, 'char').';
        class_instance_names = cellstr(class_instance_names_raw);
        for index = 1:numel(class_instance_names)
            class_instance_name = class_instance_names{index};
            model_copy.results = setfield(model_copy.results, class_instance_name, struct());
        end
        %model_copy.results = setfield(model_copy.results, class_instance, class_instance_name);
    end
    model_copy = model_copy;
    if verbose
        disp('Successfully recreated results structs:')
        for fieldname = string(fieldnames(model_copy.results))
            disp(fieldname)
        end
    end
end


function model_copy = check_inversion_class(model_copy, NCData, verbose)
    % get the name of the inversion class: either inversion or m1qn3inversion or taoinversion
    inversionGroup = netcdf.inqNcid(NCData, "inversion");
    varid = netcdf.inqVarID(inversionGroup, 'inversion_class_name');
    inversion_class = convertCharsToStrings(netcdf.getVar(inversionGroup, varid,'char'));
    if strcmp(inversion_class, 'm1qn3inversion')
        model_copy.inversion = m1qn3inversion();
        if verbose
            disp('Successfully created inversion class instance: m1qn3inversion')
        end
    elseif strcmp(inversion_class, 'taoinversion')
        model_copy.inversion = taoinversion();
        if verbose
            disp('Successfully created inversion class instance: taoinversion')
        end
    else
        if verbose
            disp('No inversion class was found')
        end
    end
    model_copy = model_copy;
end


function model_copy = walk_nested_groups(group_location_in_file, model_copy, NCData, verbose)  
    % we search the current group level for variables by getting this struct
    variables = netcdf.inqVarIDs(group_location_in_file); 

    % from the variables struct get the info related to the variables
    for variable = variables
        [varname, xtype, dimids, numatts] = netcdf.inqVar(group_location_in_file, variable);
        
        % keep an eye out for nested structs:
        if strcmp(varname, 'this_is_a_nested')
            is_nested = true;
            model_copy = copy_nested_struct(group_location_in_file, model_copy, NCData, verbose);
        elseif strcmp(varname, 'solution')
            % band-aid pass..
        else
            model_copy = copy_variable_data_to_new_model(group_location_in_file, varname, xtype, model_copy, NCData, verbose);
        end
    end

    % try to find groups in current level, if it doesn't work it's because there is nothing there
    %try
    % if it's a nested struct the function copy_nested_struct has already been called
    if logical(exist('is_nested', 'var'))
        % do nothing
    else
        % search for nested groups in the current level to feed back to this function
        groups = netcdf.inqGrps(group_location_in_file);
        if not(isempty(groups))
            for group = groups
                group_id = netcdf.inqNcid(group_location_in_file, netcdf.inqGrpName(group));
                %disp(netcdf.inqGrpNameFull(group_id))
                model_copy = walk_nested_groups(group, model_copy, NCData, verbose);
            end
        end
    end
    %catch % no nested groups here
    %end
end



function model_copy = copy_nested_struct(group_location_in_file, model_copy, NCData, verbose)
    %{
        A common multidimensional struct array is the 1xn md.results.TransientSolution struct. 
        The process to recreate is as follows:
            1. Get the name of the struct from group name
            2. Get the fieldnames from the subgroups 
            3. Recreate the struct with fieldnames 
            4. Populate the fields with their respective values
    %}

    % step 1
    name_of_struct = netcdf.inqGrpName(group_location_in_file);

    % step 2
    subgroups = netcdf.inqGrps(group_location_in_file); % numerical cell array with ID's of subgroups
    % get single subgroup's data
    single_subgroup_ID = subgroups(1);
    subgroup_varids = netcdf.inqVarIDs(single_subgroup_ID);
    fieldnames = {};
    for variable = subgroup_varids
        [varname, xtype, dimids, numatts] = netcdf.inqVar(single_subgroup_ID, variable);
        fieldnames{end+1} = varname;
    end

    % step 3
    address_in_model_raw = split(netcdf.inqGrpNameFull(group_location_in_file), '/');
    address_in_model = address_in_model_raw{2};
    
    % we cannot assign a variable to represent this object as MATLAB treats all variables as copies
    % and not pointers to the same memory address
    % this means that if address_in_model has more than 1 layer, we need to modify the code. For now, 
    % we just hope this will do. An example of a no-solution would be model().abc.def.ghi.field
    
    model_copy.(address_in_model).(name_of_struct) = struct();
    % for every fieldname in the subgroup, create an empty field
    for fieldname = string(fieldnames)
        model_copy.(address_in_model).(name_of_struct).(fieldname) = {};
    end

    % use repmat to make the struct array multidimensional along the fields axis
    number_of_dimensions = numel(subgroups);
    model_copy.(address_in_model).(name_of_struct) = repmat(model_copy.(address_in_model).(name_of_struct), 1, number_of_dimensions);
    
    % step 4
    % for every layer of the multidimensional struct array, populate the fields
    for current_layer = 1:number_of_dimensions
        % choose subgroup
        current_layer_subgroup_ID = subgroups(current_layer);
        % get all vars
        current_layer_subgroup_varids = netcdf.inqVarIDs(current_layer_subgroup_ID);
        % get individual vars and set fields at layer current_layer
        for varid = current_layer_subgroup_varids
            [varname, xtype, dimids, numatts] = netcdf.inqVar(current_layer_subgroup_ID, varid);
            data = netcdf.getVar(current_layer_subgroup_ID, varid);

            % netcdf uses Row Major Order but MATLAB uses Column Major Order so we need to transpose all arrays w/ more than 1 dim
            if all(size(data)~=1) || xtype == 2
                data = data.';
            end
            
            % set the field
            model_copy.(address_in_model).(name_of_struct)(current_layer).(varname) = data;
            %address_to_struct_in_model = setfield(address_to_struct_in_model(current_layer), varname, data)
        end
        model_copy.(address_in_model).(name_of_struct)(current_layer);
        if verbose
            fprintf("Successfully loaded layer %s to multidimension struct array\n", num2str(current_layer))
        end
    end
    model_copy = model_copy;
    if verbose
        fprintf('Successfully recreated multidimensional structure array %s in md.%s\n', name_of_struct, address_in_model)
    end
end




%{
Since there are two types of objects that MATLAB uses (classes and structs), we have to check 
which object we're working with before we can set any fields/attributes of it. After this is completed,
we can write the data to that location in the model.
%}

function model_copy = copy_variable_data_to_new_model(group_location_in_file, varname, xtype, model_copy, NCData, verbose)
    %disp(varname)
    % this is an inversion band-aid
    if strcmp(varname, 'inversion_class_name') || strcmp(varname, 'name_of_struct') || strcmp(varname, 'solution')
        % we don't need this
    else
        % putting try/catch here so that any errors generated while copying data are logged and not lost by the try/catch in walk_nested_groups function
        try
            %disp(netcdf.inqGrpNameFull(group_location_in_file))
            %disp(class(netcdf.inqGrpNameFull(group_location_in_file)))
            address_to_attr = strrep(netcdf.inqGrpNameFull(group_location_in_file), '/', '.');
            varid = netcdf.inqVarID(group_location_in_file, varname);
            data = netcdf.getVar(group_location_in_file, varid);
            
    
            % if we have an empty string
            if xtype == 2 && isempty(all(data))
                data = cell(char());
            % if we have an empty cell-char array
            elseif numel(data) == 1 && xtype == 3 && data == -32767
                data = cell(char());
            elseif isempty(all(data))
                data = []
            end
            % band-aid for some cell-char-arrays:
            if xtype == 2 && strcmp(data, 'default')
                data = {'default'};
            end
            
            % netcdf uses Row Major Order but MATLAB uses Column Major Order so we need to transpose all arrays w/ more than 1 dim
            if all(size(data)~=1) || xtype == 2
                data = data.';
            end
    
            % if we have a list of strings
            if xtype == 2
                try
                    if strcmp(netcdf.getAtt(group_location_in_file, varid, "type_is"), 'cell_array_of_strings')
                        data = cellstr(data);
                    end
                catch
                    % no attr found so we pass
                end
            end
            
            % the issm c compiler does not work with int64 datatypes, so we need to convert those to int16
            % reference this (very hard to find) link for netcdf4 datatypes: https://docs.unidata.ucar.edu/netcdf-c/current/netcdf_8h_source.html
            %xtype
            if xtype == 10
                arg_to_eval = ['model_copy', address_to_attr, '.', varname, ' = ' , 'double(data);'];
                eval(arg_to_eval);
                %disp('Loaded int64 as int16')
            else
                arg_to_eval = ['model_copy', address_to_attr, '.', varname, ' = data;'];
                eval(arg_to_eval);
            end
            
            if verbose
                full_addy = netcdf.inqGrpNameFull(group_location_in_file);
                %disp(xtype)
                %class(data)
                fprintf('Successfully loaded %s to %s\n', varname, full_addy);
            end

        catch ME %ME is an MException struct
            % Some error occurred if you get here.
            fprintf(1,'There was an error with %s! \n', varname)
            errorMessage = sprintf('Error in function %s() at line %d.\n\nError Message:\n%s', ME.stack.name, ME.stack.line, ME.message);
            fprintf(1, '%s\n', errorMessage);
            uiwait(warndlg(errorMessage));
            %line = ME.stack.line
            %fprintf(1,'There was an error with %s! \n', varname)
            %fprintf('The message was:\n%s\n',ME.message);
            %fprintf(1,'The identifier was:\n%s\n',ME.identifier);
            
            % more error handling...
        end
    end
    model_copy = model_copy;
end