Cálculo de la reflectancia aparente en imágenes de satélite con QGIS

En el artículo anterior señalé un método para determinar la reflectancia aparente en imágenes de satélite con R. A pesar de la potencia del lenguaje, que lo hace un procedimiento expedito, puede no parecer atractivo a muchos porque se maneja a través de líneas de comando; tanto en Linux como en Windows. Por esa razón, me dispuse a encontrar una vía alternativa usando la calculadora raster de QGIS; donde los valores de los números digitales (ND) obtenidos en R los empleé para calibrar el procedimiento en QGIS.

En un artículo precedente, referí que la determinación de la radiancia en el tope de la atmósfera (L) es el paso previo para la determinación de la reflectancia aparente en imágenes de satélite. La fórmula para el cálculo es la siguiente:

donde ρAS es la reflectancia aparente, d es la distancia Tierra-Sol (expresada en unidades astronómicas), Esun es la constante atmosférica extra-solar (Wm^-2µm^-1) y θz es el ángulo cenital (que se calcula a partir del ángulo de elevación del sol, θs, mediante θz = 90 – θs).

El valor de d, para la fecha de adquisición de la imagen (“1988-03-26″), lo determiné con R (comando landsat::ESdist(“1988-03-26”)). Se obtuvo un valor de 0.9971208. El Esun, constante atmosférica extra-solar (Wm^-2µm^-1), la obtuve de una tabla de Chander et al. (2009).

        Landsat 5     Landsat 7
Band 1    1983           1997
Band 2    1795           1812
Band 3    1539           1533
Band 4    1028           1039
Band 5     219.8          230.8
Band 7      83.49          84.90
Band 8      NA           1362

En mi caso, como es una banda 1 (Landsat 5), el valor de Esun = 1983. Finalmente, de la metadata se extrae que el ángulo de elevación solar es de 54.89 grados.

Todos estos valores combinados en la expresión de arriba, señalan que la constante:

K = π . d² / Esun . Cos(90 – θs) = 0.00192549860081712

es el número que se empleará para multiplicar los valores de L para así obtener los de ρAS, es decir, se multiplica dicha banda (b1_esc_rad.tif) por 0.00192549860081712 con la calculadora raster de QGIS.

Como ya disponemos de los valores de L, basta multiplicarlos por 0.00192549860081712 para verificar si éstos son iguales a los que se obtuvieron con R. Efectivamente, la transformación pixel a pixel fue hecha de la manera esperada. A modo de ejemplo, para el primer pixel:

1.492156982421875 * 0.00192549860081712 = 0.00287314618185282

que concuerda, a partir de la décima cifra decimal (precisión que se usó para estimar la constante con LibreOffice), con el valor de 0.0028731462080031633377 obtenido con R para el primer pixel (e idéntico al que se encuentra en la imagen ASCII de dicha banda).

Por tanto, con QGIS, se carga la banda (b1_esc_rad.tif), se selecciona Raster -> Calculadora Raster, se introduce la fórmula, se coloca un nombre adecuado (b1_esc_ar.tif en este caso) y click para aceptar.