NetCDF Extraction script

ExtractData.m

% DATA_EXTRACTOR is a piece of program to extract data from netCDF file
% with CF convention. 
% DATA_EXTRACTOR takes the location of the netCDF file graphically and then
% ask the user for bounding box of the data extraction. It also shows basic
% information from first file in the directory to input the to be extracted
% variable name. 
% When giving input Numbers are input as it is. However variable name must
% be enclosed in a 'singleQuote'.
% Author: Jamal Uddin Khan
% IWFM, BUET - 2015.

function [status] = ExtractData(inputdir, outputdir, outvarname)  
% Listing the files from the directory.
files = dir(fillfile(inputdir,'*.nc'));

% Showing the ncdisp of typical file
if ~isempty(files)
    ncdisp(files(1).name);
    % Loading Common values
    % Lat-lon is taken same for all values.
    file = fullfile(inputdir, files(1).name);
    lon = ncread(file, 'lon');
    lat = ncread(file, 'lat');
end

% Start Date
% It can be found from the time definition of the input files.
start_date = datenum('1949-12-1');

% lat-lon range of interest
lon_w = input('Enter longitude of the west end : ');
lon_e = input('Enter longitude of the east end : ');
lat_s = input('Enter latitude of the south end : ');
lat_n = input('Enter latitude of the north end : ');

% Here FindClosest finds the closest lat and lon that is actually in the
% file.
[lat_s, ~] = FindClosest(lat, lat_s);
[lat_n, ~] = FindClosest(lat, lat_n);
[lon_w, ~] = FindClosest(lon, lon_w);
[lon_e, ~] = FindClosest(lon, lon_e);

% Grid size is determined by subtracting the consequtive grid. 
% The grid is takes as structured and equally spaced.
grid_size = abs(lon(1) - lon(2));

% It is the variable name that is to be extracted. 
% The input must be a string. That is in single quote 'varName'
var_name = input('Enter varname to be extracted : ');

% If data needs any multiplication change here
% It is for conversion of one unit to another.
multiplier = input('Enter the multiplier if any : '); 
addition = input('Enter the value to be added : ');

% creating output file name & open file for writing
outdir = uigetdir();
outfile = strsplit(files(1).name, '_');
outfile = outfile(1: length(outfile)-2);
outfile = cell2mat(outfile);
outfile = [outdir, '/', outfile '.csv']; 
disp(['Files will be saved as - ', outfile]);

tic; % Keeping track of time.
% Create listing for lat-lon of interest
position = zeros(2, length(lon_w : grid_size : lon_e) * length(lat_s : grid_size : lat_n));
ly = lat_s : grid_size : lat_n;
lx = lon_w : grid_size : lon_e;
for i = 1 : length(ly)
    position(1, (i - 1)*length(lx) +1 : i * length(lx)) = ones(1, length(lx)) * ly(i);
    position(2, (i - 1)*length(lx) +1 : i * length(lx)) = lx;
end

% Create headers
headers = num2cell(position);
headers = [{'latitude';'longitude'} headers];
[row, col] = size(headers);

% Creating the format specifier
day_string = {'%s,'};
num_form = {'%f,'};
end_elem = {'%f\n'};
string_format = [day_string repelem(num_form, length(position)-1) end_elem];
string_format = cell2mat(string_format);

% Opening file for writing header in write mode
fid = fopen(outfile, 'w');
% Now writing header line by line using fprintf
for i = 1 : row
    fprintf(fid, string_format, headers{i, :});
end
fclose(fid);

% Opening file for writing data values
% This file will be closed at the end
fid = fopen(outfile, 'a');

% File iteration will be started here
% Iterate for file in files
for file_num = 1 : length(files)
    % Setting file name
    file = files(file_num).name;

    % Variables are time, var_name, lon, lat
    time = ncread(file, 'time'); 
    var_value = ncread(file, var_name); 
    lon = ncread(file, 'lon'); 
    lat = ncread(file, 'lat'); 

    % Temporary data storage 
    year_chunk = zeros(length(time), length(lon_w : grid_size : lon_e) * length(lat_s : grid_size : lat_n) + 1);
    year_chunk = num2cell(year_chunk);
    % Running loop to extract data
    dates = cellfun(@datestr, num2cell(time + start_date), 'UniformOutput', false);
    year_chunk(:, 1) = dates;
    for i = 1 : length(position)
        % pos is in latitude; longitude format in ith column
        varData = var_value(find(lon == position(2, i)), find(lat == position(1, i)), :) * multiplier + addition;
        % flattening prData
        varData = reshape(varData(1:length(time)), [length(time), 1]);
        year_chunk(:, i+1) = num2cell(varData);
    end

    % writing year_chunk
    [rows, cols] = size(year_chunk);
    for row = 1 : rows
        fprintf(fid, string_format, year_chunk{row, :});
    end
    
    % displaying completion message
    msg_comp = ['File ',num2str(file_num), ' of ', num2str(length(files)), ' - ', file, ' - ', 'Completed!'];
    disp(msg_comp);
end
fclose('all'); % closing all open files