GIS-Lab - Вклад [ru]

Коррекция материалов Landsat

2015-09-26T18:24:08Z

Ericsson: Исправил стилистические и пунктуационные ошибки, убрал украинизмы

{{Статья|Черновик}}
{{Аннотация|В статье описывается упрощённая технология коррекции многозональных снимков Landsat: компенсация радиометрической неоднородности режима съёмки, а также пересчёт значений теплового канала в температуру}}
== Зачем нужна коррекция ==
Материалы Landsat в том виде, как они предоставляются пользователю, уже в значительной мере скорректированы, например, с учетом кривизны поверхности Земли, особенностей рельефа снимаемой территории, применено масштабирование в пределах возможных значений регистрируемых величин и т.п. Однако все эти коррекции направлены, главным образом, на обеспечение надёжного визуального дешифрирования наземных объектов.
Если же главным при анализе космоснимков становится не форма, рисунок и взаимное расположение объектов, а их индивидуальные спектральные характеристики, то требуется дополнительная коррекция изображений.
Типичным примером такого анализа является вычисление вегетационных индексов по многозональным материалам, в первую очередь – нормализованного разностного вегетационного индекса NDVI. При интерпретации схем распределения NDVI крайне важными становятся не относительные изменения значений по одной сцене, а соотношения NDVI по одной сцене или разновременным сценам в абсолютном выражении.
Например, возьмём две сцены:
* LE71870132004212EDC01 от 2004-07-30, зарегистрированную Landsat 7.
* LT51860132004213KIS00 от 2004-07-31, зарегистрированную Landsat 5.
Эти сцены имеют взаимное перекрытие, поэтому можно сравнить, как выглядит на них один и тот же участок территории возле Кольской АЭС. Правомерно предположить, что за одни сутки (с погрешностью в 12 минут, судя по времени регистрации) в состоянии растительных сообществ не должно было произойти никаких значимых изменений. Но схемы NDVI, построенные для этих сцен, демонстрируют нечто совсем иное (рис. 1).

[[Файл:LandsatTrans1.png|center|thumb|750px|Рис. 1. Сравнение двух схем NDVI, построенных по сценам, зарегистрированными разными спутниками Landsat с промежутком в 1 сутки]]

Сравнивая эти схемы, можно решить, что северная тайга в последний день июля 2004 года испытала быстрый расцвет. Судя по значениям NDVI, которые в среднем изменились на 18%, растительность повсеместно из категории угнетённой разреженной перешла в разряд пышной сомкнутой.
Этот курьёзный случай может быть объяснён только одним образом – нельзя сравнивать несравнимое. Прежде чем проводить вычисления NDVI, надо было выполнить коррекцию значений каналов каждой сцены в соответствии с калибровочными коэффициентами каждой из них.
Дело в том, что в данном случае сцены были засняты разными спутниками, аппаратура которых была настроена по-разному. Следовательно, зарегистрированные ими изображения без дополнительной радиометрической коррекции не годятся для тонкого спектрального анализа.

== Радиометрическая коррекция ==
Значения, содержащиеся в файлах каналов многозонального изображения, представляют собой безразмерную величину. Она пропорциональна интенсивности излучения, достигающего орбиты, на которой находится спутник Landsat<ref>Landsat 7 Science Data Users Handbook [Электронный ресурс] // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsathandbook.gsfc.nasa.gov/pdfs/Landsat7_Handbook.pdf. (дата обращения: 03.09.2015). – P. 117-120.</ref>. То, что мы получаем из файла – дискретное калиброванное значение пикселя <math> \ Q </math>, – номинировано в условных относительных единицах DN (Digital Numbers – числовые значения).
Требуется получить из <math> \ Q </math> значения отражающей способности (альбедо) <math> \ \rho </math> наземных объектов, видимых на космоснимке.
Для материалов Landsat до 7 включительно, полученных из архива ГС США, это производится следующим образом (1):

<center>
<math> \rho = \frac {\pi R ~ d^2} {E \cos \theta}, \qquad (1) </math>
</center>

