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

%{
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.13\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_name = convertCharsToStrings(netcdf.getVar(resultsGroup, name, 'char'));
        model_copy.results = setfield(model_copy.results, class_instance, class_instance_name);
    end
    disp('Successfully recreated results struct')
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;
    % try to find vars in current level, if it doesn't work it's because there is nothing there
    try
        % 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);
            copy_variable_data_to_new_model(group_location_in_file,varname, xtype);
        end
    catch ME
        rethrow(ME)
    end

    % try to find groups in current level, if it doesn't work it's because there is nothing there
    try
        % 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
                %disp('found nested group!!')
                group_id = netcdf.inqNcid(group_location_in_file, netcdf.inqGrpName(group));
                %disp(netcdf.inqGrpNameFull(group_id))
                walk_nested_groups(group);
            end
        end
    catch
    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 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')
        % 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)))
            adress_to_attr = strrep(netcdf.inqGrpNameFull(group_location_in_file), '/', '.');
            
            data = netcdf.getVar(group_location_in_file, netcdf.inqVarID(group_location_in_file, varname));
            
    
            % matlab needs to know that ' ' = char()
            if xtype == 2 && isempty(all(data))
                data = cell(char());
            elseif numel(data) == 1 && xtype == 3 && data == -32767
                data = cell(char());
            end
            % band-aid for 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
            
            % 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', adress_to_attr, '.', varname, ' = ' , 'double(data);'];
                eval(arg_to_eval);
                %disp('saved int64 as int16')
            else
                arg_to_eval = ['model_copy', adress_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 e %e is an MException struct
            fprintf(1,'There was an error with %s! \n', varname)
            fprintf('The message was:\n%s\n',e.message);
            fprintf(1,'The identifier was:\n%s\n',e.identifier);
            disp()
            % more error handling...
        end
    end
end