GFC treecover, чтобы убрать лес с SRTM

[darsvid]

Задача:
Определить высоту леса на разных участках территории с неоднородным рельефом; затем — отнять эти значения от данных SRTM 1 arcsec для получения bare earth в первом приближении.

Территория:
Площадь территориии порядка 8 000 км кв. Рельеф неоднородный - есть равнинные участки, речная долина, сильнорасчлененный эрозионный рельеф.

Исходные данные:
SRTM 1 arcsec; Global Forest Change treecover 2000; топокарты 100 000. Задача состоит в том, чтобы оптимально использовать возможности данных SRTM и GFC.

Попробовали 2 варианта:

1) Устанавливать высоту леса в зависимости от % сомкнутости, которые дают данные GFC. Сначала сделали растр с градациями высот, а потом сгладили фильтром, чтобы избежать резких перепадов. Результат получился неплохой, но поскольку связь между % сомкнутости и высотой леса не всегда однозначная, некоторые участки обрабатываются некорректно (либо слишком большое значение - ямы, либо маленькое - остается лес).

2) Прошлись плавающим кругом в 27 ячеек (ориентировочное значение) и посчитали минимумы. Получилась поверхность минимальных значений, которая дает обобщенное представление о рельефе. Тогда SRTM - минимумы = высота растительности. С помощью данных GFC отмаскировали лес, убрали нереалистичные значения (например, по топокартам нет леса высотой более 25 м), опять же сгладили фильтром. Полученный таким образом слой высот вычли из данных SRTM. Такой алгоритм хорошо работает для равнинных участков, а вот на сильнорасчлененном рельефе дает сильно завышенные значения.

3)* Пока писала - нашла в одном из старых обсуждений: можно построить полигоны леса, и буферы к ним - один внутри, другой снаружи. Посчитав среднюю высоту по каждому контуру, можно получить высоту леса, как разность этих значений. Но опять же встает вопрос - как "размазать" это значение по всему контуру.

Еще нашли вот такую публикацию, но кода нет и поэтому непонятно, что и как они там считали.

Возможно вы уже делали что-то подобное - поделитесь опытом, будем признательны за любые идеи/ подсказки.

[bim2010]

Зная значения высот по ЦМР и ЦММ в одной точке можно найти их разность (разность высот SRTM и топографической карты).
gdallocationinfo позволяет дать в точке значение растра (высоту).
Есть и плагин для QGIS Point sampling tool
Для получения более интересного результата взять вместо SRTM => Глобальные мозаики PALSAR
Вместо Global Forest Change treecover 2000 => лесоустроительные планшеты.
Например по Калиниградской области взять информацию из файла породы.tif
Зная информацию о преобладающей породе в квартале …
Насколько корректно использовать SRTM 2000 года? Насколько все выросло/вырублено за 15 лет?

[darsvid]

Зная значения высот по ЦМР и ЦММ в одной точке можно найти их разность.

В том-то и состоит проблема, что ЦМР у нас нет, нам ее нужно получить

Для получения более интересного результата взять вместо SRTM => Глобальные мозаики PALSAR

Спасибо, но интересует именно SRTM

Вместо Global Forest Change treecover 2000 => лесоустроительные планшеты.

Такой возможности у нас нет

Например по Калиниградской области взять информацию из файла породы.tif
Зная информацию о преобладающей породе в квартале …

У нас Украина

Насколько корректно использовать SRTM 2000 года? Насколько все выросло/вырублено за 15 лет?

Корректно, т.к. на интересует не лес, а рельеф. Оба набора данных GFC treecover и SRTM показывают состояние местности на 2000 г.

[bim2010]

del

[gamm]

Еще нашли вот такую публикацию, но кода нет и поэтому непонятно, что и как они там считали.

как считали, там написано, правда по английски, но формулы есть. Смысл в построении переходной зоны на границе леса (картинка на слайде 24), формула в статье. А дальше - константа внутри контуров леса, либо интерполяция с контура выдела. Интерполировать можно решением уравнения Пуассона, тогда внутри не будет локальных минимумов, и все будет гладенько (задав граничные условия немного за контуром, где точно леса нет).

[darsvid]

gamm, общая идея/ смысл использованного подхода понятны.

Но поскольку я не сильна в математике (и уж точно не программист), то я не знаю как эту формулу употребить на практике. Грубо говоря, как использовать формулу:

zDSM (δx δy) = a + a δx + a δy + a δxδy + hmG + ε

чтобы получить картинку как на слайде 26 презентации.

Какие инструменты ГИС-анализа и в какой последовательности в таком случае использовать.

[gamm]

Какие инструменты ГИС-анализа и в какой последовательности в таком случае использовать.

нужен программист среднего уровня на С++. Из готовых блоков в ГИС собрать можно, но довольно муторно, с ходу даже не напишу.

[darsvid]

Ок, спасибо, поробую спросить у знакомого программиста.

Методика кажется интересной и очень хочется попробовать воспроизвести, поэтому, если у вас все же получится дать наводки относительно "муторного" варианта будет очень-очень здорово.

PS
Спасибо всем откликнувшимся, просьба к всем читающим: если у вас есть какие-либо идеи по этому поводу - поделитесь

[gamm]

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

я просто перескажу кратко, что в статье написано, и что мы делали.

1) выделить границу, оценить ее выраженность. Оставить только выраженные фрагменты. Для выделения границы использовать карту лесов, корректируя ее по месту на рельеф.

2) гоним вдоль границы локальное окно (примерные параметры - в статье), и строим модельки.

