Страница 1 из 3

Определение линейных размеров полигонов в ArcGis

Добавлено: 08 окт 2007, 11:43
Борис
Добрый день!

Есть полигональная шейп-тема.
Помогите, пожалуйста, определить размер полигонов - то есть минимальный (ширину) и максимальный (длину) размеры. Дело осложняется тем, что полигоны - не прямоугольной формы, а всякие изрезанные и пр.

Добавлено: 14 ноя 2007, 17:07
AndreyL
Если тема еще интересует, то можно попробовать вписать каждый полигон в прямоугольник минимальной площади, повернутый на определенный угол. Угол найти минимизируя площадь прямоугольника.

Добавлено: 17 ноя 2007, 23:55
Борис
Да, тема еще интересует (очень:).
AndreyL, поясните, пожалуйста, подробнее, как это сделать.
Лучше всего - в ArcGis, но есть и ArcView, и MapInfo.

Добавлено: 18 ноя 2007, 16:40
AndreyL
Есть подозрение, что придется задействовать математические пакеты (или хотя бы Excel) – решение, похоже, численное. Идея то проста: поворот группы точек (всех вертексов полигона) вокруг центра координат на заданный угол сложности не представляет (для пущей страсти можно сначала параллельно перенести ближе к центру). Далее посчитать минимумы и максимумы по Х и по Y и площадь (или периметр, если хотите) тоже не сложно. Теперь осталось подобрать нужный угол поворота, а вот это задачка численная.

Добавлено: 20 ноя 2007, 12:07
Борис
AndreyL, а не могли бы Вы бодробнее (конкретнее) описать последовательность действий. С общим ходом рассуждений согласен, но как все это вополтить... :roll: ???

Добавлено: 20 ноя 2007, 13:29
AndreyL
Вот так примерно можно на Вольфрамовской Математике (точки, естественно, для примера):
points = Table[{Random[], Random[]*1.5}, {i, 25}];

ClearAll[x]
mat = {{Cos[x], Sin[x]}, {-Sin[x], Cos[x]}};

V = (Max[Transpose[points.mat][[1]]] -
Min[Transpose[points.mat][[1]]])*(Max[
Transpose[points.mat][[2]]] - Min[Transpose[points.mat][[2]]]) //
Simplify;
rez = FindMinimum[V, {x, 0, Pi}]

rect = {{Min[Transpose[points.mat][[1]]],
Min[Transpose[points.mat][[2]]]}, {Min[Transpose[points.mat][[1]]],
Max[Transpose[points.mat][[2]]]}, {Max[Transpose[points.mat][[1]]],
Max[Transpose[points.mat][[2]]]}, {Max[Transpose[points.mat][[1]]],
Min[Transpose[points.mat][[2]]]}, {Min[Transpose[points.mat][[1]]],
Min[Transpose[points.mat][[2]]]}} /. rez[[2]];
mat2 = {{Cos[-x], Sin[-x]}, {-Sin[-x], Cos[-x]}};
rrect = rect.mat2 /. rez[[2]];

gr1 = ListPlot[points, PlotStyle -> PointSize[0.02],
DisplayFunction -> Identity];
gr2 = ListPlot[rrect, PlotJoined -> True, DisplayFunction -> Identity];
Show[gr1, gr2, AspectRatio -> Automatic, DisplayFunction -> $DisplayFunction];

Добавлено: 21 ноя 2007, 22:19
Борис
Спасибо за ответ :) , разбираюсь :? , а не могли бы Вы пояснить на конкретном примере. Последовательность действий такая: 1)есть полигон; 2)узлы полигона превращены в точечный шейп; 3)определены координаты всех узлов; 4)получены,например, следующие аттрибуты:

polygon vertice _____POINT_X__________ POINT_Y
1_____ 1_____ 545649,517515_____ 6023352,726300
1_____ 2_____ 721286,473581_____ 5968844,705450
1_____ 3_____ 642552,665689_____ 5955520,522580
1_____ 4_____ 577143,040672_____ 5980957,598980
1_____ 5_____ 522635,019824_____ 6001549,517960
1_____ 6_____ 545649,517515_____ 6023352,726300

как быть дальше? как на автомате прочесть их Математикой (ведь это только один полигон (колонка "polygon", а на самом деле их много). Принцип построения атрибут.таблицы простой: колонка номера полигона, далее номер узла в полигоне, далее координаты.

Добавлено: 22 ноя 2007, 13:18
AndreyL
Конвертора DBase-таблиц в Математику я не знаю, но можно конвертировать таблицу в текстовый файл – такое Математика нормально читает. А потом в самой Математике разбить таблицу по признаку polygon:

dat = Import[FileName, "Table"];
list = Transpose[dat][[1]];
UnVal = {};
n = Length
  • ;
    Do[If[Position[UnVal, list[]] == {}, UnVal = Insert[UnVal, list[], -1]], {i, n}]

    PointsTable = Table[{}, {i, Length[UnVal]}];
    Do[k = Position[UnVal, list[]][[1, 1]];
    PointsTable[[k]] = Insert[PointsTable[[k]], {dat[[i, 3]], dat[[i, 4]]}, -1], {i, n}]
    PointsTable


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

Добавлено: 22 ноя 2007, 18:53
geologic
Подобрать угол поворота наверно будет сложно, первое что на ум приходит - тупо поворачивать с заданной дискретностью и проверять, макс-не макс. Если же угол подобран, то просчитать остальное уже проще. Задачка интересная, кстати, со многими другими связана.

Добавлено: 23 ноя 2007, 09:25
AndreyL
geologic писал(а):Подобрать угол поворота наверно будет сложно, первое что на ум приходит - тупо поворачивать с заданной дискретностью и проверять, макс-не макс.
Да упаси меня Господи от таких подвигов. Тогда уж проще выбросить ГИС вместе с компьютером, и решать задачу на бумажной карте, вооружившись рейсшиной и калькулятором. Мы же четыре поста назад соорудили код автоматического поиска угла. Только ГИС, как и любая информационная система, на серьезные подвиги неспособна, поэтому задействуется математическое программное обеспечение, хотя алгоритм оптимизации можно и на VBA написать (если задача не единичная и делать кучу конвертов ленивей, чем писать код).

Добавлено: 23 ноя 2007, 12:26
geologic
Осенило :) Определить длинную ось можно, найдя максимальное расстояние в группе точек. Это можно сделать простым запросом SQL в реляционной БД, к таблице с колонками X, Y, IDA (индекс "бывшего" полигона, признак группы). Группа, кстати, не обязательно должна быть бывшим полигоном, любое облако годится.

