Генерализация полигонов

[KolesovDmitry]

И все-таки, мне кажется, что упрощение можно произвести и в системе координат широта-долгота....

если GRASS будет в этом случае считать расстояние хотя-бы по большой дуге (чего вроде не делается), то можно.

Любопытства ради лезем в исходники (lib/vector/diglib/prune.c), смотрим, как измеряется расстояние:

Код: Выделить всё

* The distance between a point M                ->   ->
 * and a AD segment is given                  || AM ^ AD ||
 * by the following formula :            d = ---------------
 *                                                  ->
 *                                               || AD ||
 *
 *                     ->   ->                             ->
 * When comparing   || AM ^ AD ||   and    t = thresh * || AD ||
 *
 *                     ->   ->       ->   ->
 * we call  sqdist = | AM ^ AD | = | OA ^ OM + beta | 
 *
 *                  ->   ->
 *  with     beta = OA ^ OD 

Как я понимаю, выражение

Код: Выделить всё

   
   ->     ->                            
|| AM ^ AD || 

означает обычное векторное произведение, и, соответственно, делаем вывод, что расстояния измеряются не по дуге.

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

С одной стороны -- да. Но конкретно в случае GRASS, как мне кажется, это не означает автоматически невозможность генерализации площадных объектво в системе широта/долгота:

Во-первых, GRASS -- типично UNIX'овая программа с ее философией "Пользователю лучше знать, что он хочет. Если он применяет методы там, где они неприменимы -- он сам себе злобный Буратино", тем более что можно привести примеры задач, когда анизотропность нисколько не мешает.

Во-вторых, приведенный вами аргумент не мешает GRASS произвести генерализацию линейных объектов в системе координат широта/долгота, проблема только с площадными.

Поэтому эти два соображения и наталкивают меня на мысль, что генерализовать площади должно быть можно и в областях типа широта/долгота, и проблема не в принципиальной невозможности, а, скорее, в неумении пользоваться инструментами (всякими v.generalize, v.edit, v.clean)

[gamm]

Поэтому эти два соображения и наталкивают меня на мысль, что генерализовать площади должно быть можно и в областях типа широта/долгота, и проблема не в принципиальной невозможности, а, скорее, в неумении пользоваться инструментами (всякими v.generalize, v.edit, v.clean)

а меня это наталкивает на мысль, что в GRASS (как и в большинстве программ) можно делать глупости. Но это не означает. что их нужно делать, чтобы потом "дивиться" результатам

[Александр Мурый]

То, что GRASS не упрощает (или ошибочно упрощает) векторные объекты при "широте/долготе", уже должно как бы намекать, что в данном случае мы делаем глупость. Надо посмотреть, как считает v.generalize..

По "большой дуге" считают модули v.distance:

In lat-long locations v.distance gives distances (dist and to_along) in meters not in degrees calculated as geodesic distances on a sphere.

и d.geodesic:

d.geodesic displays a geodesic line in the active frame on the user's graphics monitor. This is also known as the great circle line and traces the shortest distance between two user-specified points on the curved surface of a longitude/latitude data set.

[KolesovDmitry]

а меня это наталкивает на мысль, что в GRASS (как и в большинстве программ) можно делать глупости. Но это не означает. что их нужно делать, чтобы потом "дивиться" результатам

Надо посмотреть, как считает v.generalize.

Если я правильно разобрался, при расчете расстояний v.generalize вызывает функцию dig_distance2_point_to_line, которая описана в lib/vector/diglib/line_dist.c.

А в ней, если отбросить всякие проверки на валидность и т.п., суть проста:

Код: Выделить всё

...
    dx = x2 - x1;
    dy = y2 - y1;
    dz = z2 - z1;
...
    return (dx * dx + dy * dy + dz * dz);

Т.о. видим, что и в v.generalize расчет идет без всяких дуг.

Спойлер

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

Код: Выделить всё

t = (dx * (x - x1) + dy * (y - y1) + dz * (z - z1)) / (dx * dx +  dy * dy +  dz * dz);
...
tpx = dx * t + x1;
tpy = dy * t + y1;
tpz = dz * t + z1;
dx = tpx - x;
dy = tpy - y;
dz = tpz - z;
...

Правда, на главный вывод это не влияет: расчеты производятся не на сфере.

[gamm]

Если я правильно разобрался, при расчете расстояний v.generalize вызывает функцию dig_distance2_point_to_line, которая описана в lib/vector/diglib/line_dist.c.

А в ней, если отбросить всякие проверки на валидность и т.п., суть проста:

Код: Выделить всё

...
    dx = x2 - x1;
    dy = y2 - y1;
    dz = z2 - z1;
...
    return (dx * dx + dy * dy + dz * dz);

Т.о. видим, что и в v.generalize расчет идет без всяких дуг.

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

[KolesovDmitry]

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

А почему косинус? Я правильно понимаю, что он был выбран эвристически, просто как монотонная функция, которая с увеличением широты падает от 1 до 0? Или же он был выбран по другим соображениям?

[gamm]

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

