Ошибка в QUAC или band math (Landsat5, ENVI)

[Natalia Novoselova]

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

Не догадалась. Во вложении

L5_b1b2b3b4b5 - hdr скленного в один файл пяти каналов после изменения параметров в Edit ENVI Header
quac_b1_5_test7- hdr файла, прошедшего коррекцию QUAC. Зашла туда, там и правда что-то непонятно насчет 4-го канала указано. В сам низу почему на 4-м месте стоит ноль?

0.450000, 0.520000, 0.630000, 0.760000, 1.750000}
bbl = {
1, 1, 1, 0, 1}

Каналы файла L5_b1b2b3b4b5 - нормальные. Но после коррекции QUAC канал band4 становится, вдруг, равным нулю... Потому с ndvi и не получается

[Natalia Novoselova]

Каналы файла L5_b1b2b3b4b5 - нормальные. Но после коррекции QUAC канал band4 становится, вдруг, равным нулю... Потому с ndvi и не получается.

[ericsson]

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

Для начала - можно было такие короткие файлы не архивировать и прикладывать во вложении, а вставить в текст:
L5_b1b2b3b4b5.hdr

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

ENVI
description = {
  Create Layer File Result [Thu Oct 09 16:09:49 2014]}
samples = 7921
lines   = 6961
bands   = 5
header offset = 0
file type = ENVI Standard
data type = 1
interleave = bsq
sensor type = Landsat TM
byte order = 0
map info = {UTM, 1.000, 1.000, 228585.000, 7549315.000, 3.0000000000e+001, 3.0000000000e+001, 23, South, WGS-84, units=Meters}
coordinate system string = {PROJCS["UTM_Zone_23S",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",10000000.0],PARAMETER["Central_Meridian",-45.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]}
wavelength units = Micrometers
pixel size = {3.00000000e+001, 3.00000000e+001, units=Meters}
band names = {
 Layer (Band 1:LT52190762011111CUB00_B1.TIF), 
 Layer (Band 1:LT52190762011111CUB00_B2.TIF), 
 Layer (Band 1:LT52190762011111CUB00_B3.TIF), 
 Layer (Band 1:LT52190762011111CUB00_B4.TIF), 
 Layer (Band 1:LT52190762011111CUB00_B5.TIF)}
wavelength = {
 0.450000, 0.520000, 0.630000, 0.760000, 1.750000}

quac_b1_5_test7.hdr

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

ENVI
description = {
  QUick Atmospheric Correction Result [Thu Oct 09 16:26:39 2014]}
samples = 7921
lines   = 6961
bands   = 5
header offset = 0
file type = ENVI Standard
data type = 2
interleave = bsq
sensor type = Landsat TM
byte order = 0
map info = {UTM, 1.000, 1.000, 228585.000, 7549315.000, 3.0000000000e+001, 3.0000000000e+001, 23, South, WGS-84, units=Meters}
coordinate system string = {PROJCS["UTM_Zone_23S",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",10000000.0],PARAMETER["Central_Meridian",-45.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]}
wavelength units = Micrometers
reflectance scale factor = 10000.000000
band names = {
 QUAC (Layer (Band 1:LT52190762011111CUB00_B1.TIF):L5_b1b2b3b4b5), 
 QUAC (Layer (Band 1:LT52190762011111CUB00_B2.TIF):L5_b1b2b3b4b5), 
 QUAC (Layer (Band 1:LT52190762011111CUB00_B3.TIF):L5_b1b2b3b4b5), 
 QUAC (Layer (Band 1:LT52190762011111CUB00_B4.TIF):L5_b1b2b3b4b5), 
 QUAC (Layer (Band 1:LT52190762011111CUB00_B5.TIF):L5_b1b2b3b4b5)}
wavelength = {
 0.450000, 0.520000, 0.630000, 0.760000, 1.750000}
bbl = {
 1, 1, 1, 0, 1}

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

Сейчас будем разбираться - я этот формат вижу первый раз в жизни.

[ericsson]

Очевидно, что разница между файлами - в появлении в обработанном QUAC строк
16:

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

reflectance scale factor = 10000.000000

25,26:

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

bbl = {
 1, 1, 1, 0, 1}

Вбиваем в поиск "ENVI hdr format", находим http://www.exelisvis.com/docs/ENVIHeaderFiles.html , идем в секцию The ENVI Header Format и видим:

bbl = Lists the bad band multiplier values of each band in an image, typically 0 for bad bands and 1 for good bands.

Еще там есть ссылка на http://www.exelisvis.com/docs/enteropti ... tml#Select где объясняется, для чего служит параметр bad bands:

Use Bad Bands List to select bands to exclude from plotting or optionally omit during processing. The Bad Bands list is often used to omit the water vapor bands in hyperspectral datasets.

• In the Header Info dialog, click Edit Attributes and select Bad Bands List. The Edit Bad Bands List values dialog appears.
• All bands in the list are highlighted by default as good. Deselect any desired bands in to designate them as bad bands.
• To specify a range of bands, enter the beginning and ending band numbers in the fields next to the Add Range button. Click Add Range.
• Click OK.

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

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

[ericsson]

Вбиваю в поиск "ENVI QUAC" и нахожу http://www.exelisvis.com/docs/QUAC.html
Читаем:

Before You Begin

The input data for QUAC can be in radiance, apparent reflectance, or raw/uncalibrated units. You should apply a mask to images with large borders or background regions prior to using the image in QUAC. Input files must have at least three bands and valid wavelengths. For multispectral and hyperspectral sensors, wavelength values must be available in the ENVI header file or in separate ASCII files.
Using QUAC

• From the Toolbox, select Radiometric Correction > Atmospheric Correction Module > QUick Atmospheric Correction (QUAC). The QUick Atmospheric Correction Input File dialog appears.
• If the file is not already listed, click the Open drop-down button and select an input file. Then click Open.
• Perform any optional spatial or spectral subsetting or apply a mask.
• Click OK in the QUick Atmospheric Correction Input File dialog. The QUick Atmospheric Correction Parameters dialog appears.
• Specify a sensor type. If you do not know the sensor type of the input file or the sensor type is not available in the drop-down list, select Unknown. ENVI will handle the input file in a generic way.
• Select output to File or Memory.
• Click OK to process the file. QUAC creates a surface reflectance image, scaled into two-byte signed integers using a reflectance scale factor of 10,000.

If the mean spectra of endmembers extracted from the image has a very low response in a certain band, that band is set as a "bad band." If this occurs, check the image for bad pixels or those with background values (especially in Landsat data). These pixels bias the calculation of endmember mean spectra. If the image has bad pixels or pixels with background values, you should build a mask to ignore those pixels, then run QUAC with the mask.

You can also write a script to perform QUAC processing by using the ENVIQUACTask routine.

Собственно, как я и ожидал: в документации все расписано буквально для младшего школьного возраста (важные места я подчеркнул):
То, что появляется в 16-й строке - это значение, которое по умолчанию вбито в это средство, так что вопрос о ее смысле снимается.

А вот про строку 26 все оказывается очень просто: в первом подчеркнутом предложении вам человеческим языком сказано, что вы должны оставить в каналах только значащие данные, то есть воспользоваться маской или чем-либо еще для исключения незначащих пикселей из рассмотрения. В описании формата HDR мне это не попалось, но есть http://www.exelisvis.com/docs/EnterOpti ... tml#Enter6 (нашлось также поиском по "ENVI no data value"). Далее, в заключении, вам еще раз сказано, что если диапазон выглядит для программы слишком темным (а у Ландсата с фоновыми пикселями так и выйдет всегда - там нули), то он получит флаг bad band.

Тут снова всплывает вопрос о пользе образования: вам в голову не приходило, что работая с набором данных, где есть незначащие области, эти области из рассмотрения нужно исключать хотя бы по соображениям того, чтобы не получить там новые мусорные значения и не нагружать машину вычислением того, чего нет? Подобная догадка - сама собой разумеющаяся вещь для тех, кто задумывается о том, что стоит за методами. А для нажимателей кнопок по пошаговым инструкциям она таковой никогда не станет. Я не буду нарушать правила, высказывая сомнения в вашей профпригодности, однако позволю себе, тем не менее, шутку о том, что в своих "научных" изысканиях вам стоило бы добавить кое-кого из присутствующих на этом форуме в соавторы. В научных кругах принято указывать в списке авторов работ тех, кто их действительно выполнял.

На всякий случай, посмотрел еще кое-какие параметры в документе по формату:
data type = 2 соответствует типу данных "Integer: 16-bit signed integer".

[ericsson]

Легкий оффтопик об образовании: у меня оно инженерное. Все, что я знаю о ГИС, геодезии, картографии и так далее - результат моего банального любопытства. Хотя мне это позволяет писать технические статьи (не здесь, правда - выбираю более посещаемые площадки). Так что "профильное" образование на географических факультетах и в "НИИ Гаек" для того, чтобы обладать элементарной логикой и иметь базовые технические знания и так далее - не нужно.

[Natalia Novoselova]

Спасибо за поиск. Появился еще один (невысокий) шанс. Невысокий потому, что один раз я уже делала точно так и получалось (с полями). И также еще один нелегкий вопрос - научиться строить маски в Энви. Для меня там каждый вопросик, даже небольшой на вид, трудная в прохождении гора. Именно по причине такого хелпа (в отличие от хелпа ArcGIS)

Ну вот о чтении хелпа. Для получения дробных зачений нужно использовать float
Читаю в хелпе - приводится образе float ()

Считаю float (B1 - 0.3) - получается (B1 - канал)

Считаю логичное построение - float((B1-B2)/(B1+B2)) - не получается.


Почему? И где в хелпе искать ответ на вопрос, почему обычная логика математики там не работает и что в итоге надо делать?

[ericsson]

Исчерпывающая информация по созданию масок содержится здесь http://www.exelisvis.com/docs/ManageRas ... ml#Masking плюс, вы уже проигнорировали ссылку, которую я дал выше - http://www.exelisvis.com/docs/EnterOpti ... tml#Enter6 и уже ноете, как все сложно. Это наводит на мысль о том, что никакая помощь вам не в прок. Разубедите нас.

Для получения дробных значений чего и зачем?
Если вам нужен NDVI, почему не воспользоваться http://www.exelisvis.com/docs/NDVI.html ?

[Natalia Novoselova]

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

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

Для получения дробных значений чего и зачем?

А для Band rations? Это же общее для расчета. Но как, если не по отдельности считать, не известно. Логичное построени формулы - не работает

[ericsson]

Все что вы написали про тупик - пустое нытье. Перечитайте выше все, что я написал о проблеме с масками - это все предельно ясно и понятно, искать и читать это было быстрее, чем описывать это здесь. Вместо того, чтобы учиться понимать, вы только тренируетесь повторять. Это не то, что я имею в виду.

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

Если под Band rations вы подразумеваете band ratios http://www.exelisvis.com/docs/BandRatios.html то спрошу: а вы уверены, что на входе эта операция принимает только данные во float?

[alexandr cherepanov]

Пока в header QUAC пишет строку bbl c 0 будут ошибки расчета всего остального.

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

bbl = {
 1, 1, 1, 0, 1}

У Вас это означает, что 4 в списке канал (NIR) битый. Попробуйте создать маску для черных полей как уже писали выше.

Считаю логичное построение - float((B1-B2)/(B1+B2))

Если 4 канал полностью битый со значением 0, то при любом значение 3 канала получите по этой формуле -1.000000 или 0.000000.
И даже если избавитесь от ошибки с 4 каналом, то с такой формулой все равно на выходе будут только -1.000000, 0.000000 и 1.000000. Почему понятно?

[ericsson]

Бегло прошелся по http://www.exelisvis.com/docs/BandMath.html и вычитал следующее:

Note that, in this example, the IDL function float() is used to prevent byte overflow errors during calculation.

Выше дан пример

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

(float(b1)+float(b2)+float(b3))/3.0

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

Это - знания из области программирования, самые азы: если у нас есть, например, восьмибитная целочисленная беззнаковая переменная, то она может принимать значения 0..255, а если их две, и вы умножите одну на другую, а каждая будет иметь значение, например, 200 - вы получите 40000, что никак нельзя записать теми же восемью битами.

Потому ваше

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

(B1-B2)/(B1+B2)

может дать вам полную хрень на выходе. Чтобы этого не случилось, документация советует преобразовать типы во float внутри формулы в первую очередь, а не в последнюю, как пытаетесь поступить вы совершенно бездумно (вот зачем?).
И формула должна выглядеть, как минимум, так (вычитание не должно вызывать переполнений, если типы переменных - не беззнаковые):

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

(float(B1-B2))/(float(B1+B2))

а если это не сработает (мне лень копаться в документации по IDL, но вдруг аргументом float() не может служить выражение, а может - только переменная?), то вот так:

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

(float(B1)-float(B2))/(float(B1)+float(B2))

Если все получится - с вас мешок кофе. Потому что я ENVI никогда не использовал.

[alexandr cherepanov]

viewtopic.php?f=20&t=5620
9 пост
Оба варианты правильные, второй короче.

[ericsson]

Александр, судя по (float(b4)-b3)/(float(b4)+b3) - в ходу неявное преобразование типов? Мне стало любопытно, в онлайн-документации по IDL про это действительно не нашел ничего. Тут http://www.exelisvis.com/docs/routines-89.html только про явное.

[alexandr cherepanov]

Да в IDL тип переменной можно менять в любой момент. Плюс еще всякие фишки.
https://www.exelisvis.com/docs/Dynamic_Data_Typing.html

Using the IDL data type casting functions, like fix(), is not the only way that the data type of your Band Math result can get promoted. This is because IDL is Dynamically Typed, which means that the data type of an IDL statement is automatically promoted to the largest data type it encounters in the expression

Для расчетов все входные данные не обязательно переводить. Нужно просто следить за типом данных на входе и диапазоном значений.
Если например значения типа integer и b1=200, а b2=240, то что float(b1)-b2, что b1-float(b2) дадут на выходе одно и тоже число -40.000000. Если запишем выражение по другому float(b1-b2), тоже получим -40.000000. Разность в диапазоне допустимых значения типа integer.
А вот если входные данные будут типа byte, с первыми двумя выражениями проблем не будет, а вот float(b1-b2) даст на выходе не -40.000000, а 216.000000. Поскольку переполнения произошло до перевода в float point.

Если взять наши конкретные данные, то после QUAC/FLAASH диапазон данных должен быть 0-10000, так как это значения коэффициентов отражения умноженные на 10000. Все что за пределами можно считать ошибками и не использовать в расчетах. Тип данных integer возможные значения -32768 to + 32767.
Наша формула расчета NVDI = (b1-b2)/(b1+b2), числитель может принимать значения от -10000 до 0, знаменатель от 0 до 10000. Переполнения не происходит. Чтобы получить правильные дробные значения хватит одного float() в любом месте - либо для одного b, либо для числителя или знаменателя. Хоть даже так (b1-b2)/(b1+float(b2)).

Чтобы не думать проще использовать перевод в float() и в числителе и в знаменателе.
На выходе точно также для уменьшения занимаемого места можно оставить нужное число знаков умножив на 1000 или 10000 и перевести обратно в integer fix(b10*10000).