TALLER TELEDETECCIÓN¶

Fundamentos de la teledetección¶

Esquema de un sistema de teledetección. Fuente: V.Olaya.

1. Fuente de radiación (puede ser natural) o artificial. La primera de ellas es la forma más importante de teledetección, pues se deriva directamente de la luz solar, principal fuente de energía del plantea. El sol ilumina la superficie terrestre, que refleja esa energía en función del tipo de cubierta presente sobre ella. Este flujo reflejado se recoge por el sensor, que lo transmite posteriormente a las estaciones receptoras.
2. Objetos que interaccionan con la radiación o la emiten.
3. Atmósfera por la que se desplaza la radiación, tanto desde la fuente hasta el objeto como desde el objeto hasta el receptor. La atmósfera también interactúa con la radiación, introduciendo igualmente perturbaciones en ella.
4. Receptor que recoge la radiación una vez ésta ha sido perturbada o emitida por los objetos. El receptor va a generar como producto final una imagen en cuyas celdas se va a contener un valor que indica la intensidad de la radiación. Estos valores son valores enteros que indican el nivel de dicha radiación entre una escala definida (habitualmente valores entre 1 y 256), y se conocen dentro del ámbito de la teledetección como Niveles Digitales

La radiación electromagnética¶

La radiación electromagnética es una de las cuatro fuerzas fundamentales en la naturaleza, y deriva del campo electromagnético. Un cuerpo con carga, en reposo o que se mueve en el espacio, genera un campo electromagnético definido mediante dos vectores: el campo electrónico y el campo magnético. Así, se define la radiación electromagnética como una combinación de campos eléctricos y magnéticos oscilantes y perpendiculares entre sí que se propagan a través del espacio transportando energía de un lado a otro.

Esquema de onda electromagnética. Fuente: E. Chuvieco

Según la teoría ondulatoria, las ondas electromagnéticas se desplaza a la velocidad de la luz siguiendo un modelo armónico y continuo, y se pueden describir por dos elementos: longitud de onda (λ) y frecuencia (F). La primera hace referencia a la distancia entre dos tipos sucesivos de onda, mientras que la frecuencia designa el número de ciclos pasando por un punto fijo en una unidad de tiempo. Ambos elementos están inversamente relacionados: una mayor longitud de onda -y por tanto, una menor frecuencia- tiene asociada una menor energía de la radiación.

El espectro electromagnético¶

De lo mencionado hasta ahora se deduce que puede definirse cualquier tipo de energía radiante en función de su longitud de onda o frecuencia. Aunque la sucesión de valores de longitud de onda es continua, suelen establecerse una serie de bandas en donde la radiación electromagnética manifiesta un comportamiento similar. La organización de estas bandas de longitudes de onda o frecuencia se denomina espectro electromagnético, y comprende desde las longitudes de onda más cortas (rayos gamma, rayos X), hasta las Kilométricas.

En la siguiente imagen se muestra un esquema del espectro electromagnético y sus principales regiones de interés:

Espectro electromagnético y sus principales regiones de interés. Fuente: Wikipedia

La firma espectral¶

La firma espectral es el conjunto de reflectancias en las distintas longitudes de onda que presenta un objeto o material dado. Dicho de otro modo, la firma espectral es la ‘huella dactilar’ que caracteriza cada objeto, y que deriva directamente de sus propiedades y de cómo éstas influyen en la forma en que dicho objeto refleja la radiación incidente.

Es gracias a la información que proporcionan las firmas espectrales que es posible identificar en una imagen de satélite la naturaleza de los objetos.

Las firmas espectrales son a menudo representadas gráficamente en un eje de coordenadas, donde la longitud de onda está en el eje de las abscisas, y la reflectancia en el eje de las ordenadas.

Signaturas espectrales típicas de distintas cubiertas. Fuente: E. Chuvieco

El sensor OLI (Operational Land Imager)¶