А почему косинус? Я правильно понимаю, что он был выбран эвристически, просто как монотонная функция, которая с увеличением широты падает от 1 до 0? Или же он был выбран по другим соображениям?

Неправильно. Расстояние по горизонтали (по долготе), равное alpha в градуса, в метрах равно pi/180*alpha*R(широта), а R(широта)=R0*cos(широта), если считать Землю шаром радиуса R0.

[Denis Rykov]

Похоже, что проблемы с кодировками GRASS в QGIS решили и я решил попробовать очередную попытку упрощения данных, но споткнулся на первом же шаге - загрузке шейпа в GRASS. Использую команду:

Код: Выделить всё

v.in.ogr -o dsn=C:/regions2000_wgs.shp output=result snap=100 min_area=1

Шейп в Меркаторе. Картина следующая (исходный слой и часть, импортированная в GRASS (красным)): почему-то выгрузились только некоторые фрагменты оригинального слоя.

Спойлер

v.in.ogr -o dsn=C:/personal/next/regions2000_wgs_shp/regions2000/regions2000_wgs.shp output=export snap=100 min_area=1

Datum not recognised by GRASS and no parameters found

Over-riding projection check

Layer: regions2000_wgs

Importing map 1517 features...

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

Area size [0.0e+000], area not imported

1 feature without geometry

-----------------------------------------------------

Building topology for vector map ...

Registering primitives...
1000
1272 primitives registered

33561 vertices registered

Building areas...


1253 areas built

1254 isles built

Attaching islands...


Attaching centroids...


Number of nodes: 1262

Number of primitives: 1272

Number of points: 0

Number of lines: 0

Number of boundaries: 1272

Number of centroids: 0

Number of areas: 1253

Number of isles: 1254

Number of incorrect boundaries: 19

Number of areas without centroid: 1253

-----------------------------------------------------

Cleaning polygons, result is not guaranteed!


Building topology for vector map ...

Number of nodes: 1262

Number of primitives: 1272

Number of points: 0

Number of lines: 0

Number of boundaries: 1272

Number of centroids: 0

Number of areas: -

Number of isles: -

-----------------------------------------------------

Snap boundaries (threshold = 1.000e+002):

-----------------------------------------------------

Break polygons:

-----------------------------------------------------

Remove duplicates:

-----------------------------------------------------

Break boundaries:

-----------------------------------------------------

Remove duplicates:

-----------------------------------------------------

Clean boundaries at nodes:

-----------------------------------------------------

Change dangles to lines:

-----------------------------------------------------

Remove bridges:

-----------------------------------------------------

Building topology for vector map ...

Building areas...


1264 areas built

1258 isles built

Attaching islands...


Number of nodes: 1338

Number of primitives: 1389

Number of points: 0

Number of lines: 0

Number of boundaries: 1389

Number of centroids: 0

Number of areas: 1264

Number of isles: 1258

Number of areas without centroid: 1264

Layer: regions2000_wgs

-----------------------------------------------------

-----------------------------------------------------

Building topology for vector map ...

Registering primitives...
1000 2000
2545 primitives registered

32953 vertices registered

Building areas...


1264 areas built

1258 isles built

Attaching islands...


Attaching centroids...


Number of nodes: 2539

Number of primitives: 2545

Number of points: 0

Number of lines: 0

Number of boundaries: 1281

Number of centroids: 1264

Number of areas: 1264

Number of isles: 1258

-----------------------------------------------------

1520 input polygons

Total area: 1.464541e+017 (1264 areas)

Overlapping area: 0.000000e+000 (0 areas)

Area without category: 0.000000e+000 (0 areas)
Successfully finished

[Александр Мурый]

Шейп-файл импортируется в уже готовую область в проекции Меркатора? А если без опций "snap" и "min_area"?

Можно попробовать создать область на базе этого шейпа:

Код: Выделить всё

v.in.ogr -e -o dsn=C:/regions2000_wgs.shp output=result location=regions2000_wgs 

[Александр Мурый]

Разобрался немного. При импорте в GRASS слоя с тыщами полигонов можно использовать опцию "-c" модуля v.in.ogr ("не зачищать полигоны"). И, чтобы отсеять совсем малые площади, можно исп-ть min_area=???

Дальше зачищаем полигоны уже в созданной из шейпа области:

Код: Выделить всё

v.clean in=regions2010 out=regions2010__clean tool=bpol,rmdupl type=boundary

bpol: break (topologically clean) polygons (imported from non topological format, like ShapeFile). Boundaries are broken on each point shared between 2 and more polygons where angles of segments are different

rmdupl: remove duplicate geometry features (pay attention to categories!)

Исходя из мануала v.clean, он сохраняет топологию полигонов, в отличие от v.generalize (который прежде всего для линий).
Тогда можно сделать генерализацию так (например):

Код: Выделить всё

v.clean in=regions2010__clean out=regions2010__clean_prune tool=prune type=boundary thresh=3000

Получается примерно так:

gener.png (19.32 КБ) 7781 просмотр

Надо поэкспериментировать и с v.generalize..