Commit ab9df771 authored by Willi Rath's avatar Willi Rath
Browse files

Merge branch '3-add-SLTAC-matlab' into 'master'

Resolve "Add Matlab Example for SLTAC"

Closes #3

See merge request !8
parents d4a70a34 68fc714b
Loading
Loading
Loading
Loading
+2 −0
Original line number Diff line number Diff line
.ipynb_checkpoints/
*.nc
*.mat
+96 −0
Original line number Diff line number Diff line
% calculate monthly mean surface geostrophic velocity off Brazil from 2000-2017

clear all
close all

% define path to data
base_path = '/data/c2/TMdata/git_geomar_de_data/SLTAC_GLO_PHY_L4_REP/';
SLTAC_version = 'v1.x.x';
path_to_data = [base_path, SLTAC_version, '/data/'];

% find indices of the region
lat_lim = [-12 -4];
lon_lim = [-40 -30] + 360;
lat_nc = -89.875 : 0.25 : 89.875;
lon_nc = 0.125 : 0.25 : 359.875;
idx = find(lon_nc > lon_lim(1) & lon_nc < lon_lim(2));
idy = find(lat_nc > lat_lim(1) & lat_nc < lat_lim(2));

% choose time period
year_folder = [2000 : 2017];

% allocate monthly fields
u_geo_m = NaN(length(idx), length(idy), length(year_folder) * 12);
v_geo_m = u_geo_m;
time_m = NaN(length(year_folder) * 12, 1);

% initialize progressbar
h = waitbar(0, 'Loading and averaging ...');

% load all files of one month and average
for iy = 1 : length(year_folder)

    % find all data files for current year (determined by iy)
    filename_year = [path_to_data, num2str(year_folder(iy))];
    KK = dir(filename_year);
    KK = KK(3 : end);

    % for all file names, extract months from file name
    for id = 1 : length(KK)
        mm(id) = str2num(KK(id).name(29 : 30));
    end

    % for all months, find files for this month and average
    for im = 1:12
        kk = find(mm == im); % returns indices matching current month

        % allocate temporary data fields
        u_geo = NaN(length(idx),length(idy),length(kk));
        v_geo = u_geo;

        % load all files for current month
        for ik = 1:length(kk)
            filename_current_day = [filename_year, '/', KK(kk(ik)).name];

            time_vec(ik) = double(ncread(filename_current_day, 'time'));
            time_vec(ik) = time_vec(ik) + datenum(1950, 1, 1, 0, 0, 0);

            ugos = double(ncread(filename_current_day, 'ugos'));
            u_geo(:,:,ik) = ugos(idx, idy);

            vgos = double(ncread(filename_current_day, 'vgos'));
            v_geo(:,:,ik) = vgos(idx, idy);
        end

        % average and store in final arrays
        u_geo_m(:, :, (iy - 1) * 12 + im) = nanmean(u_geo, 3);
        v_geo_m(:, :, (iy - 1) * 12 + im) = nanmean(v_geo, 3);
        time_m((iy - 1) * 12 + im) = nanmean(time_vec);

        % update status bar
        waitbar(((iy - 1) * 12 + im) / (length(year_folder) * 12), h)
    end
    clear mm KK
end

% close status bar
close(h)

% save as matfile
lon = lon_nc(idx) - 360;
lat = lat_nc(idy);
time_vec = time_m;
ugeo_mon = u_geo_m;
vgeo_mon = v_geo_m;
save example_07_SLTAC_monthly_mean_uv_off_Brazil lon lat time_vec ugeo_mon vgeo_mon

% and plot
figure(1);
m_proj('mercator', 'long', [-40 -30], 'lat', [-12 -4]);
hold on
m_contourf( ...
    lon_nc(idx) - 360, lat_nc(idy), nanmean(u_geo_m, 3)', ...
    [-.5 : 1 / 63 : .5], 'edgecolor', 'none')
m_gshhs_h('patch', [.5 .5 .5]);
colorbar
saveas(gcf, 'example_07_SLTAC_monthly_mean_uv_off_Brazil.png', 'png')
+45.5 KiB
Loading image diff...