Точки экстремума из srtm в GRASS?

[gimran]

Здравствуйте!
Как получить изолинии из srtm в GRASS понятно, а вот как получить точки вершин и точки низин из srtm?
Подскажите, пожалуйста, где в мануале искать.

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

Спец.модуля для подобных изысканий в составе GRASS для нет.

Но: есть старый скрипт <r.localmax>, написанный когда-то Jachym Cepicky на Perl, который находит локальный максимум для растра (т.е., по идее, вершины).

Вот сам модуль:

r.localmax.zip

(622 байт) 400 скачиваний

Работает так:

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

r.localmax -i dem  -o dem.localmax -m 5

, где после "-i" идёт имя ЦМР, после "-o" — имя выходного растра; после "-m" — размер матрицы в ячейках (только нечётные числа).
На выходе команда для <r.mapcalc>, которую ему надо подсунуть (чем больше размер матрицы, тем больше размер команды и тем реже в итоге будут точки).

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

Пример получения вершин для набора данных Spearfish:

1)

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

r.localmax -i elevation.10m  -o localmax -m 9
localmax = if(elevation.10m == max(elevation.10m[-4,-4],elevation.10m[-4,-3],elevation.10m[-4,-2],elevation.10m[-4,-1],elevation.10m[-4,0],elevation.10m[-4,1],elevation.10m[-4,2],elevation.10m[-4,3],elevation.10m[-4,4],elevation.10m[-3,-4],elevation.10m[-3,-3],elevation.10m[-3,-2],elevation.10m[-3,-1],elevation.10m[-3,0],elevation.10m[-3,1],elevation.10m[-3,2],elevation.10m[-3,3],elevation.10m[-3,4],elevation.10m[-2,-4],elevation.10m[-2,-3],elevation.10m[-2,-2],elevation.10m[-2,-1],elevation.10m[-2,0],elevation.10m[-2,1],elevation.10m[-2,2],elevation.10m[-2,3],elevation.10m[-2,4],elevation.10m[-1,-4],elevation.10m[-1,-3],elevation.10m[-1,-2],elevation.10m[-1,-1],elevation.10m[-1,0],elevation.10m[-1,1],elevation.10m[-1,2],elevation.10m[-1,3],elevation.10m[-1,4],elevation.10m[0,-4],elevation.10m[0,-3],elevation.10m[0,-2],elevation.10m[0,-1],elevation.10m[0,0],elevation.10m[0,1],elevation.10m[0,2],elevation.10m[0,3],elevation.10m[0,4],elevation.10m[1,-4],elevation.10m[1,-3],elevation.10m[1,-2],elevation.10m[1,-1],elevation.10m[1,0],elevation.10m[1,1],elevation.10m[1,2],elevation.10m[1,3],elevation.10m[1,4],elevation.10m[2,-4],elevation.10m[2,-3],elevation.10m[2,-2],elevation.10m[2,-1],elevation.10m[2,0],elevation.10m[2,1],elevation.10m[2,2],elevation.10m[2,3],elevation.10m[2,4],elevation.10m[3,-4],elevation.10m[3,-3],elevation.10m[3,-2],elevation.10m[3,-1],elevation.10m[3,0],elevation.10m[3,1],elevation.10m[3,2],elevation.10m[3,3],elevation.10m[3,4],elevation.10m[4,-4],elevation.10m[4,-3],elevation.10m[4,-2],elevation.10m[4,-1],elevation.10m[4,0],elevation.10m[4,1],elevation.10m[4,2],elevation.10m[4,3],elevation.10m[4,4]),1,null())

2) Подсовываем команду <r.mapcalc>:

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

r.mapcalc 'localmax = if(elevation.10m == max(elevation.10m[-4,-4],elevation.10m[-4,-3],elevation.10m[-4,-2],elevation.10m[-4,-1],elevation.10m[-4,0],elevation.10m[-4,1],elevation.10m[-4,2],elevation.10m[-4,3],elevation.10m[-4,4],elevation.10m[-3,-4],elevation.10m[-3,-3],elevation.10m[-3,-2],elevation.10m[-3,-1],elevation.10m[-3,0],elevation.10m[-3,1],elevation.10m[-3,2],elevation.10m[-3,3],elevation.10m[-3,4],elevation.10m[-2,-4],elevation.10m[-2,-3],elevation.10m[-2,-2],elevation.10m[-2,-1],elevation.10m[-2,0],elevation.10m[-2,1],elevation.10m[-2,2],elevation.10m[-2,3],elevation.10m[-2,4],elevation.10m[-1,-4],elevation.10m[-1,-3],elevation.10m[-1,-2],elevation.10m[-1,-1],elevation.10m[-1,0],elevation.10m[-1,1],elevation.10m[-1,2],elevation.10m[-1,3],elevation.10m[-1,4],elevation.10m[0,-4],elevation.10m[0,-3],elevation.10m[0,-2],elevation.10m[0,-1],elevation.10m[0,0],elevation.10m[0,1],elevation.10m[0,2],elevation.10m[0,3],elevation.10m[0,4],elevation.10m[1,-4],elevation.10m[1,-3],elevation.10m[1,-2],elevation.10m[1,-1],elevation.10m[1,0],elevation.10m[1,1],elevation.10m[1,2],elevation.10m[1,3],elevation.10m[1,4],elevation.10m[2,-4],elevation.10m[2,-3],elevation.10m[2,-2],elevation.10m[2,-1],elevation.10m[2,0],elevation.10m[2,1],elevation.10m[2,2],elevation.10m[2,3],elevation.10m[2,4],elevation.10m[3,-4],elevation.10m[3,-3],elevation.10m[3,-2],elevation.10m[3,-1],elevation.10m[3,0],elevation.10m[3,1],elevation.10m[3,2],elevation.10m[3,3],elevation.10m[3,4],elevation.10m[4,-4],elevation.10m[4,-3],elevation.10m[4,-2],elevation.10m[4,-1],elevation.10m[4,0],elevation.10m[4,1],elevation.10m[4,2],elevation.10m[4,3],elevation.10m[4,4]),1,null())'

3) Векторизуем точечный растр:

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

GRASS 6.4.3svn (spearfish60):~/bin > r.to.vect in=localmax out=localmax feat=point

4) Если надо, получаем высоты для точек.
Сначала добавляем новое поле "z" в таблицу атрибутов:

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

v.db.addcol localmax col="z double"

Добавляем в атрибуты высоты с ЦМР:

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

v.what.rast vector=localmax raster=elevation.10m column=z

Итог:

spearfish_localmax.png (279.43 КБ) 5518 просмотров

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

Насчёт локальных минимумов ЦМР: если заменить в выложенном выше скрипте все "max" на "min", то можно вытащить некое подобие минимумов

Пример:

spearfish_localmin.png (288.98 КБ) 5513 просмотров

На картинках видны недостатки этого способа: на расчленённых участках точек много, на ровных — мало. Сделать сеть точек более равномерной можно, видимо, использованием минимального размера матрицы (3) и последующей фильтрацией точек.

-------------------------------------------------------------------------------
P.S. в одну строку пункты №1 и №2 из поста выше выглядят так (в Linux+bash):

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

r.mapcalc "$(r.localmax -i elevation.10m  -o localmin -m 9)"

[gimran]

Спасибо! Попытаюсь установить и изучить.

[gamm]

Спец.модуля для подобных изысканий в составе GRASS для нет.

а локальное направление потока по D8 там есть? ну тогда "низины" (pit) - это точки, в которых потока нет. А вершины - это низины на минус-рельефе (умножить его на -1)