Мы вместо (2) делали дистанционное преобразование от границы (1) внутрь контуров леса, и использовали модель зависимости высоты леса от расстояния до границы - были варианты сигмоида, "половинки" синуса, atanh, и т.д. - что лучше подойдет. Подбиралось фиксированное расстояние выхода на асимптоту, менялась только высота. Полученные (в локальных окнах) высоты леса интерполировались, и гладко сшивались с модельками около границы, там главная проблема на склонах, нужно поймать градиент склона, и его учесть (например, добавив к описанному линейную модель от координат). Но у вас вроде в основном равнина степная, проблем особых не должно быть.

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

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

Теоретически, все это можно реализовать скриптом, если прицепить код к локальному окну. Даже более того, для построения регрессии (линейная от координат и расстояния до границы) хватает вычисления среднего в локальном окне, и растровой арифметики, но скрипт будет длинный и с большим числом промежуточных растров (и нужно будет определители 4 порядка явными формулами считать).

Тема очень интересная, даже чисто методически. Было бы отлично "коллективным разумом" соорудить что-то общеупотребимое (например, скрипт). У меня, естественно, первым делом на ум приходит GRASS (r.mapcalc) + Python (или R?).

Может, выберем какой-то тестовый участок (или darsvid сможет предоставить данные)?
Например: SRTM 3'' + SRTM 1'' + ASTER GDEM + Global Forest Change + Landsat + отдешифрированный лес + топокарта (что-то их этого лишнее, наверняка).

[darsvid]

Александр,

спасибо за проявленный интерес к теме, данные (15 Мб) доступны по ссылке:

https://drive.google.com/file/d/0B1vjyq ... sp=sharing

файлы:
Hansen_GFC2014_treecover2000.tif - оригинальные данные GFC в проекции 4326, значения 0 устанолены как nodata
srtm1arcsec_32636.tif - SRTM 1'' перепроецированные в UTM с разрешением 30 м
treecover_4326.tif - маска лес (1)/ не лес (0) сделанная по оригинальным данным GFC. В данном случае лесом было назначено все в диапазоне 1-100%, но можно брать и другие значения (например, только больше 70%)
treecover_32636.tif - та же маска в UTM, значения пересчитаны Nearest Neighbor, параметры грида (экстент, разрешение) подогнаны под рельеф
ну и собственно шейп границ территории - суббасейн Десны 6 порядка (данныt HydroSHEDS).

Территория хороша тем, что там много леса и рельеф весьма разнообразный, но если слишком много - могу обрезать куски поменьше.

Если получится коллективными усилиями создать рабочее решение - обещаю все это оформить в виде статьи.

[gamm]

Тема очень интересная, даже чисто методически. Было бы отлично "коллективным разумом" соорудить что-то общеупотребимое (например, скрипт). У меня, естественно, первым делом на ум приходит GRASS (r.mapcalc) + Python (или R?).

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

[Boris]

Тема интересная, и мой вопрос для собственного развития, во второй части обсуждения я не увидел - используются топокарты или нет? Если речь шла о советских 100-ках, то на них с очень большой долей вероятности был именно рельеф, а не общие высоты.
В этом случае задача разбивается на две:
а) совместить карты и SRTM с приемлемой точностью, а есть ли лес сейчас ил его нет - это уже брать из других источников. Считать, что на картах леса - нет.
б) поднять рельеф с топокарт, веруя, что он навран не более, чем на 3-5 метров.

[bim2010]

darsvid хочу сначала поздравить вас и Александра с выходом новой книги QGIS By Example
Выкладываю ссылки на работы John Gallant и соавторов.
Gallant J Better quality DEMs from global elevation products
Gallant JC, Read AM and Dowling TI Building the national one second digital elevation model for Australia
стр 260
J.C. Gallant The ground beneath your feet: digital elevation data for today and tomorrow
Jenet Austin, John Gallant Comparison of Intermap 5 m DTM with SRTM 1 second DEM
J. C. Gallant, A. Read Enhancing the SRTM Data for Australia
John C. Gallant, Michael F. Hutchinson Producing digital elevation models with uncertainty estimates using a multi-scale Kalman filter
John Gallant Adaptive smoothing for noisy DEMs

Ссылки на готовые модули для расчета bare earth:
Module DTM Filter SAGA GIS
Wichmann, V.; Conrad, O.; Jochem, A. LiDAR Point Cloud Processing with SAGA стр26
Yury Ryabov DSM to DEM Conversion
Svidzinska, D. Methods of Geoecological Research: A Geoinformational Tutorial with the Open Source GIS SAGA.
стр 25

может пригодится:
Remove noisy spikes from LiDAR data using SAGA GIS and SRTM data
Convert a Land Cover Classification from Raster to Vector in QGIS

Готовая обработка Bare-Earth SRTM University of Bristol:
Development of a ‘bare-earth’ SRTM DEM product
download Bare-Earth SRTM University of Bristol

Вопрос: для расчета bare earth достаточно ли убрать леc? Ведь есть еще строения и т.п.

[dime1]

Команда Hydrosheds, оказывается, не смогла решить вопрос с лесом для hydrologically conditioned elevation.
Вот что пишут в спецификации:

The most significant characteristic is likely the fact that the elevation values of SRTM, being a radar-derived product, are influenced by the vegetation cover. In areas of low relief, these small deviations from the true surface elevation cancause significant errors in the derived river courses and flow directions.

Далее они просто берут основные реки в векторе и формируют искусственную долину рек процедурами Stream burning-Molding of valley courses-Sink filling-Carving through barriers, решив частичную задачу получения синтетических водотоков.
Потому сравнивая их модифицированные данные с сырыми SRTM1sec (по маске леса с GFC) можно увидеть, что где были лесные массивы (возвышенности относительно окружающего рельефа), там они и остались.