A bordo del Landsat 8, se encuentra el sensor OLI, puesto en órbita en Febrero del 2013. Las bandas espectrales capturadas por este sensor son muy parecidas a las del sensor ETM+ (a bordo del Landsat 7), aunqué se han añadido dos de nuevas: un nuevo canal en azul visible (banda 1) diseñado específicamente para observar la calidad del agua en lagos someros y zonas costeras, así como para detectar aerosoles, y otro canal en el infrarrojo (banda 9) para determinar la presencia de nubes -fundamentalmente cirrus-.

Asimismo, junto cada escena del Landsat 8 se incluye también una banda de valoración de calidad (Quality Assurance band), que ofrece información respecto anomalías en la toma de datos por problemas en el instrumental o la presencia de elementos como nubes, agua y nieve.

El sensor TIRS (Thermal Infrarred Sensor)¶

El sensor TIRS capta información acerca la temperatura de la superficie terrestre a través de dos bandas del infrarroje térmico (banda 10 y banda 11). Permite distinguir entre la temperatura de la superficie terrestre y la temperatura atmosférica. Tienen una resolución de 100 mts.

Las bandas del satélite Landsat 8¶

Comparativa longitud de onda de las bandas del Landsat 8 y Landsat 7:

Combinación de bandas¶

Las bandas multiespectrales de una imagen de satélite pueden visualizarse mediante distintas composiciones en color. Basta por ello aplicar cada uno de los tres colores primarios (rojo, verde y azul) a una banda distinta de la imagen, seleccionada con el criterio y el orden que se estime más oportuno. El proceso permite visualizar simultáneamente, imágenes de distintas regiones del espectro, lo que facilita la delimitación visual de algunas cubiertas. La elección de las bandas para realizar la composición, y el orden de los colores destinados a cada una, dependen del sensor sobre el que se trabaje y de la aplicación última del proyecto.

Encontraréis algunas combinaciones posibles para Landsat 8 en: http://landsat.usgs.gov/L8_band_combos.php

Split-window algorithm¶

Para el cálculo de la temperatura de la superficie, se utilizará el algoritmo Split-Window (SW). Se trata del algoritmo mas utilizado para el cálculo de LST debido a su simplicidad y robustez.

Este algoritmo se basa en el hecho que la radiación absorbida por la atmósfera es proporcional a la diferencia de brillo entre las mediciones simultáneas en dos longitudes de onda diferentes, correspondientes a las dos bandas del sensor TIRS.

LST = TB₁₀+ C₁(TB₁₀-TB₁₁) + C₂(TB₁₀-TB₁₁)2+ C₀+(C₃+C₄W) (1-ε) + (C₅+C₆W) ∆ ε

donde:

LST - Land Surface Temperature (Kelvin)

C₀ to C₆- Valores del coeficiente para SW (Skokovic et al, 2014; Sobrino et al, 1996; Shaouhua Zhao et al, 2009)

TB₁₀y TB₁₁- temperatura de brillo de la banda 10 y 11 (K)

∆ - valor medio LSE (Land Surface Emissivity) de las bandas del TIR

W - contenido de vapor de agua en la atmósfera

∆ ε - Diferencia en LSE

Convertir los ND a valores de radiancia (TOA - Top of athmosphere)¶

A partir de los datos medidos por el sensor, puede obtenerse la energía reflejada, ya que la radiancia espectral medida por éste es consecuencia de la reflexión de la radiación electromagnética en las cubiertas. Esta reflexión se codifica con un valor numérico, que se denomina ND, de acuerdo a los coeficientes de calibración específicos para cada sensor. Dada que estos coeficientes son conocidos, puede realizarse el proceso inverso, obteniendo así los valores de radiancia espectral detectado por el sensor a partir de los ND.

Para aplicar la fórmula del SW algorithm, es imprescindible conocer estos valores de radiancia. Las bandas del sensor OLI y TIRS pueden convertirse a valores de radiancia espectral del siguiente modo:

Lλ = MLQcal + AL

donde:

Lλ = TOA spectral radiance (Watts/( m2 * srad * μm))

ML = Band-specific multiplicative rescaling factor from the metadata (RADIANCE_MULT_BAND_x, where x is the band number)

AL = Band-specific additive rescaling factor from the metadata (RADIANCE_ADD_BAND_x, where x is the band number)

Qcal = Quantized and calibrated standard product pixel values (DN)

