Количественные характеристики качества ЦМР

[eucrariano]

Вот имеем ЦМР, создали ее из отвекторизированных горизонталей. Путем подбора параметров в v.surf.rst получили более-менее приемлимый результат. В даный момент оцениваю его визуально, то-есть, чтобы не было ступенек и прочих артефактов на поверхности. Есть ли способ получить какие-то количественные характеристики "качества" модели, чтобы можно было сравнить цифры и сказать: вот, например, модель 3, лучше модели 1 и 2, поскольку ... и пошли показатели? Или можно было сравнить модели из грасса и серфера или других источников. Вот здесь http://grass.wz.cz/ch01s02.html#v.surf.rst рассчитывают девиации и т.п., но я что-то не очень понимаю, как это все происходит.
Спасибо.

Из ответа amuriy: "Вот на GRASS-Wiki кое-что по теме: http://grass.osgeo.org/wiki/RST_Spline_Surfaces и Contour_lines_to_DEM"

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

Модуль v.surf.rst имеет опцию -с -- "выполнить процедуру кросс-валидации без аппроксимации растра".
Как я понимаю, кросс-валидация -- это один из главных методов оценки качества ЦМР при интерполяции.

Вот на Google Books версия статьи Hofierka J., Cebecauer T., Suri M., 2007. Optimisation of Interpolation Parameters Using Cross-validation, на которую ссылается модуль.

Из методов оценки ЦМР нашел в сети в статьях еще "Independent check points" и "Surveyed points", то бишь "независимые проверочные точки" и "точки наблюдений на местности".

[eucrariano]

Та же книжица, в более удобном для прочтения виде http://www.springerlink.com/content/978 ... 9&locus=50 , даже превью пропущенных страниц есть, правда в препоганом нечитабельном качестве

[gamm]

Модуль v.surf.rst имеет опцию -с -- "выполнить процедуру кросс-валидации без аппроксимации растра".
Как я понимаю, кросс-валидация -- это один из главных методов оценки качества ЦМР при интерполяции.

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

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

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

