Для того, чтобы тренд был явно выражен мы выбрали только последние несколько лет (да-да, это тот самый случай, когда учёные манипулируют с данными, никогда так не делайте в своей научной работе):
При работе с временными рядами часто требуется их визуализировать. Однако, для того, чтобы отложить по оси Х даты нужно немного специальных навыков работы с датами. Начнём мы с более старого сценария работы, который до сих пор может быть актуален, если у вас версия матлаба старее 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 с сайта , а именно Monthly mean NAO index since January 1950 - .
Широты и долготы принято измерять в градусах, в одном градусе 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
Конечно, если у вас данные записаны не в одной переменной, а в двух разных, то проблем тоже не возникнет:
Если ваши данные представлены в формате "градусы - минуты - секунды", то вам нужна будет другая функция dms2degrees. Работает она практически также, только на вход ей нужна теперь матрица состоящая из трёх столбцов (градусы, минуты, секунды). Грубо говоря, вам нужно будет записать dm2degrees([lat_d, lat_m, lat_s]).
Есть и обратные функции degrees2dm и degrees2dms, которые конвертируют всё обратную сторону.
Чтобы не пропустить новые материалы с рецептами по работе с океанологическими данными, подпишитесь на канал в Telegram:
В сегодняшней заметке мы переведём долготы, записанные в формате [-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, которая избавит вас от написания каких либо математических операций. Например,
Бываю и более непривычные форматы, когда избегают всяческих резких переходов в долготах, и градусы отсчитываются от минус бесконечности до плюс бесконечности, т.е. после 359 идёт не 0, а 360, 361 и т.д. Данная функция с лёгкостью справляется и с этой задачей:
Ну и напоследок следует сказать, что при картировании ваших данных может случиться, что вам нужно, чтобы долготы записывались в формате [0 360] (по умолчанию на картах долготы записываются в [-180 190]), то делается при помощи функции mlabelzero22pi.
В следующий решим ещё одну проблему, связанную с форматом записи координат. Чтобы не пропустить новые материалы, подпишитесь на канал в Telegram:
Сегодня мы поговорим о том, какую цветовую шкалу выбрать для отображения океанологических данных.
Я неоднократно натыкался в сети на негативные отзывы по поводу радужных цветовых шкал, к коим относится и цветовая шкала jet (colormap jet), которая до 2014 года использовалась в MATLAB по умолчанию для построения различного рода визуализаций.
Например, рассказывается, как радужная шкала может ввести вас и ваших читателей в заблуждение (особенно это касается интерпретации градиентов визуализируемых характеристик). На самом деле в интернете ещё много подобного рода заметок и даже статей в научных журналах. Не случайно в 2014 году в MATLAB решили изменить colormap, используемый по умолчанию с jet на parula, озвучив в качестве объяснения .
Поводом для моей сегодняшней заметки стал "заказ" коллеги, который, читая материалы koldunov.ru попросил рассказать, как можно использовать цветовую шкалу, отличную от той, что предлагает MATAB. Он признался, что ему не нравится ни одна матлабовская шкала, кроме разве что jet, да и то с натяжкой. Мне, честно говоря тоже больше всего нравится (или уже нравилась) jet, может быть потому-что она красочная, может быть потому что часто используется и к ней все привыкли. Однако при подготовки ответа на вопрос коллеги, я в очередной раз окунулся в пучину негативных отзывов о радужной шкале и кажется меня эти доводы окончательно убедили в том, что неправильный выбор шкалы может заставить на изображениях видеть то, чего на самом деле нет и не видить того, что на самом деле есть.
В одной из наших предыдущих заметок мы визуализировали данные SST на карте. И в той заметке зашла речь о том, что при нанесении данных на карту, используя метод 'texturemap' в функции geoshow, области с пропущенными значениями (NaN) закрашиваются цветом. Например, используя те же данные SST, построим карту для немного другог района:
Раньше такая конструкция работала, но теперь, по неведомой мне причине MATLAB начал ругаться: "Error using setm (line 41) Expected H to be a handle to a MATLAB graphics object." Просит сделать указатьель (насколько я понимаю, как на русский язык переводится слово Handle, или ещё называют "дескриптор", что быть может даже более прпвильно), не хочет обращаться к карте при помощи gca. Ладно, сделаем через этот самый hadle.
k = worldmap([42 78],[-50 180])
setm(k,'MapProjection','mercator')
Теперь получилось построить подложку карты, осталось нанести данные:
Как мы и предупреждали, значения NaN имеют тёмно синий цвет, что не есть хорошо, так как неотличимы от значений с низкими температурами. Одно из решений мы обсуждали в прошлый раз - заменить 'texturemap' на 'surface':
Часто приходится на карте указывать район нашего исследования, например, если берегов нет или по их очертаниям сразу не определить район. Когда вы показываете маленький район на более крупной и узнаваемой карте, то для зрителя воспринимать ваш доклад, например, будет комфортнее.
Сегодня мы продолжаем цикл заметок посвящённых визуализации в MATLAB и построим на карте прямоугольник, многоугольник,... короче говоря, полигон.
Мы будем идти по стопам прошлой заметки, где мы строили топографическую карту. Поэтому следующий код должен вам оказаться понятен:
Сегодня мы построим топографическую карту (карту рельефа) для всего мира, а также для отдельных районов, используя средства MATLAB. Нам как океанологам, конечно, важна её батиметрическая часть, но строить мы будем и "сухопутную" часть тоже. Сегодня ограничимся только картой, где топография отлита цветом. Изолинии будем рисовать в другой раз.
Начнём с того, что используем родной файл с топографией, который встроен в MATLAB.
load topo
В результате этой команды откроются несколько переменных, главная из которых topo, в которой и содержатся нужные нам значения высот
В прошлый раз (советуем прочитать тот материал, прежде чем читать дальше) для построения карты мы использовали матрицы с широтами и долготам, равные по размерам матрице с данными. Здесь же нет даже векторов, из которых мы эти матрицы можем сделать. Но на самом деле, при равномерной координатной сетке мы можем обойтись и без них. Поэтому сейчас мы построим карту немного другим способом (отличным от того, которым пользовались ранее).
Если мы знаем размер матрицы с топографией и знаем пределы этой карты (по широте и долготе), то можно сказать, что наша координатная сетка нам известна. Функция georefcells по пределам широт и долгот, а также размеру матрицы topo создаст в матлабе референсную ячейку, которую впоследствии мы сможем использовать для построения:
Для того, чтобы цвета выглядели более привычно для географа, лучше воспользоваться функцией demcmap, которая создаст специальную цветовую политру (colormap), подходящую для отображения рельефа.
Если внимательно посмотреть на последние заметки нашего блога, то можно подумать, что он превратился в чисто фотографический. Однако недавно мы хорошенько обдумали перспективы развитя наших сайтов, не побоялись и поставили всё с ног на голову. Поскольку фотографическая тематика постов отныне будет жить только на , этот сайт вновь станет выглядеть так, как он выглядел в самом начале. Да-да, начиналось всё не с фотографических обзоров, а с заметок, типа: Вязание крючком игрушки Бендера (Футурама) Взрослая загадка на смекалку — расшифровать послание инопланетянам Работа с 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()