Conversión a valores de reflectividad aparente y temperatura de Brillo a la altura de la atmósfera (TOA)¶

A partir de los valores de radiancia, utilizando las bandas térmicas fácilmente se puede derivar la temperatura a la altura del sensor. Las bandas del TIRS se pueden convertir a temperatura de brillo utilizando las constantes que figuran en el archivo de metadatos. De este modo:

donde:

T = At-satellite brightness temperature (K)

Lλ = TOA spectral radiance (Watts/( m2 * srad * μm))

K1 = Band-specific thermal conversion constant from the metadata (K1_CONSTANT_BAND_x, where x is the band number, 10 or 11)

K2 = Band-specific thermal conversion constant from the metadata (K2_CONSTANT_BAND_x, where x is the band number, 10 or 11)

Abrid el archivo de metadatos de la imagen (.met) con un editor de texto, y comprobad que figuran todos los parámetros para el cálculo de la reflectividad: date, sun_elevation, product_date, gain.

Corrección atmosférica¶

Este método, no obstante, no tiene en consideración las influencias atmosféricas y asume un terreno plano y una observación vertical, lo cual supone una simplificación de la realidad. De ahí que debamos hablar de reflectividad aparente, pues el valor calculado representa sólo la reflectividad captada por el sensor, pero no la realmente medible en la superficie.

En cuanto a la observación vertical, es asumible para la mayor parte de los sensores de interés ambiental (Landsat, IRS, MOS-MESSR...), pero debe considerarse cuando la adquisición no es vertical, como ocurre con el SPOT-HRV.

Por otro lado, la influencia atmosférica no afecta por igual a los dos componentes del cálculo de la reflectividad (energía reflejada e incidente), ya que el espesor de la atmósfera que atraviesan es distinto. Además, hay que considerar que a la radiancia solar directa hay que añadir la difusa, procedente de otros objetos vecinos. Por ello, la radiancia que recibe el satélite no es la misma que la que sale del suelo, que es la que interesa, dado que pretendemos medir la reflectividad de la cubierta, no la influida por la atmósfera.

Existen distintos procedimientos para abordar las influencias atmosféricas, y determinar la radiancia del suelo, no la que recibe el satélite. Desde GRASS, disponemos del módulo i.landsat.toar que permite aplicar el método de corrección DOS (Dark Object Substraction), y i.atcorr, para aplicar el método 6S (Second Simulation of Satellite Signal in the Solar Spectrum).

Por su sencillez, utilizaremos el método DOS, que consiste en suponer que el mínimo valor de ND debe corresponder a las zonas oscuras presentes en la imagen y que en ausencia de efecto atmosférico ese valor debería ser cero. Conforme a esta suposición, se deduce que las diferencias entre el cero y los valores mínimos de los histogramas de las distintas bandas afectadas se deben al incremento de la radiancia absorbida por el sensor como consecuencia de la radiación difusa de la atmósfera. El procedimiento de corrección consiste en restar de todos los ND de cada banda, el ND mínimo de ella, de modo que se haga coincidir con el cero del origen del histograma.

Existen algunas variantes del método de corrección atmosférica DOS (DOS1, DOS2, DOS3 y DOS4)

http://gis.stackexchange.com/questions/40531/how-to-determine-aerosol-model-value-for-i-atcorr-in-grass

Cálculo de valores de emisividad de la superficie terrestre¶

Otro de los parámetros del SW algorithm que es necesario obtener, es el de la emissividad de la superficie terrestre. Para entender que es la emisividad, hay que hacer referencia previamente a la emitancia.

A partir de Índice de Vegetación de Diferencia Normalizada (NDVI) es posible derivar el valor de emisividad de las cubiertas del suelo.

El índice de Vegetación de Diferencia Normalizada

El Índice de Vegetación de Diferencia Normalizada, también conocido como NDVI por sus siglas en inglés, es uno de los índices más importantes y más ampliamente utilizados.

La fórmula para la obtención de este índice es la siguiente:

NDVI = (IR – R) / (IR + R)

donde R e IR son las reflectancias correspondientes al rojo (Banda 4) e infrarrojo (Banda 5) respectivamente.

