% Очищаем командную строку (Command Window)
clc
% Очищаем Workspace (удаляем все переменные)
clear all

%% Скачиваем себе на компьютер данные реанализа, которые лежат на сервере (ничего не меняем!)
% Данные скачичаются в ту папку, которая открыта в окне слева (Current Folder)
% Проделать эту операцию нужно один раз, потом можно просто закомментировать этот блок
furl ='http://kibel.srcc.msu.ru:8080/share.cgi?ssid=5743fc3a82a14733a325cf8c68398db2&fid=5743fc3a82a14733a325cf8c68398db2&filename=surface_t2m_precip.nc&openfolder=forcedownload&ep='
fname = 'p-levels_t_u_v_geopotential.nc'
fpath= websave(fname, furl);

%% Открываем данные
% Прописываем путь к файлу
% Знак точка с запятой ";" в конце строки блокирует вывод полученного в ходе 
% выполнения команды результата в командной строке
filepath = "C:\";
% Прописываем название файла
filename = "p-levels_t_u_v_geopotential.nc";
% Просматриваем содержимое файла NetCDF командой ncdisp (переменные, их размерность и др.)
% Информация отображается в командной строке 
ncdisp (filepath + filename);

% Выбираем необходимые нам переменные
% Загружаем в Workspace u- и v- компоненты скорости ветра  командой ncread
% В скобках указываем название переменной - то, как она называется в файле
% В Workspace должны появиться переменные u_data и v_data в виде 4-х мерного массива
u_data = ncread(filepath + filename, 'u');
v_data = ncread(filepath + filename, 'v');

% Можно посмотреть размер массива командой size и вывести в переменную u_dimensions
% Размер массива u_data: 11440*721*3*3 (долготы*широты*барические уровни*время)
u_dimensions = size(u_data);

% Аналогичным образом загружаем векторы с долготами, широтами и временем
% А также барические уровни (700, 850 и 1000 гПа)
lonn = ncread(filepath + filename, 'longitude');
latt = ncread(filepath + filename, 'latitude');
time = ncread(filepath + filename, 'time');
p_level = ncread(filepath + filename, 'level');

%% Преобразовываем данные
% Вычисляем полный вектор скорости через u- и v- компоненты
wind = sqrt(u_data .^ 2 + v_data .^ 2);
% Находим среднее значение скорости ветра за три летних месяца командой
% mean. В скобках указывается, к какой размерности массива применяется
% операция. В нашем случае для времени - это 4. 
wind_ave_time = mean(wind, 4);

% Данные о времени даны в количестве часов с 00:00 1 января 1900 г. -> переводим их в понятные нам даты
% В Matlab'е есть команда datenum для работы с данными формата даты
% Она переводит заданную дату в количество дней с 00:00 1 января 0000 г.
% создадим переменную t0 - количество дней с 00:00 1 января 0000 г. по 00:00 1 января 1900 г.
t0 = datenum('1900-01-01 00:00:00', 'yyyy-mm-dd HH:MM:SS');
% переведем даты из нашего файла из формата "часы с 00:00 1 января 1900 г."
% в "дни с с 00:00 1 января 1900 г."
% time_days_since1900year - количество дней с 00:00 1 января 1900 г.по
% даты, указанные в файле. Точка перед знаком деления обозначает, что
% данная арифметическая операция применяется отдельно к каждому члену
% массива
time_days_since1900year = time ./ 24;
% сложим две переменные между собой, чтобы получить time_days_since0year - количество дней с 00:00 1 января 0000 г. 
% по даты, указанные в файле
time_days_since0year = time_days_since1900year + t0;
% конвертируем с помощью команды datetime из формата "дни с с 00:00 1 января 0000 г." в нормальный вид, 
% указав при этом, что стартуем именно с 0000 года флагом 'datenum'
% если нажать на переменную time_real, можно увидеть, за какие даты даны
% данные
time_I_understand = datetime(time_days_since0year, 'ConvertFrom', 'datenum');

%% Отрисовка данных на карте
% Переворачиваем командой flipud вектор широт, чтобы они начинались с -90
latt1 = flipud(latt);
% Строим координатную сетку с размерностью долгот и широт
[lon,lat] = meshgrid(lonn, latt1);

% Открываем окно для нашего рисунка во весь экран
figure('units','normalized','outerposition',[0 0 1 1])
% Указываем проекцию и координаты
m_proj ('Equidistant Cylindrical','lon',[0 360], 'lat', [-90 90]);
% Выбираем переменную на первом барическом уровне (700 гПа) 
wind_one_lev = wind_ave_time(:,:,1);
% Транспонируем (знак апострофа "'") -> меняем местами
% размерности), так как сетка имеет размерность широты*долготы
wind_transp = wind_one_lev';
% Отрисовываем переменную, заодно переворачивая её командой flipud, чтобы
% данные начинались с 90° ю.ш.
m_contourf(lon, lat, flipud(wind_transp),'LineColor','none', 'LevelStep', 1);
% Отрисовываем линию побережья
m_coast('Color', 'k','LineWidth', 1)
% Наносим линии сетки
m_grid
% Обозначаем цветовую шкалу
colormap(turbo)
% Выводим цветовую шкалу на рисунок
c = colorbar;
% Подписываем шкалу
c.Label.String = 'Wind, m/s';
% Поворачиваем её на -90 градусов
c.Label.Rotation = -90;
% Сдвигаем её, чтобы она не заезжала на шкалу после поворота
c.Label.Position = [3.5 -12.5 0];
% Определяем размер шрифта подписи
c.Label.FontSize = 12;
% Ограничиываем размах значений шкалы
clim([0 25])
% Добавляем название рисунка
% Команада datestr переводит значение переменной из формата даты в формат
% строки
% По умолчанию название будет написано жирным шрифтом -> команда 'FontWeight','Normal'
% делает подпись обычной
title("Mean Wind Speed" + " at " + int2str(p_level(1)) + "hPa; " + datestr(time_I_understand(1), 'mmmm') + " - " + datestr(time_I_understand(3), 'mmmm, yyyy') ,'FontWeight','Normal','FontSize',16)