Да, хорошо, что вспомнили про Вуда, вот ссылка на его PhD (надо будет таки почитать труды зарубежного коллеги, а то нехорошо). Кстати, в его диссере активно используется GRASS (потом он ещё написал модуль r.param.scale, и в его Landserf'е много чего с рельефом можно делать.

Только вот не пойму, как в GRASS быть с топологическим "каркасом".. Тальвеги (водотоки) выделить можно, водораздельные линии (слегка потанцевав с бубном) тоже можно. Только вот как применить всё это при интерполяции..?

[eucrariano]

Вот статья в тему: http://skagit.meas.ncsu.edu/~helena/gms ... _et_al.pdf
И вообще-то, здесь много чего еще есть: http://skagit.meas.ncsu.edu/~helena/gmslab/papers/ и просто http://skagit.meas.ncsu.edu/~helena/ - все в открытом доступе

[eucrariano]

Ну вот, умельцем был, умельцем остался. Запустил процедуру кроссвалидации. 39400 сек, почти 11 часов считало. Построило векторную карту crossvalid, в таблице атрибутов, открываемой грассом - 10000 значений саt, если dbf файл открыть Кальком, то показывает все 171354 записи. Исходный регион имеет размеры 1725 на 1758. "The program writes the values of parameters used in computation into the comment part of history file elev as well as the following values which help to evaluate the results and choose the suitable parameters: minimum and maximum z values in the data file (zmin_data, zmax_data) and in the interpolated raster map (zmin_int, zmax_int), rescaling parameter used for normalization (dnorm), which influences the tension."
Типичный ламерский вопрос: а где ж-то искомый файл elev. Не могу найти
Или что-то не так делаю? Скрипты из http://www.grassbook.org/examples_menu3rd.php также запускал:
Computing tension/smoothing t130s090...
ERROR: Vector map <precip3d@> not found
ERROR: Vector map <data_cv_t130s090> not found
и так много десятков/сотень раз. И правда, ведь откуда карте precip3d или data_cv_t130s090 взятся, ведь при кроссвалидации только одна crossvalid создается. Короче, начитался я всего, а толку мало, еще и англицкий хромает немного. Что скажете по поводу?

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

Будем разбираться поочерёдно

Запустил процедуру кроссвалидации. 39400 сек, почти 11 часов считало ... Исходный регион имеет размеры 1725 на 1758

Какое разрешение выставлено в свойствах региона перед интерполяцией?

а где ж-то искомый файл elev

Как я понимаю, кросс-валидация -- очень умная штука для оценки качества будущей модели, собственно модель (тот самый elev) она не создаёт. Без этой опции должна получаться ЦМР (elev=..).

Скрипты из http://www.grassbook.org/examples_menu3rd.php также запускал:
Computing tension/smoothing t130s090...
ERROR: Vector map <precip3d@> not found
ERROR: Vector map <data_cv_t130s090> not found
и так много десятков/сотень раз.

Не мучайте себя и процессор Скрипты со страницы http://www.grassbook.org/examples_menu3rd.php (написанные как примеры к книге, которую, кстати, очень рекомендуется почитать любому грасс-юзеру) работают только на "учебном" наборе данных (точнее, наборах, т.к. их несколько). Основной (самый свежий) набор North Carolina (NC, USA) data set. По сути это заменитель старого-доброго набора Spearfish, который активно исп-ся в мануалах. ОЧЕНЬ рекомендую экспериментировать не на живых людях рабочих данных, а на тренировочных наборах.

Скрипт http://www.grassbook.org/scripts3rd/gra ... ion_rst.sh для проверки кросс-валидации, откуда как раз вектор precip3d (тот, что "not found"), из дополнит. набора данных http://skagit.meas.ncsu.edu/~helena/gra ... kia.tar.gz. Если хочется исп-ть скрипт для своих данных, в нём надо заменить значения переменных:
INMAP="наша векторная карта (слой)"
ZCOL="поле с высотами для интерполяции (из атрибутивной таблицы нашей векторной карты)"
OUTFILESTATS="выходной файл с вычисленной статистикой". Вот тут и главная трудность: понять, что там написано

[eucrariano]

Хорошо. По первому пункту: g.region -3
projection: 0 (x,y)
zone: 0
north: 1719
south: 0
west: 0
east: 1755
top: 1.00000000
bottom: 0.00000000
nsres: 0.99652174
nsres3: 1
ewres: 0.99829352
ewres3: 1
tbres: 1
rows: 1725
rows3: 1719
cols: 1758
cols3: 1755
depths: 1
cells: 3032550
3dcells: 3016845

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

Судя по g.region -3p, нужен объём (в смысле 3D-растр)? Или интерполяция 2D?
Потому как

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

top: 1.00000000
bottom: 0.00000000

-- это "мощность" 3D-растра в 1м (ну или что там за условн. единицы..) -- как-то странно..

А можно характеристики вектора? "v.info вектор_для_интерполяции"
Перед интерполяцией делалось "g.region vect=вектор_для_интерполяции res=.."?

[eucrariano]

v.info:
Title: |
| Map scale: 1:1 |
| Map format: native |
| Name of creator: eucrariano |
| Organization: |
| Source date: Sat Jun 26 20:26:35 2010 |
|----------------------------------------------------------------------------|
| Type of Map: vector (level: 2) |
| |
| Number of points: 26 Number of areas: 0 |
| Number of lines: 163 Number of islands: 0 |
| Number of boundaries: 0 Number of faces: 0 |
| Number of centroids: 0 Number of kernels: 0 |
| |
| Map is 3D: No |
| Number of dblinks: 1 |
| |
| Projection: Unknown projection |
| N: 1722.35269916 S: -3.92105348 |
| E: 1761.03663001 W: -5.16043004 |
| |
| Digitization threshold: 0 |
| Comments: |
| |
+----------------------------------------------------------------------------+
WIND:
proj: 0
zone: 0
north: 1719
south: 0
east: 1755
west: 0
cols: 1758
rows: 1725
e-w resol: 0.99829352
n-s resol: 0.99652174
top: 1
bottom: 0
cols3: 1755
rows3: 1719
depths: 1
e-w resol3: 1
n-s resol3: 1
t-b resol: 1

Изучаемая область была создана Sat Jun 26 20:26:35 2010. Словами Хазанова: "А тогда я был еще глупее". Характеристики при создании области выставил по умолчанию, а некоторые "по наитию" изменил. По видимому это и были top: 1.00000000 bottom: 0.00000000. Вот так и еду на этих данных. Перед интерполяцией делались эсперименты типа "g.region vect=вектор_х res=..". В результате была получена слишком зернистая модель (видимо разрешение было слишком большое. И в конце перед интерполяцией сделал "g.region vect=вектор_для_интерполяции res= по видимому 1". В истории bash не сохранилось, я там много чего еще набивал, хранятся последние 20 команд, так ведь? Ну вот и все. Ну и на будущее вопрос. Масштаб 10000, горизонтали через 5 м, полугоризонтали через 2,5 м. Какое разрешение вертикальное создавать? Я, как видно из вышеперечисленного, его установил в 1. Потом при 3д визуализации, я "поднимаю" ползунком рельеф, а то он слишком плоск. Это, наверное, не так правильно, как может показаться на первый взляд . Вот и получается чудо-юдо рыба-кит. Спасибо.

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

Параметры top и bottom нужны только при создании 3D-растров (т.е. volume). В случае обычной ЦМР (как здесь), т.к. она 2D, надо задать только горизонтальное разрешение. В NVIZ рельеф показывается в т.н. "2.5D", но никак не в "настоящем 3D". А увеличивать вертикальный масштаб в NVIZ -- нормально (при горном рельефе у меня стоит не больше 1.0)

Максимальное количество команд в "истории команд" баша можно задать через переменную окружения HISTSIZE (у меня сейчас, например, стоит 10000), размер файла команд HISTFILESIZE. Это же относится и к переменным внутри GRASS-окружения.

Масштаб 10000, горизонтали через 5 м, полугоризонтали через 2,5 м

Мне кажется, гориз.разрешение 1-2 м -- нормально.

Видно, что вектор не 3D, т.е. интерполировать можно только с исп-ем атрибутов (высот). Хотя, может, интерп-я пойдёт быстрее, если задать высоты изолиниям и точкам согласно атрибут.таблице (v.to.3d) и делать v.surf.rst -z ...

[eucrariano]

Пробовал v.to.3d и v.surf.rst -z . . Вроде-бы все по прежнему, на первый взгляд скорость не изменилась. Поскольку проблем хватает еще с 2 и 2,5д, то с 3д будем разбираться позднее. Пока, что нашел вот 2-ю и 3-ю редакцию книги (в свободном доступе), которую

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

Будем изучать. Хотелось бы еще прочитать Robert J. Peckham Gyozo Jordan Digital Terrain Modelling по теме.

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

По поводу качества ЦМР: т.к. GRASS можно связать с R (spgrass6), то логично было бы считать всякую статистику именно в R. Вот нашёл статейку в книге, упомянутой выше ("Digital Terrain Modelling") GRASS And R -- Advanced GIS and Statistical Tools for DEM Analysis (Radoslav Bonk)