Построение точек на карте с береговой линией по заданным координатам при помощи MATLAB.

7 ноября 2018

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

Точки зададим вручную, также, как делали это в прошлый раз (точки находятся в губе Чупа, Белое море).

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]

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

worldmap([66.22 66.37],[33.60 34])
setm(gca,'MapProjection','mercator')
geoshow('landareas.shp','FaceColor',[0.5 0.5 0.5])

Результат, мягко скажем, не очень хороший. Нет даже намёка на существование здесь губы Чупа. Всё дело в том, что файл landareas.shp содержит не очень подробную информацию о береговой черте. Для больших районов это годится, но для небольшых — явно маловато. Давайте оценим масштаб бедствия и посмотрим как наша береговая линия выглядит для всего Белого моря:

worldmap([64 68],[32 45])
setm(gca,'MapProjection','mercator')
geoshow('landareas.shp','FaceColor','none')

Кроме файла landareas.shp, который мы использовали, есть ещё одна береговая линия, встроенная в Матлаб. Загрузить её можно, написав:

load coastlines

В итоге загрузится два вектора с широтами и с долготами береговой линии:

Name             Size            Bytes  Class     

  coastlat      9865x1             78920  double              
  coastlon      9865x1             78920  double

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

Теперь давайте посмотрим, насколько лучше новая береговая черта:

worldmap([64 68],[32 45])
setm(gca,'MapProjection','mercator')
geoshow(coastlat,coastlon,'Color', 'black')

береговая линия для Белого моря в MATLAB

Что ж, кажется, стало только хуже. Ясно, что нам нужна другая береговая линия. Неплохим выбором будет GSHHG (A Global Self-consistent, Hierarchical, High-resolution Geography Database), которую можно загрузить с их сайта.

Скачиваем файл с бинарными данными gshhg-bin-2.3.7.zip и распаковываем его, внутри мы найдём несколько файлов. Нас будет интересовать береговая линия в самом высоком разрешении (файл gshhs_f.b)

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

shorelines = gshhs('gshhs_f.b');

Однако, если вам не нужно строить карту всего Мирового океана, для того, чтобы процесс не был слишком долгим, лучше указать пределы нашей карты при помощи двух векторов:

shorelines = gshhs('gshhs_f.b',[64 68],[32 45]);

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

worldmap([64 68],[32 45])
setm(gca,'MapProjection','mercator')
geoshow([shorelines.Lat], [shorelines.Lon],'Color', 'black')

подробная береговая линия gshhs для Белого моря в MATLAB

Не сложно заметить, насколько эта береговая линия более подробна. Теперь построим интересующий нас район Белого моря (губа Чупа):

worldmap([66.22 66.37],[33.60 34])
setm(gca,'MapProjection','mercator')
geoshow(shorelines,'FaceColor', [0.5 0.5 0.5])

Теперь нужно нанести на карту точки при помощи функции scatterm. Сделаем точки заполненные (‘filled’) чёрным цветом (‘MarkerFaceColor’,[0 0 0]):

worldmap([66.22 66.37],[33.60 34])
setm(gca,'MapProjection','mercator')
geoshow(shorelines,'FaceColor', [0.5 0.5 0.5]) 
scatterm(lat,lon,'filled','MarkerFaceColor',[0 0 0])

Основной недостаток получившейся карты — невозможность понять определить координаты точек, так как присутствует всего одна линия координатной сетки вдоль широты. Поэтому добавим координатные оси и подписи к ним вручную (сделаем это и для широт и для долгот, хотя в нашем случае это было необходимо только для широт)

worldmap([66.22 66.37],[33.60 34])
setm(gca,'MapProjection','mercator')
geoshow(shorelines,'FaceColor', [0.5 0.5 0.5]) 
scatterm(lat,lon,'filled','MarkerFaceColor',[0 0 0])
setm(gca,'mlinelocation',[33.6:0.1:34.0])
setm(gca,'mlabellocation',[33.6:0.1:34.0])
setm(gca,'plinelocation',[66.22:0.04:66.38])
setm(gca,'plabellocation',[66.22:0.04:66.38])

Установка координатной сетки и подписей в geoshow

Чтобы добиться большей схожести с той картой, которую мы строили в Python сделаем линии координатной сетки чёрного цвета:

setm(gca,'GColor',[0 0 0])

чёрная координатная сетка geoshow matlab

Ну что ж, с задачей мы справились, но напоследок добавим ещё одну деталь, которая возможно вам пригодится — добавим подписи с номерами точек. Для этого создадим массив с номерами точек:

points = {'1','2','3','4','5','6','7','8','9','10'};

А затем, при помощи функции textm добавим точки на карту.

textm(lat,lon,points)

Поскольку текст мы нанесли в те же координаты, что и точки, то они наложились друг на друга, поэтому немного увеличим значения широт и долгот. Окончательно код будет выглядеть так:

worldmap([66.22 66.37],[33.60 34])
setm(gca,'MapProjection','mercator')
geoshow(shorelines,'FaceColor', [0.5 0.5 0.5]) 
scatterm(lat,lon,'filled','MarkerFaceColor',[0 0 0])
setm(gca,'mlinelocation',[33.6:0.1:34.0])
setm(gca,'mlabellocation',[33.6:0.1:34.0])
setm(gca,'plinelocation',[66.22:0.04:66.38])
setm(gca,'plabellocation',[66.22:0.04:66.38])
setm(gca,'GColor',[0 0 0])
textm(lat+0.0025,lon+0.0025,points)

P.S. Через несколько месяцев после публикации заметки я вдруг обнаружил, что MATLAB опять продемонстрировал свою неспособность автоматически подбирать логичные настройки. Вы наверное более наблюдательны и сразу увидели, что у нас несколько одинаковых широт «66.3». Причина ясна — MATLAB округляет значения координат до десятых. Исправить ситуацию легко, нужно использовать ещё один параметр:

setm(gca,'PLabelRound',-2)

где «-2» означает, что округлять нужно до второго знака после запятой. Если бы мы поставили «0», то округление происходило бы до целых, «1» — до десятков и т.п.

Все карты в заметке я перестраивать не буду, приведу только один новый пример: