Вырезание полигона из данных

Задача — вырезать район из матрицы данных.

В наличии — матрица с данными и две матрицы такого же размера что и матрица с данными, в одной находятся широты в другой долготы.

?нструменты — Matlab и GoogleEarth

Чтобы не отрываться от реальности, скажу что данные — это концентрация льда по SSM/I , NASA Team алгоритм. Решение грязное, медленное и нудное, правильного и красивого мне найти пока не удалось )

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

Сначала нарисуем полигон. Для этого лучше всего подходит Google Earth. Выбираем там «полгигон» и начинаем рисовать, ставя точечки через определённые промежутки. Промежутки не должны быть слишком большими, и вот почему. Попробуйте поставить точки на расстоянии где ни будь градусов 50 по долготе и увидите что линия между ними не прямая а изогнутая, потому что земля, к сожалению, круглая. Маталабовская же функция, которой мы будем вырезать не настолько умная и проведёт между точками просто прямую, что в результате может вылиться в очень большую ошибку. Поэтому точечки ставим друг от друга недалеко, но и не мельчим. Координаты каждой точки записываем.

Следующий этап. Нам нужно найти эти точки в наших данных. То есть такие ячейки матрицы которые находились бы ближе всего к точкам которые мы выбрали. Делается это при помощи нехитрого скрипта

lat = 81;
lon = 15;
di = 0.25;

[xx,yy] = find ((latitude > (lat-di)) & (latitude < (lat+di)) & (longitude > (lon-di)) & (longitude < (lon+di)))
fff = (latitude (xx,yy))
ggg = (longitude (xx,yy))

где

lat — ваша широта (с минутами в десятых)

lon — долгота (с минутами в десятых)

di — окрестность в которой искать точки

latitude — двумерный массив (матрица) широты

longitude — двумерный массив долготы

xx — строчка(и)

yy — столбец(цы)

fff и ggg — показывают что за широты и долготы стоят за полученными столбцами и строчками.

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

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

xall71 = [xall71 xx];
yall71 = [yall71 yy];

save xall71 xall71
save yall71 yall71

Дальше надо собственно вырезать полигон. ?спользовать будем функцию poly2mask. В качестве аргументов мы должны дать ей два вектора с индексами строчек и столбцов, которые у нас уже имеются (xall71 и yall71) а также размеры матрицы из которой мы впоследствии будем вырезать полигон (в нашем случае 448 на 304). Poly2mask создаст маску, по форме нашего полигона, которую мы затем применим к нашим данным и посмотрим на вырезанный кусочек.

%создаём маску как было описано. poly2line должна провести между вашими точками прямые и выбрать всё что ограничено этими прямыми. Этим точкам присваиваются значения 1

mmask = poly2mask(yall71,xall71, 448, 304);

%находим индексы всех точек внутри полигона (у них значение 1)
iindex = find(mmask == 1 );

%вырезаем полигон из данных (в нашем случае данные это map1)

polygone = map1(iindex);

%чтобы посмотреть на наш полигон перемножаем маку и данные, превращая все точки вне полигона в 0

cuted_polygone = map1 .* mmask;

%отрисовываем получившийся полигон

figure
imagesc(cuted_polygone)

Главное в общем то было получить переменную polygone. Проделав то же самое с широтой и долготой вы сможете создать какой ни будь ASCII типа широта\долгота\значение .