где <math> \ R </math> – интенсивность излучения от объекта, достигшего орбиты Landsat, <math> \ d </math> – расстояние между Землёй и Солнцем, <math> \ E </math> – коэффициент светимости для каждого канала, <math> \ \theta </math> – высота стояния Солнца над горизонтом в момент съёмки.

Интенсивность излучения объекта, достигающего орбиты <math> \ R </math> вычисляется по формуле (2):

<center>
<math> R = \left ( R_{\mbox{max}} - R_{\mbox{min}} \right ) \frac {Q - 1} {254} + R_{\mbox{min}}, \qquad (2) </math>
</center>

Калибровочные коэффициенты <math> \ R_{\mbox{min}} </math> и <math> \ R_{\mbox{max}} </math> должны быть взяты из текстового файла метаданных с именем «*_mtl.txt», который поставляется вместе со сценой в одном архивном файле.

Если открыть этот файл для просмотра в текстовом редакторе, то там нужно найти строки, содержащие параметры, как показано на рис. 2. Из всего набора параметров следует выбрать коэффициенты, соответствующие определённому каналу.

[[Файл:LandsatTrans2.png|center|thumb|750px|Рис. 2. Содержимое файла с метаданными сцены Landsat 7, необходимыми для получения интенсивности излучения на сенсоре]]

Значение интенсивности излучения на сенсоре <math> \ R </math> имеет размерность Вт/м²/ср/мкм, т.е. мощность излучения, падающая на единицу площади земной поверхности сквозь телесный угол в 1 стерадиан, взятого относительно единицы длины волны излучения.

Расстояние между Землёй и Солнцем <math> \ d </math> может быть вычислено по приближённой формуле:

<center>
<math> d = 1 - 0{,}01668 \cos \left ( i \frac {2 \pi} {365} \right ) </math>
</center>

где <math> \ i </math> – порядковый номер в году дня получения изображения. Его можно определить по дате регистрации сцены, которая также находится в файле метаданных (рис. 3):

[[Файл:LandsatTrans3.png|center|thumb|750px|Рис. 3. Метаданные сцены: дата регистрации (DATE_ACQUIRED)]]

Если для вычислений использовать MS Excel, то в нём для получения значения <math> \ i </math> удобно использовать формулу <code> =Ref-ДАТА(ГОД(Ref);1;1)</code>, где Ref – ссылка на ячейку, где записана дата регистрации.

Величина <math> \ d </math> измеряется в астрономических единицах (а.е.).

Значения светимости <math> \ E </math> для каждого канала могут быть взяты из следующей таблицы (табл. 1):

<center>
'''Табл. 1. Коэффициенты светимости для каналов Landsat, Вт/м²/мкм'''

{| class="wikitable"
|-
! Каналы !! 1 !! 2 !! 3 !! 4 !! 5 !! 7 !! 8
|-
| Landsat 4 || 1957 || 1825 || 1557 || 1033 || 214.9 || 80.72 || –
|-
| Landsat 5 || 1957 || 1826 || 1554 || 1036 || 215.0 || 80.67 || –
|-
| Landsat 7 || 1970 || 1842 || 1547 || 1044 || 225.7 || 82.06 || 1369
|}
</center>

И последний, необходимый для вычисления альбедо, параметр <math> \ \theta </math> – высота стояния Солнца над горизонтом в момент съёмки – может быть взят из файла с метаданными сцены (рис. 4). Значение это в файле приведено в градусах. Перед подстановкой в формулу его надо перевести в радианы.

[[Файл:LandsatTrans4.png|center|thumb|750px|Рис. 4. Метаданные сцены: высота Солнца (SUN _ELEVATION)]]

Для материалов Landsat 8, взятых из архива ГС США, расчётные формулы имеют гораздо более простой вид:

<center>
<math> \rho = \frac {R} { \sin \theta}, \qquad (3) </math>
</center>

где

<center>
<math> R = 2 \cdot 10^{-5} Q - 0{,}1. \qquad (4) </math>
</center>

После вычисления по представленным формулам значения альбедо <math> \ \rho </math> по разным каналам могут быть использованы, например, для определения величины NDVI:

<center>
<math> \mbox{NDVI} = \frac {\rho_{_{\mbox{NIR}}} - \rho_{_{\mbox{RED}}}} {\rho_{_{\mbox{NIR}}} + \rho_{_{\mbox{RED}}}}. </math>
</center>

Здесь NIR – номер канала ближнего инфракрасного диапазона, RED – номер канала красного диапазона. Определить какие номера каналов соответствует эти диапазонам можно с помощью табл. 2.

<center>
'''Табл. 2. Соответствие между диапазонами электромагнитного излучения и номерами каналов спутников разных поколений миссии Landsat<ref>The Multispectral Scanner System // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=3227. (дата обращения: 13.09.2015).</ref><ref>The Thematic Mapper // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=3229. (дата обращения: 13.09.2015).</ref><ref>The Enhanced Thematic Mapper Plus // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=3225. (дата обращения: 13.09.2015).</ref><ref>Operational Land Imager // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=5107. (дата обращения: 13.09.2015).</ref><ref>Thermal Infrared Sensor // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=5112. (дата обращения: 13.09.2015).</ref>'''

{| class="wikitable"
|-
! Диапазон !! Landsat 1 !! Landsat 2 !! Landsat 3 !! Landsat 4 !! Landsat 5 !! Landsat 7 !! Landsat 8
|-
| RED || 5 || 5 || 5 || 3 || 3 || 3 || 4
|-
| NIR || 6 || 6 || 6 || 4 || 4 || 4 || 5
|-
| Тепловой || – || – || 8 || 6 || 6 || 61, 62 || 10, 11
|}
</center>

== Определение температуры земной поверхности ==
По значениям тепловых каналов (табл. 2) можно определить температуру подстилающей поверхности. Теоретически точность оценки температуры около 0,5°С, однако дымка в атмосфере занижает значения на несколько градусов.

Исходными данными для определения температуры служат значения интенсивности излучения <math> \ R </math>, пришедшего на сенсор спутника и зарегистрированного соответствующим тепловым каналом, в соответствии с методикой, описанной в предыдущем разделе (формулы 2, 4).

Для Landsat 4-5 температура <math> \ T </math> вычисляется по следующей формуле:

<center>
<math> T = \frac {K_2} { \ln{ \left( \frac {K_1} {R} + 1 \right) }} - 273{,}15 \ (^{\circ} \mbox{C}), \qquad (5) </math>
</center>

Здесь <math> \ K_1 </math> и <math> \ K_2 </math> – калибровочные константы, значения которых можно взять из табл. 3.

<center>
'''Табл. 3. Калибровочные константы для вычисления температуры, Вт/м²/мкм'''

{| class="wikitable"
|-
! Диапазон !! Landsat 4 !! Landsat 5 !! Landsat 7
|-
| <math> \ K_1 </math> || 671.62 || 607.76 || 666.09
|-
| <math> \ K_2 </math> || 1284.30 || 1260.56 || 1282.71
|}
</center>

Пример результатов определения значений температуры подстилающей поверхности по данным теплового канала показан на рис. 5.

[[Файл:LandsatTrans5.png|center|thumb|750px|Рис. 5. Тепловая аномалия города Воронеж в марте 1985 г. Значения шкалы – °С]]

Материалы Landsat 7 содержит два тепловых канала: 61 и 62, которые отличаются настройками усиления сигнала. Для вычисления температуры можно использовать в формуле (5) любой из них или оба – они дают почти одинаковые значения (±0,1°С).

Почти также и с помощью той же формулы (5) производится вычисление температуры по данным Landsat 8, за исключением того, что значения констант <math> \ K_1 </math> и <math> \ K_2 </math> надо брать отдельно для каждого теплового канала данной сцены из её файла метаданных (рис. 6).

[[Файл:LandsatTrans6.png|center|thumb|750px|Рис. 6. Метаданные сцены Landsat 8: калибровочные константы для вычисления температуры]]

Значения температуры, определённые по каналам 10 и 11 (они различаются охватываемыми интервалами длин волн теплового диапазона) Landsat 8, отличаются друг от друга на 1,5-3°С. В общем случае можно усреднить их.