Los valores resultantes de este índice se encuentran dentro del intervalo (-1,1), indicando los valores altos la presencia de vegetación.

(NDVI < 0), correspondientes a agua o cubiertas artificiales

(0 < NDVI < 0,2), correspondientes a suelo desnudo o vegetación muerta

(0,2 < NDVI < 0,4), correspondientes a vegetación dispersa o poco vigorosa

(0,4 < NDVI < 0,6), correspondientes a vegetación abundante y vigorosa

(NDVI > 0,6), correspondientes a vegetación muy densa y vigorosa,

Fractional Vegetation Cover (FVC) El FVC es un índice que permite estimar la fracción de superficie ocupada por vegetación, y se obtiene a partir del NDVI. Obtener el FVC es necesario para hallar los valores de LSE.

Land Surface Emissivity

Finalmente, podréis obtener los valores correspondientes a LSE para las bandas del sensor TIRS (Band 10 y Band 11), teniendo en consideración que:

donde:

Fuente: Skokovic et al, 2014; Sobrino et al, 1996; Shaouhua Zhao et al, 2009.

Aplicación del Split-Window Algorithm¶

Recordad que el SW Algorithm que vamos a aplicar, es:

LST = TB₁₀+ C₁(TB₁₀-TB₁₁) + C₂(TB₁₀-TB₁₁)2+ C₀+(C₃+C₄W) (1-ε) + (C₅+C₆W) ∆ ε

En este momento, ya únicamente nos faltan los parámetros correspondientes a:

∆ - valor medio LSE (Land Surface Emissivity) de las bandas del TIR = (LSE₁₀+ LSE₁₁) / 2

∆ ε - Diferencia en LSE = LSE₁₀- LSE₁₁

