Обрезка Landsat

[Nule]

Здравствуйте!
Задача - пакетно обрезать снимки Landsat, с целью убрать артефакты с боков.
Т.е. надо наложить маску такого вида:

Причем маска должна быть "универсальной", т.е. не зависеть от координат обрабатываемой сцены.
Можно ли это реализовать в GRASS?
Спасибо!

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

Может, и можно. Нужны примеры снимков с "артефактами" (устно не совсем понятно).

[Nule]

Может, и можно. Нужны примеры снимков с "артефактами" (устно не совсем понятно).

Примером может служить абсолютно любой снимок landsat.

[ericsson]

Здравствуйте!
Причем маска должна быть "универсальной", т.е. не зависеть от координат обрабатываемой сцены.

А кто сказал, что маска может быть универсальной? Вот browse-файлы из архива GLCF, одного спутника и одного класса обработки. Похоже, но не совпадает:
ftp://ftp.glcf.umd.edu/glcf/Landsat/WRS ... browse.jpg
ftp://ftp.glcf.umd.edu/glcf/Landsat/WRS ... browse.jpg
ftp://ftp.glcf.umd.edu/glcf/Landsat/WRS ... browse.jpg

Я бы при решении задачи обрезки мусора сделал бы следующее:
Для каждой сцены извлек бы незначащие точки из каждого канала в композите (можно даже из уменьшенной в два или четыре раза копии), инвертировал, нашел бы пересечение получившихся масок значащих данных, инвертировал лезультат, (увеличил бы во столько раз во сколько уменьшал в начале), добавил бы некоторый threshold и по получившемуся - обрезал бы композит. GRASS не владею, но раз уж это можно сделать в ImageMagick (коим пользуюсь для обычной графики в т.ч. в подобных задачах), то и в GRASS, наверняка, можно.

[Nule]

В ERDAS IMAGINE есть универсальная маска в инструментах Mosaic Pro (маску можно нарисовать вручную) и Mosaic Express (можно задать, сколько процентов срезать). Но это инструменты для создания мозаик и использовать их для обрезания отдельного снимка неудобно.
Я надеялся, что может есть аналогичные, но более гибкие инструменты в других пакетах.

[ericsson]

Ну так а что может быть более гибким, чем пара-тройка булевых операций над растрами?

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

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

Здорово было бы увидеть этот метод в действии с помощью IM. Хотя бы на примере приведённого выше увеличенного куска снимка. А самое интересное -- совмещение модулей GRASS + утилит IM.

[ericsson]

Постараюсь сделать в ближайшие дни.

[ericsson]

Ну собственно, как обещал...
Требуемые средства - listgeo, geotifcp (из Libgeotiff), convert (свежий ImageMagick).

1) Сделать дамп метаданных любого канала в отдельный текстовик через listgeo

2) Для каждого канала, который предполагается использовать в композите, выполняем команду такого вида:

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

convert L71169028_02820060620_B10.TIF -alpha set -channel RGBA -fuzz 5% -fill none -floodfill +0+0 black L71169028_02820060620_B10.withalpha.TIF

Параметр fuzz отвечает за то чтобы захватить не совсем черные края.
На выходе получаем тот же файл с убитыми метаданными (привязкой) но двухканальный - восьмибитный greyscale, плюс маска там, где был черный ±5%.
Примечание: если используются "дурные" данные с вырезанными полосами и прочим безобразием (которые с поврежденного спутника, если я не путаюсь в причине их возникновения), можно использовать другой код:

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

convert L71169028_02820060620_B10.TIF  -transparent black  L71169028_02820060620_B10.withalpha.TIF

Он пропишет прозрачность черным точкам не только по краям растра, но и изолированным незначащим точкам среди значащих.

3) Отделяем маску командой

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

convert.exe L71169028_02820060620_B10.withalpha.TIF -alpha extract L71169028_02820060620_B10.alpha.TIF

И получаем однобитные маски в виде отдельных файлов.

4) Попарно (1+2)+3 складываем маски командами типа

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

convert L71169028_02820060620_B10.alpha.TIF L71169028_02820060620_B20.alpha.TIF -compose Darken -composite L71169028_02820060620_B10_B20.alpha.TIF

convert L71169028_02820060620_B10_B20.alpha.TIF L71169028_02820060620_B30.alpha.TIF -compose Darken -composite L71169028_02820060620_B10_B20_B30.alpha.TIF

Наложение растров Darken эквивалентно логическому пересечению, если незначащие части - черные, значащие - белые.

5) Убиваем получившиеся файлы *_withalfa.TIF - они больше не нужны.

6) Через geotifcp вернуть файлу L71169028_02820060620_B10_B20_B30.alpha.TIF привязку из дампа, который был сделан на первом шаге. Далее эту маску можно использовать где угодно вместе с исходными растрами.

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

Спасибо большое, ericsson. Будем пробовать

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

Создание подобной маски для всех каналов собственно в GRASS (взято с отсюда):

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

BASE=L71074092_09220040924

 

g.region rast=$BASE.1

r.series in=`g.mlist pat="$BASE.[0-8]*" sep=,` -n out=$BASE.thresh method=threshold thresh=1

r.mapcalc "$BASE.mask = if(isnull($BASE.thresh))"

g.remove $BASE.thresh

Проверено -- работает, причём довольно шустро.
На заметку: надо следить за именами растров (КО: они могут быть не как в примере).

P.S. А вот как вообще обрезать "бахрому" на автомате?

[ericsson]

Ну еще можно вооружиться полигонами сетки path/row (где-то в недрах архива GLCF такие шейпы были), извлекать соответствующие номера из имен файлов или сопутствующих метаданных и обрезать по этим полигонам.

[Nule]

Странно, казалось бы для такой массовой вещи как ландсат должна быть куча разнообразнейших тулзовин.
Что-нибудь типа: находим самый верхний пиксель, отличный от noData, потом нижний, правый и левый, по ним строим векторный прямоугольник, уменьшаем его на 10%, обрезаем по шейпу.

[rhot]

Создание подобной маски для всех каналов собственно в GRASS (взято с отсюда):

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

BASE=L71074092_09220040924

 

g.region rast=$BASE.1

r.series in=`g.mlist pat="$BASE.[0-8]*" sep=,` -n out=$BASE.thresh method=threshold thresh=1

r.mapcalc "$BASE.mask = if(isnull($BASE.thresh))"

g.remove $BASE.thresh

Проверено -- работает, причём довольно шустро.
На заметку: надо следить за именами растров (КО: они могут быть не как в примере).

P.S. А вот как вообще обрезать "бахрому" на автомате?

Объясните, пожалуйста, что это за метод threshold?

Странно, казалось бы для такой массовой вещи как ландсат должна быть куча разнообразнейших тулзовин.

А зачем? Таких задач бывает так же много, как пользователей ландсат. ИМХО, программа не должна быть перегружена всякими примочками и тулзовинами. GRASS в этом смысле очень позитивно смотрится; исходники - чуть больше 20 Мб.
Хочешь, качай и устанавливай доп. модули, а хочешь, сам пиши на питоне - свобода

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

Объясните, пожалуйста, что это за метод threshold?

Честно говоря, для меня это тоже не до конца прояснилось. Кроме мануала, пояснения от автора модуля (Glynn Clements) есть, например, здесь.