Как нарисовать прямоугольник на карте в MATLAB

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

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

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

[Z, refvec] = etopo('etopo1_ice_c_f4.flt',1, [40 65],[-65 10]);
worldmap ([40 65],[-65 10])
setm(gca,'MapProjection','mercator')
geoshow(Z, refvec, 'DisplayType', 'texturemap');
demcmap(Z)

Получилась такая карта:

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

Теперь давайте нарисуем прямоугольник. На самом деле, рисование прямоугольника — это частный случай рисования полигона любой формы. Смысл в том, что вы указывает два вектора, в первом векторе содержатся широты всех точек вашего многоугольника, а во втором — их долготы. Само построение производится через geoshow с параметром ‘DisplayType’ установленным на ‘polygon’:

worldmap ([40 65],[-65 10])
setm(gca,'MapProjection','mercator')
geoshow(Z, refvec, 'DisplayType', 'texturemap');
demcmap(Z)
geoshow ([42 62 62 42 42],[-60 -60 -40 -40 -60],...
    'DisplayType','polygon')

Что ж, вряд ли вам понравится вариант, в котором весь выделенный район закрашен. Поэтому сделаем его полупрозрачным через параметр ‘FaceAlpha’. Кроме того, сразу изменим его цвет на нейтральный серый у увеличим толщину линии:

worldmap ([40 65],[-65 10])
setm(gca,'MapProjection','mercator')
geoshow(Z, refvec, 'DisplayType', 'texturemap');
demcmap(Z)
geoshow ([42 62 62 42 42],[-60 -60 -40 -40 -60],...
    'DisplayType','polygon','FaceAlpha',0.5,...
    'FaceColor',[0.5 0.5 0.5],'linewidth',2)

Само-собой, при желании вы можете выставить ‘FaceAlpha’ на 0 и тем самым вообще избавиться от заполнения, оставив только границы района. 

Двигаемся дальше. У описанного выше подхода есть важный недостаток. Мы специально начали с построение в проекции Меркатора, давайте теперь посмотрим, как будет картинка выглядеть в другой проекции (оставим ‘eqdconic’, которую Матлаб устанавливает по умолчанию):

worldmap ([40 65],[-65 10])
geoshow(Z, refvec, 'DisplayType', 'texturemap');
demcmap(Z)
geoshow ([42 62 62 42 42],[-60 -60 -40 -40 -60],...
    'DisplayType','polygon','FaceAlpha',0.5,...
    'FaceColor',[0.5 0.5 0.5],'linewidth',2)

Хорошо видно, что точки на карте соединены прямыми линиями, в то время как параллели изогнуты. Это выглядит не очень симпатично (за исключением конечно, случаев, когда ваш район действительно так выглядит). Для таких случаев MATLAB предлагает «простой» выход — указывать точки вашего полигона почаще. Указывать вручную такое количество точек лениво и, на наше счастье, есть функция которая сделает это за нас:

[lat_a, lon_a] = outlinegeoquad([42 62],[-60 -40],1,1);

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

Ну и ещё пару слов нужно сказать о некоторых нюансах построения полигона. Область, которая на карте закрашена лежит справа от линий, которые вы проводите. Т.е. если вы очерчиваете многоугольник по ходу часовой стрелки, то закрасится внутренность многоугольника, если же против часовой, то закрасится внешняя часть. Это легко увидеть, зеркально отразив наши векторы с широтами и долготам точек:

worldmap ([40 65],[-65 10])
geoshow(Z, refvec, 'DisplayType', 'texturemap');
demcmap(Z)
geoshow (flip(lat_a), flip(lon_a),...
    'DisplayType','polygon','FaceAlpha',0.5,...
    'FaceColor',[0.5 0.5 0.5],'linewidth',2)

Вероятно, такой вариант может показать более приемлемым, так как логично, что мы именно интересующий нас район оставляем незакрашенным. Однако, в этом варианте рамка рисуется с внутренней и с внешней стороны полигона (т.е.если мы хотим, чтобы рамка района была толстой, то рамка карты тоже станет толстой). Я сильно не копал «официальное» решение этой проблемы, возможно оно простое, но первое решение, которое мне пришло в голову показалось тоже весьма несложным. Нужно просто построить два полигона — один с внешним закрашиванием, но с нежирной рамкой, а второй с внутренним закрашиванием, но с параметром ‘FaceAlpha’ равным нулю и толстой рамкой:

worldmap ([40 65],[-65 10])
geoshow(Z, refvec, 'DisplayType', 'texturemap');
demcmap(Z)
geoshow (flip(lat_a), flip(lon_a),...
    'DisplayType','polygon','FaceAlpha',0.5,...
    'FaceColor',[0.5 0.5 0.5])
geoshow (lat_a, lon_a,...
    'DisplayType','polygon','FaceAlpha',0,'linewidth',3)

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