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

3 ноября 2018

Сегодня мы рассмотрим несколько готовых рецептов для визуализации океанологических данных. Эта заметка именно с рецептами и не претендует на полноценное руководство, которое для Матлаба существует и достаточно объёмное (около 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)

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

sst = permute(sst, [2 1]);

(но это только если смущает, нас не смущает и мы не будем этого делать).

Да, на полученной картинке мы видим распределение температуры и очертания континентов, но картой такую визуализацию матрицы пока что совсем нельзя назвать. Что ж, идём дальше. А дальше нам понадобятся данные о широтах и долготах для нашей координатной сетки, поэтому «вынимаем» из файла и их тоже:

lat = ncread('A20182442018273.L3m_MO_SST_sst_9km.nc','lat');
lon = ncread('A20182442018273.L3m_MO_SST_sst_9km.nc','lon');

Теперь посмотрим на все наши переменные (либо взглянув в окно Workspace, либо, набрав в командной строке команду whos):

whos
  Name         Size                 Bytes  Class    

  lat       2160x1                  17280  double              
  lon       4320x1                  34560  double              
  sst       4320x2160            74649600  double

Итого, у нас есть двумерная матрица sst и два одномерных вектора lat и lon. Однако, для построения карт нам нужно будет сделать из lat и lon двумерные матрицы такого же размера, что и sst. Дело в том, что функции, при помощи которых мы будем строить карты, требуют, чтобы все три матрицы были одинакового размера Т.е. для каждой ячейки матрицы sst должны существовать такие же ячейки в матрицах lat и lon, содержащие координаты нужной нам ячейки. Всё это можно сделать лёгким движением руки при помощи функции ndgrid:

[lon,lat]=ndgrid(lon,lat);

Обратите внимание, что, если вы переворачивали матрицу функцией permute (см. выше), то вам нужно будет поменять в коде lat и lon местами, вот так:

[lat,lon]=ndgrid(lat,lon);

Теперь начинаем строить саму карту и прежде всего нам нужна «подложка» с какой-нибудь картографической проекцией. На самом деле MATLAB позволяет делать это разными способами, мы используем пожалуй, самый интуитивно понятный:

worldmap('World')

В качестве картографической проекции для карты мира по умолчанию выбирается проекция Робинсона (‘robinson’), но это, конечно, можно будет поменять.

Что ж, подложка у нас есть, осталось нанести на неё наши данные. Сделаем мы это при помощи функции geoshow, в которой в качестве аргументов укажем наши три двумерные матрицы lat, lon и sst:

geoshow(lat,lon,sst)

Визуализация карты в MATLAB через geoshow с параметром 'image'

Такая карта, состоящая только из белого и чёрного цветов нас вряд ли устроит. Дело в том, что у функции geoshow аргумент ‘displaytype’, отвечающий за тип построения, установлен по умолчанию на ‘image’. Этот тип нам подходит, пожалуй, меньше всего, поэтому мы поменяем его на ‘texturemap’, которые наиболее похож на функцию imagesc. Обратите внимание, что MATLAB поймёт, даже если мы вместо ‘texturemap’ напишем ‘texture’:

geoshow(lat,lon,sst,'displaytype','texture')

Визуализация карты в MATLAB через geoshow с параметром 'texture' или ' texturemap'

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

sst = ncread('A20182442018273.L3m_MO_SST_sst_9km.nc','sst');
lat = ncread('A20182442018273.L3m_MO_SST_sst_9km.nc','lat');
lon = ncread('A20182442018273.L3m_MO_SST_sst_9km.nc','lon');
[lon,lat]=ndgrid(lon,lat);
worldmap('World')
geoshow(lat,lon,sst,'displaytype','texture')

Задача выполнена, однако очень формально и обычно этого не достаточно. Теперь давайте пофантазируем, какие ещё видоизменения полученной карты нам могут пригодится в работе. Прежде всего следует отметить, что океанологи чаще работают с отдельными районами, а не со всем Мировым океаном. Ограничить район достаточно просто, для этого функции worldmap нужно задать в качестве аргументов два вектора с пределами по широте и пределами по долготе.

worldmap([40 80],[-40 20])
geoshow(lat,lon,sst,'displaytype','texture')

Ограничиваем район построения при визуализации данных на карте в MATLAB

Обратите внимание, что проекция поменялась, теперь она равнопромежуточную коническую (‘eqdconic’).

Глядя на такую карту, сложно понять, где очертания континентов и островов, а где пропуски в данных, так как и то и другое обозначается у нас в данных при помощи NaN. Нанести очертания берега можно разными способами и мы ещё рассмотрим их в дальнейших заметках. Сейчас же мы используем встроенный в MATLAB shp-файл с береговой линией, нас он вполне устроит:

worldmap([40 80],[-40 20])
geoshow(lat,lon,sst,'displaytype','texture')
geoshow('landareas.shp')

Построение береговой черты (материки, океаны) в MATLAB

Теперь я сделаю ещё несколько видоизменений на свой вкус, прокомментировав сделанное прямо в коде:

worldmap([40 80],[-40 20])
geoshow(lat,lon,sst,'displaytype','texture') 
%установим границы цветовой шкалы от -2 до 25 градусов:
caxis([-2 25])
%изменим цвет континентов на серый, указав три значения RGB:
geoshow('landareas.shp','FaceColor',[0.5 0.5 0.5])
%изменим цветовую карту на jet:
colormap jet
% добавим сбоку шкалу с цветами:
colorbar

Изменяем цветовую шкалу и диапазон значений на карте MATLAB

Полученную карту можно раскритиковать за то, что невозможно отличить значения с низкой температурой от NaN (пропусков в данных). Самый простой способ избежать этого — выбрать другой тип отображения данных на карте (заменим ‘texturemap’ на ‘surface’):

worldmap([40 80],[-40 20])
geoshow(lat,lon,sst,'displaytype','surface') 
caxis([-2 25])
geoshow('landareas.shp','FaceColor',[0.5 0.5 0.5])
colormap jet
colorbar

Визуализация карты в MATLAB через geoshow с параметром 'surface'

На данных картинках ‘texture’ не отличим от ‘surface’, но это потому что матрицы у нас слишком большие, если бы сетка будет не такой густой, то ‘texture’ выдаст нам изображение, состоящее из прямоугольников, а ‘surface’ создаст сглаженное изображение:

Визуализация карты в MATLAB через geoshow сравнение 'texture' или 'surface'

Всего в geoshow существуюет пять типов отображения: ‘surface’, ‘contour’, ‘mesh’, ‘texturemap’, ‘image’. С тремя из них, мы уже познакомились. Что касается двух оставшихся, то, прежде чем их попробовать, давайте условимся, что для построения мы будем брать не все данные, а только каждое десятое значение. Сделаем мы это из-за того, что наши матрицы достаточно большого размера, поэтому в случае с ‘mesh’ мы без масштабирования просто не разглядим результат, а в случае с ‘contour’ процесс построения будет долгим.

Что ж, давайте попробуем. Параметр ‘mesh’ представит наши данные на карте в виде сетки:

worldmap([40 80],[-40 20])
geoshow(lat(1:10:end,1:10:end),lon(1:10:end,1:10:end),...
    sst(1:10:end,1:10:end),'displaytype','mesh') 
caxis([-2 25])
geoshow('landareas.shp','FaceColor',[0.5 0.5 0.5])
colormap jet
colorbar

Визуализация карты в MATLAB через geoshow с параметром 'mesh'

Для построения с параметром ‘contour’ сразу обозначим шаг прорисовки изолиний (‘LevelStep’,1) иначе по умолчанию MATLAB построит слишком мало изолиний:

worldmap([40 80],[-40 20])
geoshow(lat(1:10:end, 1:10:end),lon(1:10:end, 1:10:end),...
    sst(1:10:end, 1:10:end),'displaytype','contour','LevelStep',1)
caxis([-2 25])
geoshow('landareas.shp','FaceColor',[0.5 0.5 0.5])
colormap jet
colorbar

Визуализация карты в MATLAB через geoshow с параметром 'contour'

По умолчанию изолинии имеют цвет, это слегка непривычно. Сделать их чёрными можно, присвоив параметру ‘LineColor’ значение ‘black’. Дополнительно заполним цветом промежутки между изолиниями (‘Fill’,’on’):

worldmap([40 80],[-40 20])
geoshow(lat(1:10:end,1:10:end),lon(1:10:end,1:10:end),...
   sst(1:10:end,1:10:end),'displaytype','contour','Fill','on',...
    'LineColor','black','LevelStep',1)
caxis([-2 25])
geoshow('landareas.shp','FaceColor',[0.5 0.5 0.5])
colormap jet
colorbar

имитируем conourf на карте в MATLAB, а именно contiurfm в функции geoshow, получаем закрашенные изолинии

Ну и пожалуй последнее на сегодня видоизменение, которое внесём в нашу карту — поменяем проекцию, например, на проекцию Меркатора. Конечно, выбор для океанолога не случаен, это проекция является одной из самых популярных, так как используется в навигационных картах и хороша для отображения данных, связанных с направлением (направления ветра, направления течений). Тонкая настройка карты выполняется при помощи функции setm. Необходимо помнить, что задавать проекцию следует в самом начале, сразу после worldmap:

worldmap([40 80],[-40 20])
setm(gca,'MapProjection','mercator')
geoshow(lat,lon,sst,'displaytype','surface') 
caxis([-2 25])
geoshow('landareas.shp','FaceColor',[0.5 0.5 0.5])
colormap jet
colorbar

изменение картографической проекции в geoshow MATLAB, меркатор

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