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

%{
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)
    fprintf('NetCDF42C v1.1.14\n');

    % make a model framework to fill that is in the scope of this file
    global model_copy;
    model_copy = model();

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

        % Open the given netCDF4 file
        global NCData;
        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");
            make_results_subclasses();
        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");
            check_inversion_class();
        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
            walk_nested_groups(group_id);
        end
        
        % Close the netCDF file
        netcdf.close(NCData);
        disp('Model Successfully Copied')
    else
        fprintf('File %s does not exist.\n', filename);
    end
end


function make_results_subclasses()
    global model_copy;
    global NCData;
    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
    disp('Successfully recreated results structs:')
    for fieldname = string(fieldnames(model_copy.results))
        disp(fieldname)
    end
end


function check_inversion_class()
    % get the name of the inversion class: either inversion or m1qn3inversion or taoinversion
    global model_copy;
    global NCData;
    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();
        disp('Successfully created inversion class instance: m1qn3inversion')
    elseif strcmp(inversion_class, 'taoinversion')
        model_copy.inversion = taoinversion();
        disp('Successfully created inversion class instance: taoinversion')
    else
        disp('No inversion class was found')
    end
end



function walk_nested_groups(group_location_in_file)
    global model_copy;
    global NCData;    
    % 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;
            copy_nested_struct(group_location_in_file)
        elseif strcmp(varname, 'solution')
            % band-aid pass..
        else
            copy_variable_data_to_new_model(group_location_in_file, varname, xtype);
        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))
                walk_nested_groups(group);
            end
        end
    end
    %catch % no nested groups here
    %end
end



function copy_nested_struct(group_location_in_file)
    global model_copy;
    global NCData;
    %{
        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);
        fprintf("Successfully saved layer %s to multidimension struct array\n", num2str(current_layer))
    end
    fprintf('Successfully recreated multidimensional structure array %s in md.%s\n', name_of_struct, address_in_model)
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 copy_variable_data_to_new_model(group_location_in_file, varname, xtype)
    global model_copy;
    global NCData;
    %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('saved int64 as int16')
            else
                arg_to_eval = ['model_copy', address_to_attr, '.', varname, ' = data;'];
                eval(arg_to_eval);
            end
            
            full_addy = netcdf.inqGrpNameFull(group_location_in_file);
            %disp(xtype)
            %class(data)
            fprintf('Successfully saved %s to %s\n', varname, full_addy);

        catch Me %e 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
end