Как построить вертикальные профили солёности и температуры воды на одном графике в MATLAB
Сегодня мы построим графики распределения значений солёности и температуры воды по глубине.
Работать мы будем с данными World Ocean Atlas 2013 V2. Вместо того, чтобы скачивать .nc файлы целиком, для удобства загрузим с сервера только нужные нам переменные. Делается это при помощи протокола OPENDAP.
Нам потребуются пять переменных: температура (запишем в переменную t_global), солёность (s_global), широта (lat), долгота (lon) и глубина (z). Запускаем на закачку прямо в MATLAB:
t_global = ncread('https://data.nodc.noaa.gov/thredds/dodsC/woa/WOA13/DATAv2/temperature/netcdf/decav/1.00/woa13_decav_t00_01v2.nc','t_an');
lat = ncread('https://data.nodc.noaa.gov/thredds/dodsC/woa/WOA13/DATAv2/temperature/netcdf/decav/1.00/woa13_decav_t00_01v2.nc','lat');
lon = ncread('https://data.nodc.noaa.gov/thredds/dodsC/woa/WOA13/DATAv2/temperature/netcdf/decav/1.00/woa13_decav_t00_01v2.nc','lon');
z = ncread('https://data.nodc.noaa.gov/thredds/dodsC/woa/WOA13/DATAv2/temperature/netcdf/decav/1.00/woa13_decav_t00_01v2.nc','depth');
s_global = ncread('https://data.nodc.noaa.gov/thredds/dodsC/woa/WOA13/DATAv2/salinity/netcdf/decav/1.00/woa13_decav_s00_01v2.nc', 's_an');
Итак, мы загрузили необходимые данные. В переменных t_global и s_global, содержатся данные о среднемноголетнем распределении температуры и солёности в Мировом океане на разных глубинах. Нам потребуется только одна точка на карте, в которой мы будем строить изменчивость характеристик с глубиной. Выберем любую, на свой вкус:
t = squeeze(t_global(lon==-20.5,lat==85.5,:));
s = squeeze(s_global(lon==-20.5,lat==85.5,:));
В океанологии принято строить профили распределения характеристик так, чтобы глубина была отложена по оси Y, поэтому в функции plot на первое место ставим t (или s), а на второе z: plot(t,z). Сразу давайте договоримся, что мы дополнительно будем увеличивать толщину линии (‘linewidth’) до значения 2, иначе линии как-то совсем тонкие получатся и нам их плохо будет видно. Итого:
plot(t,z,'linewidth',2)
приведёт к такому результату:

Следует также отметить, что океанологу привычнее смотреть на подобные графики если ноль глубин находится вверху, т.е. нужно наш профиль перевернуть. Сделать это можно добавив знак минуса к переменной z при построении:
plot(t,-z,'linewidth',2)

Напомню, что перед нами стоит задача сравнить распределение по глубине значений температуры и солёности воды. Один из способов это сделать — построить рядом с графиком температуры график распределения солёности, используя функцию subplot:
subplot 121
plot(t,-z,'linewidth',2)
ylabel 'Depth (m)'
title 'Temperature ({\circ}C)'
subplot 122
plot(s,-z,'linewidth',2)
title 'Salinity'

Однако, есть ещё один классический способ отображения профилей температуры и солёности, который часто можно встретить в литературе, а именно — поместить обе линии на одном графике. Самый простой путь — сделать это через hold on:
plot(t,-z,'linewidth',2)
hold on
plot(s,-z,'linewidth',2)

