Как построить вертикальные профили солёности и температуры воды на одном графике в 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');
Читать далее «Как построить вертикальные профили солёности и температуры воды на одном графике в MATLAB»

Notebook interface внутри MATLAB

Многие учёные, работающие с данными любят писать код не в виде голых скриптов, а используя так называемый notebook interface (типичный пример - Jupyter Notebook). Удобно это прежде всего потому, что ваш код становится более понятным, комментарии к коду хорошо читаются, а главное, прямо между строк можно вставлять результаты выполнения кода, например, полученные графики и построения.

В 2016 году MATLAB тоже создал свой вариант такого интерфейса, но на самом деле я пока не часто вижу, что им кто-то пользуется. Одна из причин (именно она побудила меня написать эту заметку) - не все знают об этой возможности MATLAB, хотя найти её совсем не сложно, она буквально постоянно маячит перед глазами (если у вас, конечно, установлен свежий MATLAB).

Итак, Notebook interface в MATLAB называется Live Editor, в которм можно создавать "живые скрипты" (Live Script). Создать такой скрипт можно также, как вы создаёте обычный скрипт (m-файл). На вкладке HOME для этого отведена специальная кнопка "new live script'. Если кнопки не окажется, поищите её под кнопкой "New":

Читать далее «Notebook interface внутри MATLAB»

Наносим линейный тренд на график временного хода в MATLAB

Сегодня мы нанесём линию тренда с доверительными интервалами на график временного ряда.

Работать мы будем с теми же данными (среднемесячные значения индекса NAO), с которыми работали в предыдущей заметке про добавление дат на график временного хода.

nao = load('norm.nao.monthly.b5001.current.ascii.txt');

Для того, чтобы тренд был явно выражен мы выбрали только последние несколько лет (да-да, это тот самый случай, когда учёные манипулируют с данными, никогда так не делайте в своей научной работе):

dates = datetime(nao(750:end,1),nao(750:end,2),15);
nao=nao(750:end,3);

Построим график для того, чтобы представлять, с чем имеем дело (параметр 'k' делает линию чёрной, а 'LineWidth' устанавливает толщину линии):

plot(dates, nao,'k','LineWidth',2)
Индекс NAO

Тренд здесь хорошо заметен. Если вам нужно удалить линейный тренд из данных, то сделать это можно при помощи простой функции detrend:

Читать далее «Наносим линейный тренд на график временного хода в MATLAB»

Добавление дат на график временного ряда в MATLAB

При работе с временными рядами часто требуется их визуализировать. Однако, для того, чтобы отложить по оси Х даты нужно немного специальных навыков работы с датами. Начнём мы с более старого сценария работы, который до сих пор может быть актуален, если у вас версия матлаба старее 2014 года. Также полезно будет его изучить, если вы будете разбираться в старых кодах.

Итак, наиболее классическим является вариант представления даты в виде, так называемого, serial date number. Например, 14 октября 1992 года в этой сиситеме будет равняться 727851:

datenum(1992, 10, 14)

ans =

      727851

Не сложно догадаться, что функция datenum как раз таки и переводит "человеческую" дату в этот формат. Само число при этом (в нашем случае число 727851) - это количество дней прошедших с первого января нулевого года ('Jan-1-0000 00:00:00').

Вообще говоря, вышеупомянутое обращение к функции datenum - не единственно возможное. Можно, например, так:

datenum('12-jun-17','dd-mmm-yy')

ans =

      736858

Аргумент 'dd-mmm-yy' даёт матлабу пояснение - в каком виде представлена наша дата. Видов этих может быть много, так год можно представить из четырёх цифр (yyyy) или из двух (yy), месяц записать полностью буквами (mmmm), трёхбуквенным сокращением (mmm), двумя цифрами (mm) или одной буквой (m) и т.д. Полный список вариантов можно найти в хэлпе к функции datenum. Иными словами, вы можете собирать дату как конструктор, главное рассказать матлабу, что вы имеете ввиду, например:

datenum('12--:06_+2017','dd--:mm_+yyyy')

ans =

      736858

Ну и на всякий случай скажем, что если вы хотите перевести дату обратно к текстовому виду, то сделать это можно так:

datestr(728447)

ans =

    '02-Jun-1994'

а если к раскидать дату в виде чисел и записать в матрицу, то так:

datevec(728447)

ans =

        1994           6           2           0           0           0

Давайте попробуем с этим поработать на конкретном примере. Для этого скачаем данные с ежемесячными значениями индекса NAO с сайта noaa, а именно Monthly mean NAO index since January 1950 - Ascii format for downloading.

Загрузим скачанные данные в MATLAB:

nao = load('norm.nao.monthly.b5001.current.ascii');
Первый столбец - годы, второй - месяцы, третий - значения индекса NAO

Сам временной ряд выглядит так:

Читать далее «Добавление дат на график временного ряда в MATLAB»

Как в MATLAB перевести координаты из формата «градусы-минуты» в градусы.

В прошлый раз мы с вами переводили долготы из формата [0 360] в формат [-180 180]. Сегодня мы разберёмся ещё с одной проблемой, которая периодически встречается в работе океанолога да и вообще любого географа. 