Определить азимут длинной оси можно, но механизм сложнее. Вычислить вторую ось как перпендикуляр... Здесь мысль тормозит пока :(

Добавлено: 23 ноя 2007, 14:43
AndreyL
geologic писал(а):Осенило :) Определить длинную ось можно, найдя максимальное расстояние в группе точек.
В принципе Вы правы, если нужно найти прямоугольник с максимальной длиной. Можно поставить задачу найти прямоугольник с минимальной шириной. Можно искать прямоугольник с минимальной длиной или максимальной шириной (типа минимакса и максимина). Мы пока ищем прямоугольник с минимальной площадью. Это вопрос – что мы должны понимать под длиной и шириной. И кроме Бориса на него никто не ответит.

Борис! Я тут пока накидал нечто, проводящее расчеты после загрузки и разбора таблицы, надеюсь, разберетесь и подправите по Вашим требованиям:
mm = Length[PointsTable]
vertTab = Table[{}, {i, mm}];
rectTab = Table[{}, {i, mm}];

Do[
points = PointsTable[[mmi]];
rot = Transpose[points.mat];
V = (Max[rot[[1]]] - Min[rot[[1]]])*(Max[rot[[2]]] - Min[rot[[2]]]) //
Simplify;
rez = FindMinimum[V, {x, 0, Pi}];

rot = Transpose[points.mat] /. rez[[2]];
rect = {{Min[rot[[1]]], Min[rot[[2]]]}, {Min[rot[[1]]],
Max[rot[[2]]]}, {Max[rot[[1]]], Max[rot[[2]]]}, {Max[rot[[1]]],
Min[rot[[2]]]}};

rrect = rect.mat2 /. rez[[2]];

perim = 2*(Max[rot[[1]]] - Min[rot[[1]]] + Max[rot[[2]]] - Min[rot[[2]]]);
midPoint = {Sum[rrect[[i, 1]], {i, Length[rrect]}]/Length[rrect],
Sum[rrect[[i, 2]], {i, Length[rrect]}]/Length[rrect]};
angle =
If[(Max[rot[[2]]] - Min[rot[[2]]]) > (Max[rot[[1]]] - Min[rot[[1]]]),
x/Pi*180, x/Pi*180 + 90] /. rez[[2]];
angle = If[angle < 0, angle + 360, angle];
vertTab[[mmi]] =
Table[{UnVal[[mmi]], rrect[[i, 1]], rrect[[i, 2]]}, {i, Length[rrect]}];
rectTab[[mmi]] = {UnVal[[mmi]], midPoint[[1]], midPoint[[2]], rez[[1]],
perim, angle},
{mmi, mm}]


Кстати, этот код, похоже, универсальный в отношении того, что мы понимаем под длиной и шириной – все будет зависеть от целевой функции V. Например, для того решения, что предложил geologic целевая функция будет:
V = -(Max[rot[[1]]] - Min[rot[[1]]]) // Simplify;

Добавлено: 23 ноя 2007, 15:02
geologic
Ну, выборка одно дело, расчёть - другое. Задача ставилась рассчитать для всех длину и ширину. С длиной разобрались, что это просто - у вас свой подход сложился, у меня свой. А вот насчет ширины...

Я бы прежде поставил вопрос "что такое ширина", соответственно. В том смысле, что самому интересно услышать. задача и правда часто стоит, и не только для замкнутых контуров, измерить ширину, мол.

Добавлено: 24 ноя 2007, 23:53
Борис
Всем спасибо, погружаюсь в мир вольфрама 8)
а про ширину --- да, ширина - это (в бытовом плане ( :? ) ) то, что перпенидикулярно длине. Задачка, думаю, и впрямь может быть довольно широкоупотребимая, так что пытаюсь сам и призываю всех желающих создать скрипт (VBA или с коннектом на ту же Математику)
:wink:

Добавлено: 26 ноя 2007, 14:06
AndreyL
Борис писал(а):пытаюсь сам и призываю всех желающих создать скрипт (VBA или с коннектом на ту же Математику)
Коннект на Математику не нужен, если уж писать скрипт, то писать. Математика используется, когда задача единичная – решил и забыл. Или для исследования задачи, чтобы потом перенести это в тот же скрипт.
Если сделаете оболочку, то готов помочь с поиском оптимального значения, поскольку задачу я исследовал пока отвечал Вам.