== Литература ==
<references />

== См. также ==
* [http://gis-lab.info/qa/ss.html Технические характеристики инструментов ДЗЗ и их носителей]
* [http://gis-lab.info/qa/landsat-glovis.html Получение бесплатных космических снимков Landsat TM,ETM+ через Glovis]
* [http://gis-lab.info/qa/earthexplorer-work.html Работа с архивом материалов ДЗЗ через EarthExplorer]
* [http://gis-lab.info/qa/dn2radiance.html Конвертация данных TM, ETM+ в показатели излучения на сенсоре]
* [http://gis-lab.info/qa/dn2temperature.html Конвертация данных Landsat TM/ETM+ в значения температуры - Теория]
* [http://gis-lab.info/qa/earthsundist.html Определение расстояния от Земли до Солнца в момент получения снимка]
* [http://gis-lab.info/qa/ndvi.html NDVI - теория и практика]
[[Участник:Константин Силкин|Константин Силкин]] 22:34, 24 сентября 2015 (MSD)

Формат ASCIIGRID с примерами

2012-01-27T19:32:41Z

Ericsson:

{{Статья|Черновик}}
Формат ASCIIGRID это форма представления регулярно-сетевых данных (растров). Исходно этот формат использовался как обменный в ГИС ARC/INFO (файлы в этом формате создавались командой GRIDASCII), сейчас его основная часть поддерживается всеми основными ГИС (QGIS, ArcGIS). Файл в формате ASCIIGRID представляет собой текстовый файл в котором содержимое разделяется на две части: 1) заголовок, где записывается общая информация о растре, его размеры, географическая привязка и т.д., 2) собственно, секция данных.

Данные в этом формате удобны прежде всего удобством быстрого визуального анализа содержимого: много-мало данных, какие значения, есть ли nodata и т.д.. Разумеется, в этом формате не рекомендуется хранить большие растры.

Поддерживается чтение и запись в этот формат с помощью GDAL.

=== Заголовок ===

<code>ncols</code> - количество колонок (целочисленное); 
<code>nrows</code> - количество рядов (целочисленное); 
<code>xllcorner</code> - х-координата левого нижнего угла, долгота (десятичные градусы DD.D(D), DD:MM:SS не поддерживается) или смещение на восток (метры) (целочисленное или с плавающей точкой) 
<code>yllcorner</code> - y-координата левого нижнего угла, широта или смещение на север 
<code>cellsize</code> - линейный размер квадратной ячейки (целочисленное или с плавающей точкой) 
<code>nodata_value</code> - значение "нет данных", этот параметр не обязателен, он нужен для того, чтобы программное обеспечение автоматически делало все ячейки с этим значением прозрачными (целочисленное или с плавающей точкой). Часто используемое значение -9999, нецифровые значения, например "NA" использовать нельзя 

=== Данные ===

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

=== Примеры ===
Все примеры можно копировать и вставлять в текстовые файлы с расширением *.asc, они должны работать. Чтобы получить такие же иллюстрации в QGIS, нужно: 1. добавить растр, 2. зайти в свойства\Стиль и изменить Цветовая карта на Псевдоцвета, 3. зайти в свойства\Прозрачность и удалить прозрачный пиксел со значением -9999, а также заменить в значении no data -9999 на 0.
==== Стандартный заголовок ====

[[Файл:Raster example 8x6.gif|400px|Пример растрового изображения в псевдоцветах]]

Это же изображение в табличном виде:
{|class="wikitable"
|-
| -9999   ||-9999   ||-2080   ||-1909   ||-9999   ||-601    ||-9999   ||-1362   
|-
| -510||-387||-350||-182||-9999||-9999||-9999||27
|-
| 13 || -92 || -9999 || -33 || 54 || 28 || -182 || -750
|-
| 59 ||-165 ||18 ||167 ||-182 ||-9999 ||-608 ||666
|-
| -62 ||-9999 ||-805 ||-9999 ||-9999 ||-9999 ||-9999 ||-497
|-
| -9999 ||-9999 ||332 ||1240 ||-9999 ||-9999 ||-9999 ||-9999
|}

И непосредственно в формате ASCIIGRID
<pre>ncols 8
nrows 6
xllcorner -40.5
yllcorner 25.25
cellsize 10
NODATA_value -9999
-9999 -9999 -2080 -1909 -9999 -601 -9999 -1362
-510 -387 350 -182 -9999 -9999 -9999 27
13 -92 -9999 -33 54 28 -182 -750
59 -165 18 167 -182 -9999 -608 666
-62 -9999 -805 -9999 -9999 -9999 -9999 -497
-9999 -9999 332 1240 -9999 -9999 -9999 -9999</pre>

==== Неквадратная ячейка ====

В примере приведенном выше, ячейка растра квадратная (т.е. высота равна широте). Если необходима прямоугольная ячейка, где высота не равна широте, вместо параметра cellsize можно использовать DX и DY, например:

<pre>
ncols 8
nrows 6
xllcorner -40.5
yllcorner 25.25
dx 11.6
dy 5.025
NODATA_value -9999
-9999 -9999 -2080 -1909 -9999 -601 -9999 -1362
-510 -387 350 -182 -9999 -9999 -9999 27
13 -92 -9999 -33 54 28 -182 -750
59 -165 18 167 -182 -9999 -608 666
-62 -9999 -805 -9999 -9999 -9999 -9999 -497
-9999 -9999 332 1240 -9999 -9999 -9999 -9999</pre>
ASCIIGRID файлы с параметрами DX и DY вместо CELLSIZE не поддерживаются многими пакетами GIS (в т.ч. ArcGIS), но поддерживаются QGIS.

Вот как будет выглядеть файл выше, где cellsize заменен на DX, DY:

[[Файл:Raster example 8x6 dx-dy.gif|400px|Пример растрового изображения в псевдоцветах]]

==== Привязка по центру ячейки, а не углу ====
Альтернативой привязки по углу, т.е. использования <code>xllcorner</code> и <code>yllcorner</code> является привязка по центру ячейки, в этом случае нужно использовать <code>xllcenter</code> и <code>yllcenter</code>. Результатом, при неизменных координатах будет сдвиг растра на половину ячейки вниз и влево (по диагонали).

<pre>
ncols 8
nrows 6
xllcenter -40.5
yllcenter 25.25
cellsize 10
NODATA_value -9999
-9999 -9999 -2080 -1909 -9999 -601 -9999 -1362
-510 -387 350 -182 -9999 -9999 -9999 27
13 -92 -9999 -33 54 28 -182 -750
59 -165 18 167 -182 -9999 -608 666
-62 -9999 -805 -9999 -9999 -9999 -9999 -497
-9999 -9999 332 1240 -9999 -9999 -9999 -9999</pre>

Если взять один и тот же растр и заменить xllcenter,yllcenter на xllcorner,yllcorner, а сверху наложить сетку границ пикселей с первого растра, то это будет выглядеть так:

[[Файл:Raster example 8x6 xllcenter.gif|400px|Пример сдвига]]

=== Практика ===
Таким образом, чтобы создать растр размером m*n не используя никаких инструментов кроме текстового редактора и табличного процессора, нужно сделать следующее.

# Создать таблицу нужного размера в табличном процессоре (OpenOffice Calc, MS Excel и т.д.) или любом редакторе. Скопировать и вставить в текстовый документ. Заменить табуляции на пробелы.
# Прописать заголовок файла в описанном выше формате. Еще один пример для матрицы 360*180, представляющий массив с разрешением 1 градус на весь мир:

<pre>ncols 360
nrows 180
xllcorner -180.000000000000
yllcorner -90.000000000000
cellsize 1.000000000000
nodata_value 0
12 321 1 12 3 4 1 23 4 1 234 1234 1 2 3 4 2 1 2 3 ...</pre>

# Поменять/задать расширение текстового файла на asc.
# Полученный файл можно открыть в QGIS и экспортировать в любой другой растровый формат или работать прямо с ним.