Широты и долготы принято измерять в градусах, в одном градусе 60 минут, а в одной минуте 60 секунд, секунды обычно записывают уже десятичной дробью, что намного удобнее. А ещё удобнее взять и сразу записать всё просто в градусах, а десятичной дробью записывать всё что меньше градуса и никаких минут и секунд.

Т.е. вместо 20°30' записывать 20.5°, или вместо 40°59'59'' записывать 40.9997°. Если вы не штурман, а датасайнтист, то в большинстве случаев так намного удобнее. Однако данные иногда приходят нам с координатами, записанными в разных форматах. Конечно, с пересчётом из одного формата в другой справится даже школьник, однако в MATLAB эти лишние телодвижения никчему, так как есть родная функция, которая сделает всё за вас. Вернее говоря, это несколько функций. Первая, dm2degrees переводит координаты из формата "градусы-минуты" в градлусы. Для её работы необходимо иметь матрицу, состоящую из двух столбцов, в первом столбце указаны градусы (только целые числа), во втором минуты (можно с десятыми долями), например:

lat = [50 34; 51 45.3; 53 21.89]

lat =

   50.0000   34.0000
   51.0000   45.3000
   53.0000   21.8900

В итоге, скормив эту матрицу функции dm2degrees, вы получите:

dm2degrees(lat)

ans =

   50.5667
   51.7550
   53.3648

Конечно, если у вас данные записаны не в одной переменной, а в двух разных, то проблем тоже не возникнет:

lat_d = [50; 51; 53]

lat_d =

    50
    51
    53

lat_m = [34; 45.3; 21.89]

lat_m =

   34.0000
   45.3000
   21.8900

dm2degrees([lat_d lat_m])

ans =

   50.5667
   51.7550
   53.3648

Если ваши данные представлены в формате "градусы - минуты - секунды", то вам нужна будет другая функция dms2degrees. Работает она практически также, только на вход ей нужна теперь матрица состоящая из трёх столбцов (градусы, минуты, секунды). Грубо говоря, вам нужно будет записать dm2degrees([lat_d, lat_m, lat_s]).

Есть и обратные функции degrees2dm и degrees2dms, которые конвертируют всё обратную сторону.

 telegram

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

Как конвертировать долготы из формата [0 360] в [-180 180] в MATLAB

В сегодняшней заметке мы переведём долготы, записанные в формате [-180 180] в формат [0 360] и обратно.

Как известно, долготы обычно принято отсчитывать от нулевого меридиана. На запад от него - западная широта, на восток - восточная. Если буквенные обозначения использовать неудобно (а их при работе с данными неудобно использовать практически всегда), то восточные долготы записываются как положительные числа, а западные - как отрицательные. Такой формат записи будем обозначать [-180 180].

Однако иногда отрицательных значений избегают и долготы отсчитывают не на запад и восток от нулевого меридиана, а только на восток. Такой формат географических координат будем называть [0 360] и бывают случаи, когда данные приходят именно нём. Что же делать, если вам хочется перевести формат [0 360] в [-180 180]? Процедура несложная:

lon_new=rem((lon+180),360)-180;

где rem - функция, которая находит остаток от деления (lon+180)/360. Такая функция есть во многих языках программирования, но конкретно этот пример записан для MATLAB.

Кроме того, в MATLAB есть специальная функция wrapTo180, которая избавит вас от написания каких либо математических операций. Например,

wrapTo180([0 50 175 195 355])

ans =

     0    50   175  -165    -5

Если же вам хочется конвертировать координаты в формат [0 360], то есть аналогичная функция - wrapTo360:

wrapTo360([-10 -5 5 180 185])

ans =

   350   355     5   180   185

Бываю и более непривычные форматы, когда избегают всяческих резких переходов в долготах, и градусы отсчитываются от минус бесконечности до плюс бесконечности, т.е. после 359 идёт не 0, а 360, 361 и т.д. Данная функция с лёгкостью справляется и с этой задачей:

wrapTo360([-4000 -500 -5 89 567])

ans =

   320   220   355    89   207
wrapTo180([-4000 -500 -5 89 567])

ans =

   -40  -140    -5    89  -153

Ну и напоследок следует сказать, что при картировании ваших данных может случиться, что вам нужно, чтобы долготы записывались в формате [0 360] (по умолчанию на картах долготы записываются в [-180 190]), то делается при помощи функции mlabelzero22pi.

 telegram

В следующий решим ещё одну проблему, связанную с форматом записи координат. Чтобы не пропустить новые материалы, подпишитесь на канал в Telegram: https://t.me/koldunovaleksey

 

Как визуализировать океанологические данные на карте в MATLAB

Сегодня мы рассмотрим несколько готовых рецептов для визуализации океанологических данных. Эта заметка именно с рецептами и не претендует на полноценное руководство, которое для Матлаба существует и достаточно объёмное (около 1000 страниц). Скачать его сегодня (в ноябре 2018 года) можно по ссылке.