Однако, как видите, результат в нашем случае будет совершенно негодным. Такое может «прокатить» только в редких случаях, когда значения температуры и солёности воды очень близки. Значит нам нужно построить два графика так, чтобы на нём присутствовали две независимые оси X (обычно, одну располагают сверху, а вторую снизу).
Итак, для начала построим график температуры:
plot(t,-z,'r','linewidth',2)
Затем при помощи get запишем в переменную gp информацию о том, где расположены оси этого графика:
gp = get(gca,'pos');
Теперь создаём новые оси в том же самом месте (‘pos’,gp), делаем их бесцветными (‘color’,’none’) и располагаем подписи наверху (‘XAxisLocation’,’top’):
axes('pos',gp,'color','none','XAxisLocation','top');
Вводим команду hold on, чтобы последующее построение строилось в указанных осях и рисуем распределение солёности:
plot(point1_s,-z,'b','linewidth',2)
Итоговый код целиком и результат будут выглядеть так:
plot(t,-z,'r','linewidth',2)
gp = get(gca,'pos');
axes('pos',gp,'color','none','XAxisLocation','top');
hold on
plot(s,-z,'b','linewidth',2)

Добавим подписи и зададим разный цвет для верхней и нижней оси X:
plot(t,-z,'r','linewidth',2)
ylabel 'Depth (m)'
xlabel 'Temperature ({\circ}C)'
set(gca,'XColor','r')
gp = get(gca,'pos');
axes('pos',gp,'color','none','XAxisLocation','top','XColor','b');
hold on
plot(s,-z,'b','linewidth',2)
xlabel('Salinity')

Если вы внимательно посмотрите на получившуюся картинку, то заметете, что на верхней синей шкале оси X (шкала солёности) видны красные дублирующие отметки нижней шкалы (температуры)

Убрать это можно, выключив контур, которым по умолчанию очерчены графики (box off):
plot(point1_t,-z,'r','linewidth',2)
ylabel 'Depth (m)'
xlabel 'Temperature ({\circ}C)'
set(gca,'XColor','r','box','off')
gp = get(gca,'pos');
axes('pos',gp,'color','none','XAxisLocation','top','XColor','b');
hold on
plot(point1_s,-z,'b','linewidth',2)
xlabel('Salinity')

Иногда бывает важно получить возможность оценить по графику значения построенных характеристик, в таком случае удобно добавить на картинку сетку (grid on). Однако, поскольку разбивка шкал температуры и солёности не совпадает, то придётся либо строить две разные сетки:
plot(t,-z,'r','linewidth',2)
grid on
ylabel 'Depth (m)'
xlabel 'Temperature ({\circ}C)'
set(gca,'XColor','r','box','off')
gp = get(gca,'pos');
axes('pos',gp,'color','none','XAxisLocation','top','XColor','b');
hold on
plot(s,-z,'b','linewidth',2)
grid on
xlabel('Salinity')

либо подобрать пределы диапазонов (‘XLim’) температуры и солёности на графиках так, чтобы шкалы были разделены на одинаковое количество частей. Цвет сетки при этом сделаем чёрным (set(gca,’GridColor’,’k’)):
plot(t,-z,'r','linewidth',2)
grid on
ylabel 'Depth (m)'
xlabel 'Temperature ({\circ}C)'
set(gca,'XColor','r','box','off','XLim',[-2.5 2.5])
gp = get(gca,'pos');
axes('pos',gp,'color','none','XAxisLocation','top',...
'XColor','b','XLim',[31 36]);
hold on
plot(s,-z,'b','linewidth',2)
grid on
set(gca,'GridColor','k')
xlabel('Salinity')

Ну и последний штрих. Есть мнение, что строить такие профили линиями — не совсем корректно. Давайте модернизируем нашу картинку так, чтобы вместо линий значения температуры и солёности обозначались точками. Для этого поменяем ‘r’ на ‘.r’ и ‘b’ на ‘.b’:
plot(t,-z,'.r','linewidth',2)
grid on
ylabel 'Depth (m)'
xlabel 'Temperature ({\circ}C)'
set(gca,'XColor','r','box','off','XLim',[-2.5 2.5])
gp = get(gca,'pos');
axes('pos',gp,'color','none','XAxisLocation','top',...
'XColor','b','XLim',[31 36]);
hold on
plot(s,-z,'.b','linewidth',2)
grid on
set(gca,'GridColor','k')
xlabel('Salinity')


Чтобы не пропустить новые материалы с рецептами по работе с океанологическими данными, подпишитесь на канал в Telegram: https://t.me/koldunovaleksey