ref_B10+1.378*( ref_B10 - ref_B11)+0.183*( ref_B10 - ref_B11 ) ^2-0.268+ (54.300-2.238*0.013 )*(1- LSE_valorMedio)+(-129.20+16.40*0.013)* ( LSE_B10 - LSE_B11

Obtener valores en ºC¶

Los valores de temperatura que habéis obtenido, están en grados Kelvin. Para obtener los valores en ºC, recordad que 0 Kelvin = -273.15 Centígrados

Descargar código de los ejercicios¶

Debéis de ejecutar el script ‘obtener_datos_taller.sh’ que tenéis en el escritorio.

Pygrass Library (grass.script)¶

Pygrass es una API orientada al objeto para trabajar con GRASS, que permite a los programadores trabajar con GRASS de una forma mas fácil y natural con un lenguaje de script.

Abrir grass70

$ grass70

En la consola de GRASS

help(grass)

Nos muestra la ayuda y los módulos que dispone. El código de estos de encuentra en /usr/lib/grass70/etc/python/grass/script/

Al ejecutar en la consola de GRASS, nos muestra la del módulo ‘grass.utils’.

help(grass.utils)

Al ejecutar en la consola de GRASS, nos muestra la del módulo ‘grass.utils’.

help(grass.utils)

La consola nos permite escribir código y rutinas que ejecuten procesos y interactuar con la interfaz gráfica.

En la consola de GRASS. Vemos como carga la capa LST en el árbol de capas.

AddLayer('LST')

En la consola de GRASS. Obtenemos la lista de imágenes que tenemos.

grass.list_strings('raster')

En la consola de GRASS. Obtenemos información acerca de la imagen ‘LST’.

grass.raster.raster_info('LST')

Mas información en:

http://grass.osgeo.org/grass71/manuals/libpython/pygrass_index.html

Además de la consola de python de GRASS también podemos trabajar con GRASS desde scripts de python.

$ python nuestro_script.py

Y que este ejecute cógigo que llame a los módulos de GRASS.

Variables de entorno¶

Antes de empezar a programar con python y GRASS desde scripts hay que definir y inicializar las variables de entorno y variables de GRASS así nuestros scripts de python puedan funcionar de forma correcta.

#!/bin/bash

#http://grasswiki.osgeo.org/wiki/GRASS_and_Python
#http://grass.osgeo.org/grass70/manuals/variables.html

export GISBASE=/usr/lib/grass70
export GRASS_VERSION="7.0.0"

GISRC="$HOME/.grassrc7"

export GISRC
export GRASS_PYTHON=python

export PATH="$GISBASE/bin:$GISBASE/scripts:$PATH"
export LD_LIBRARY_PATH="$GISBASE/lib"
export GRASS_LD_LIBRARY_PATH="$LD_LIBRARY_PATH"
export PYTHONPATH="$GISBASE/etc/python:$PYTHONPATH"


#Generamos el GISRCRC
MYGISDBASE=$HOME/Documentos/grassdata
MYLOC=catalunya
MYMAPSET=taller
echo "GISDBASE: $MYGISDBASE" > "$GISRC"
echo "LOCATION_NAME: $MYLOC" >> "$GISRC"
echo "MAPSET: $MYMAPSET" >> "$GISRC"

http://grass.osgeo.org/grass70/manuals/variables.html

Ejecutar script que inicializa nuestras variables de entorno y de GRASS.

$ source grass-env.sh

Ejemplo que muestra las variables de entorno de GRASS.

"""Muestra las variables de entorno de GRASS"""

import grass.script as grass

env = grass.gisenv()

print env 

Para ejecutar el script

$ python ejemplo_01_entorno.py

Si os fijáis en el código, las variables de GRASS se almacenan en el archivo /home/user/.grassrc7

GISDBASE: /home/user/Documentos/grassdata
LOCATION_NAME: catalunya
MAPSET: taller

python¶

En la consola

$ python

Para importar la librería de GRASS

import grass.script as grass

Ejemplo que muestra los límites del location.

import os

import grass.script as grass

#Print region extent
r = grass.read_command("g.region", flags='p' )
print r

Para ejecutar el script

$ python ejemplo02_location.py

Sintaxi¶

http://grass.osgeo.org/grass70/manuals/libpython/script_intro.html#syntax

Consola de bash o CMD de windows

r.info map=NDVI -g --verbose

En python

grass.read_command('r.info', map='NDVI', flags='g', verbose=True)

>>> grass.run_command('g.remove', type='raster', name='raster1', flags='f')
0

>>> grass.read_command('g.region', flags='p')
'projection: 1 (UTM)\nzone:       31\ndatum:      wgs84\nellipsoid:  wgs84\nnorth:      4741215\nsouth:      4503885\nwest:       370485\neast:       603915\nnsres:      30\newres:      30\nrows:       7911\ncols:       7781\ncells:      61555491\n'

>>> grass.parse_command('g.region', flags='p')
{'ewres:      30': None, 'ellipsoid:  wgs84': None, 'nsres:      30': None, 'rows:       7911': None, 'north:      4741215': None, 'cells:      61555491': None, 'projection: 1 (UTM)': None, 'zone:       31': None, 'cols:       7781': None, 'west:       370485': None, 'datum:      wgs84': None, 'east:       603915': None, 'south:      4503885': None}

Ejercicios¶

En la carpeta ‘/home/user/Documentos/taller_grass/workspace/ejercicios’ tenéis los archivos que debéis de completar.

• ejercicios_00.py
• ejercicios_01.py
• ejercicios_02.py
• ejercicios_03.py
• ejercicios_04.py
• ejercicios_05.py
• ejercicios_06.py
• ejercicios_07.py
• ejercicios_08.py

Los ejercicios son los mismos que habéis realizado con la interfaz gráfica de GRASS pero con una imagen distinta. En el ejercicio 08 obtendremos datos estadísticos de la imagen ‘LST’ y ‘LST2’ y los cruzaremos con la capa de ‘barrios’ de Girona.

Si todo funciona correctamente podemos ver los datos en un gráfico al ejecutar en la consola de bash.

$ python ejercicios_limpiar.py

Ejemplo que os indica que debéis de implementar el código del ejercicio.

##EJERCICIO## ##EJECUTAR_COMANDO## ##OBTENER TEMPERATURA DE BRILLO## ##PARÁMETROS##

Antes de empezar debés de limpiar todo lo que habéis hecho menos la capa ‘LST’ para ello tenéis el archivo ‘ejercicios_limpiar.py’

En la consola de bash

$ python ejercicios_limpiar.py