В прошлой заметке мы с вами открыли netCDF файл в MATLAB, теперь настало время визуализировать данные на карте. Если в прошлый раз мы выдумали гипотетическое название файла, то сегодня откроем самые настоящие данные. Для примера скачаем файл с температурой поверхности океана (SST) с сайта https://oceancolor.gsfc.nasa.gov/

В нашем файле содержатся среднемесячные значения температуры поверхности воды за сентябрь 2018 года, если хотите точно такой же, то вот ссылка.

Теперь настало время выудить из этого файла данные о температуре:

sst = ncread('A20182442018273.L3m_MO_SST_sst_9km.nc','sst');

Полученную переменную sst легко визуализировать простыми средствами (как вы бы визуализировали любую другую матрицу), например через функцию imagesc:

imagesc(sst)
визуализация данных в MATLAB

Если вас смущает, что матрица "лежит на боку", то можете переориентировать её, написав в командной строке:  

Читать далее «Как визуализировать океанологические данные на карте в MATLAB»

Как открыть netCDF файл в MATLAB

В данной заметке я покажу два способа открыть NetCDF-файл в MATLAB, один способ более продвинутый и более сложный, а второй более простой. 

Несколько лет назад брат написал статью в своём блоге koldunov.net, где показал, как открыть в MATLAB файлы формата netCDF.

Он использовал функцию netcdf и всю процедуру необходимо было выполнять в несколько шагов. Не буду полностью повторять его статью, но вкратце напомню примерный порядок действий.

Для начала нужно открыть netCDF файл, в результате чего в переменную ncid будет записано число, которое будет являеться неким идентификатором, при помощи которого мы сможем вдальнейшем обращаться к нашему файлу:

ncid = netcdf.open('temperature_and_salinity_2018.nc','NC_NOWRITE');

Зная название (информацию о файле можно получить при помощи функции ncdisp) нужной вам переменной (например, температуры) внутри netCDF-файла, определить её ID (который на самом деле является просто порядковым номером переменной внутри файла; счёт начинается с нуля)

varid = netcdf.inqVarID(ncid,'temperature')
varid =

     3

Затем выудить по полученному ID данные для нужной нам характеристики:

data = netcdf.getVar(ncid,3);

Вот и всё, мы открыли наши данные. Но этих процедур часто бывает недостаточно, потому как для уменьшения размеров netCDF файлов прибегают к различного рода хитростям, например приводят все данные к целочисленным значениям и вычитают какое-нибудь число, если все данные больше этого самого числа.

Читать далее «Как открыть netCDF файл в MATLAB»

Построение точек на карте по заданным координатам при помощи Basemap в Python.

Видео аватар

Если внимательно посмотреть на последние заметки нашего блога, то можно подумать, что он превратился в чисто фотографический. Однако недавно мы хорошенько обдумали перспективы развитя наших сайтов, не побоялись и поставили всё с ног на голову. Поскольку фотографическая тематика постов отныне будет жить только на koldunov.com, этот сайт вновь станет выглядеть так, как он выглядел в самом начале. Да-да, начиналось всё не с фотографических обзоров, а с заметок, типа:
Вязание крючком игрушки Бендера (Футурама)
Взрослая загадка на смекалку — расшифровать послание инопланетянам
Работа с GMT
Изготовление трафарета
и т.д.

Сегодня мы вернёмся к теме картирования. Задача стоит следующая.
В отдельных точках акватории проводились измерения океанологических хараткеристик. Для каждой точки нам известны географические координаты (широта и долгота). Необходимо нанести в виде кружков положение точек на карте. Делать всё мы будем в Python (я обычно при этом использую ipython notebook, так удобнее).

Для начала импортируем всё, что нам понадобится:

%pylab inline
import matplotlib.pylab as plt
from mpl_toolkits.basemap import Basemap

затем: создадим две переменные с широтами и долготами. В данном случае проще и быстрее было создать их вручную:

lat = [66.297,66.299,66.298,66.295,66.301,66.304,66.288,66.289,66.286,66.289]
lon = [33.640,33.660,33.690,33.747,33.829,33.908,33.891,33.839,33.781,33.740]

Ну и теперь само построение карты:

m = Basemap(projection='merc',llcrnrlat=66.22,urcrnrlat=66.37,\
            llcrnrlon=33.60,urcrnrlon=34,resolution='f')
figsize(10,15)

x, y = m(lon,lat)

m.drawcoastlines()
m.fillcontinents(color='gray',lake_color='white')
m.drawmapboundary(fill_color='white')
m.drawparallels(np.arange(66.22,66.37,.04),labels=[1,0,0,0], fontsize =14)
m.drawmeridians(np.arange(33.60,34.,.1),labels=[0,0,0,1], fontsize =14)

m.scatter(x,y,20,marker='o',color='k')
plt.title("Location of the measurement points", fontsize =14)
plt.show()

plot_points
Некоторые пояснения к тому, что мы сделали (я не стал это делать в коде, чтобы не перегружать его): Читать далее «Построение точек на карте по заданным координатам при помощи Basemap в Python.»