%% Отрисовка данных на карте циклом для всех высот и сохранение рисунков
% Цикл от 1 до длины вектора с точками по времени, то есть от 1 до 3 с
% шагом 1
for i = 1:length(p_level)
    % Открываем окно для нашего рисунка во весь экран
    figure('units','normalized','outerposition',[0 0 1 1])
    % Указываем проекцию и координаты
    m_proj ('Equidistant Cylindrical','lon',[0 360], 'lat', [-90 90]);
    % Выбираем переменную на первом барическом уровне (700 гПа) 
    wind_one_lev = wind_ave_time(:,:,i);
    % Транспонируем (знак апострофа "'") -> меняем местами
    % размерности), так как сетка имеет размерность широты*долготы
    wind_transp = wind_one_lev';
    % Отрисовываем переменную, заодно переворачивая её командой flipud, чтобы
    % данные начинались с 90° ю.ш.
    m_contourf(lon, lat, flipud(wind_transp),'LineColor','none', 'LevelStep', 1);
    % Отрисовываем линию побережья
    m_coast('Color', 'k','LineWidth', 1)
    % Наносим линии сетки
    m_grid
    % Обозначаем цветовую шкалу
    colormap(turbo)
    % Выводим цветовую шкалу на рисунок
    c = colorbar;
    % Подписываем шкалу
    c.Label.String = 'Wind, m/s';
    % Поворачиваем её на -90 градусов
    c.Label.Rotation = -90;
    % Сдвигаем её, чтобы она не заезжала на шкалу после поворота
    c.Label.Position = [3.5 -12.5 0];
    % Определяем размер шрифта подписи
    c.Label.FontSize = 12;
    % Ограничиываем размах значений шкалы
    clim([0 25])
    % Добавляем название рисунка
    % Команада datestr переводит значение переменной из формата даты в формат
    % строки
    % По умолчанию название будет написано жирным шрифтом -> команда 'FontWeight','Normal'
    % делает подпись обычной
    title("Mean Wind Speed" + " at " + int2str(p_level(i)) + "hPa; " + datestr(time_I_understand(1), 'mmmm') + " - " + datestr(time_I_understand(3), 'mmmm, yyyy') ,'FontWeight','Normal','FontSize',16)
    
    % Сохраняем рисунок в формате png с разрешением 300 dpi
    print("Mean_Wind_Speed_at" + int2str(p_level(i)) + "hPa_" + datestr(time_I_understand(1), 'mmmm') + "-" + datestr(time_I_understand(3), 'mmmm_yyyy'), '-dpng', '-r300');
end


%% Строим вертикальный разрез скорости ветра, осредненной по всем широтам, долготам и за весь временной период
% По времени мы уже осредняли, необходимо повторить операцию осреднения применительно к долготам и широтам
% Осредняем по всем долготам. Теперь у нас не 1440 значений по долготе, а 1. Размер массива 1*721*3.
% Команда squeeze убирает размерность, которая равна 1 => массив становится размером 721*3 
wind_ave_lon = squeeze(mean(wind_ave_time, 1));
% Аналогично осредняем по всем широтам
wind_ave_lat = squeeze(mean(wind_ave_lon, 1));

%% Отрисовывем полученные значения на графике командой plot
% Указываем, что будет по оси X - скорость ветра, и что по оси Y -
% барические уровни
plot(wind_ave_lat, p_level, 'Color', [190/256 0 0], 'LineWidth', 2)
% Переворачиваем ось Y, чтобы значение барических уровней шли по убыванию
% (от 1000 к 700 гПа)
set(gca, 'YDir','reverse')
% Подписываем оси
xlabel('Wind speed, m/s')
ylabel('Pressure level, hPa')
% Накладываем сетку
grid on
% Ограничиваем размах значений по оси X
xlim([0 7])
% Добавляем название к графику
title("Mean Wind Speed for " + datestr(time_I_understand(1), 'mmmm') + " - " + datestr(time_I_understand(3), 'mmmm, yyyy') ,'FontWeight','Normal','FontSize',12)

% Сохраняем график
print("Vertical_plot_Mean_Wind_Speed_for_" + datestr(time_I_understand(1), 'mmmm') + "-" + datestr(time_I_understand(3), 'mmmm_yyyy'), '-dpng', '-r300');