Download Procesamiento de Datos de NDVI para la obtención de

Procesamiento de Datos de NDVI para la
obtención de Precipitación de Lluvias en
la Cuenca del Altiplano
Yarlequé, C. P 1, Posadas, D. A.2 y Quiroz,
R.2
Centro Internacional de la Papa2,
División de Manejo de Recursos Naturales
PO Box 1558, Lima 12.
Universidad Nacional del Callao,
Septiembre, 2004
Resumen
El pronóstico de Precipitación de lluvias, es
un tema muy investigado en el Perú (y el
Mundo). El presente trabajo nos muestra el
análisis de Precipitación por intermedio de
datos de NDVI (índice de vegetación),
siendo ambos eventos periódicos y
proporcionales. Se presenta una aplicación,
con datos reales de precipitación y NDVI,
ambos en la Cuenca del Altiplano Peruano,
Puno, durante un dominio de 5 años de
tomas de datos, en distintas escalas
(muestras diarias y cada 10 días,
respectivamente).
La
técnica
de
Transformadas de Wavelets es usada al
aplicarse en los datos de NDVI, para
obtener la misma resolución temporal
(escala) que la de precipitación, y la técnica
de la Transformada de Fourier aplicada
para obtener, la periodicidad de los datos.
Ambas técnicas en conjunto nos producen
una función similar a la que describe la
precipitación de lluvias.
I.
INTRODUCCIÓN
La cuantificación de la precipitación pluvial es
de suma importancia para estimar la
disponibilidad de agua para uso domestico,
agricultura, generación de energía y otros.
Para ello se utilizan estaciones meteorológicas,
que por su costo no es posible tener el numero
mínimo necesario para cubrir adecuadamente
la variación espacial en una región o país.
Debido a que el crecimiento de las plantas en
tomas de secano responde directamente a la
cantidad y distribución de la lluvia, la
evaluación de este fenómeno puede ser usado
para cuantificar la precipitación.
Los Índices de Vegetación v.s. el Índice
Vegetación normalizado de la diferencia
(NDVI) pueden ser utilizados para dar
seguimiento a la variación en “Verdosidad”
producido por el crecimiento de las plantas y
por ende cuantificar la precipitación. El
procesamiento del NDVI para cuantificar la
proporcionalidad y periodicidad de la
precipitación es descrito en este articulo.
Palabras Clave:
Precipitación; Índice de
Vegetación; Transformadas; Wavelets;
Fourier.
1
Bachiller - CIP cyg_fisic@yahoo.com ,
cyg_pollux@hotmail.com
2
CIP. a.posadas@cgiar.org , r.quiroz@cgiar.org
Figura 1. Técnica de adquisición de datos
(digitales); Remote Sensing.
Toda transformada (matemática) nos
produce un cambio de dominio (variables
independientes) en la función analizada,
para
poder
percibir
las
diversas
características de las funciones en otro
dominio (por ejemplo; una función que da
valores de temperatura, se puede aplicar la
transformada de Fourier para obtener la
frecuencia en que ocurre la intensidad de
cada temperatura), pero la función sigue
siendo la misma. En el presente trabajo se
utilizó las Transformadas de Fourier y
Wavelets, para el procesamiento de datos
de NDVI, decadales (muestreo cada 10
días), para obtener una función de igual
clase
como
la
de
Precipitación.
Posteriormente, comparamos el resultado
dado con datos reales de precipitación de
lluvia.[8]
II.
PRECIPITACIÓN, NDVI Y
TRANSFORMADAS
Los campos de cultivos, bosques,
pastizales, etc. nos brindan por intermedio
de Percepción Remota, datos para la
medición de Biomasa (Vegetación, etc.).
Estos campos pueden ser interpretados
geoespacialmente; es decir la imagen que es
tomada por los sensores de los satélites u
otros, lo podemos obtener en un formato tal
que nos brinde una composición matricial
del mismo. Así es posible analizar los
campos de biomasa (imágenes en general, ó
señales en 2D) utilizando herramientas
físico-matemáticas.
La vegetación saludable absorbe mas luz
roja, que es utilizada en la fotosíntesis y
refleja radiación infrarroja (IR), la
vegetación estresada y el suelo reflejan la
luz roja y la radiación IR en la misma
cantidad. Así, se tiene que la magnitud
relativa de luz roja y radiación IR censada
desde un cultivo provee un indicador del
vigor del mismo. El cultivo en las partes del
campo con suelos pobres tiende a madurar
tempranamente. Ambos, el bajo rendimiento y
la madurez temprana pueden ser detectados
convirtiendo la fotografía aérea en una imagen
digital (datos discretos) y calculando una
cantidad llamada “Índice de Diferencia
Normalizada de Vegetación (NDVI)”. Esta es
una proporción de valores de luz roja é IR.
(IR − R )
NDVI =
(1)
(IR + R )
El cálculo de NDVI para un píxel dado
siempre resulta en un numero dentro del rango
de –1 a +1. Hojas que no son verdes dan
valores cercanos a cero (0); El cero indica que
no existe vegetación y valores cercanos a +1
(0.8~0.9) indican la más alta densidad posible
de hojas verdes.
Las imágenes de un campo dado en radiación
IR y rojo; son procesadas según (1), para
poder obtener una con valores de NDVI; ya
que estos valores muestran mejor información
respecto a la existencia de biomasa.
La Transformada de Fourier (TF), es la
herramienta que nos brinda el análisis de una
cierta señal con un comportamiento periódico
(frecuencias, periodos, amplitud, ...). Es decir;
nos brinda información adicional de la señal
respecto a la frecuencia de repeticiones de
eventos que describe la misma. En un conjunto
de datos de Precipitación, nos puede mostrar
directamente su frecuencia característica para
poder predecir comportamientos (estándar)
posteriores. Además, brinda una buena técnica
de filtrado; exonerando frecuencias no
deseadas (altas frecuencias, ruidos, etc.),
alisando así la imagen.
∞
F(w) =
∫ f(t) e
-iwt
dt
(2)
−∞
Podemos mostrar la ecuación discreta
equivalente, si tomamos el cambio de variable:
w=2πn/T, donde; n=1, 2, 3...N. y T=periodo
(para el caso discreto, T=N).
1 N −1
-i 2πnt N
F (n) = ∑ f (t ) e
= an + ibn (3)
N t =0
2
La Transformada Wavelet (TW), es un
operador que reescala una señal de mayor a
menor píxeles, ó viceversa haciendo variar
su resolución espacial (tamaño). Al
reescalarse una señal se pierde o se gana
detalles de ella (bajas y altas frecuencias,
respectivamente), debido a la compresión o
elongación (espacial) de la señal. Este
reescalamiento es diferente a otros, por la
caracterización especifica que se da a la
señal, en cada proceso de escalamiento. Al
igual que la TF, la TW puede realizar un
proceso de filtrado al reproducir la señal en
otra escala con la ayuda de funciones bases,
denominados Wavelets (o mas conocidos
como Wavelet Madre). Este Wavelet
(escogido) convoluciona con la señal,
produciendo una nueva señal reescalada y
con mayor o menor detalles.
Wf (λ , t ) =
∞
∫ f (u)ψ λ
,t
(u )du
(4)
−∞
∞
=
∫ f (u)
−∞
1
λ
ψ(
u −t
λ
)du
(5)
ψ λ ,t (u ) : Son llamados Wavelet. λ, es el
valor de la variable escala, relacionada a las
mediadas de la variable t (tiempo). Son
equivalentes a filtros de distintas clases.
La serie de Fourier, es utilizada en los
procesos de reconstrucción de señales:
N
a
f (t ) = 0 + ∑ (a n cos nt + bn sen nt ) (6)
2 n =1
III.
197 imágenes de agregados decadales de
NDVI, provenientes de los satélites SPOT (4 y
5) - vegetation (VGT1 y VGT2),
correspondientes al mismo periodo de los
datos de precipitación, fue utilizada en el
análisis. Las imágenes fueron corregidas
geométricamente y se extrajo la información
contenida en las siguientes coordenadas:
DATOS RECOPILADOS
Los datos procesados son de dos clases:
precipitación de lluvia y NDVI. Los datos
diarios de precipitación fueron tomados por
la estación climatológica, Mazo Cruz,
perteneciente al SENAMHI, ubicado en las
coordenadas; Lat. 16°44’ ”S”, Long. 69°42’
“W”, Alt. 4100 msnm. Se proceso el
periodo cubierto desde el 01/01/98 hasta el
31/12/02. Una base de datos conteniendo
Figura 2. Datos de NDVI (área) tomados
desde el 11/04/98 hasta el 11/09/03;
Superpuestos con los datos de Precipitación
de lluvia, tomados en la estación Mazo Cruz
(punto), desde 1/01/98 hasta el 31/12/02.
Las coordenadas de la zona geográfica de las
197 imágenes de NDVI son:
X0
=70°14’5.78’’; Y0 = 16°3’55.36’’; Yf =
17°26’42.10’’; Xf = 68°51’19.04’’; tomados
desde el 11/04/98 al 11/09/03.
IV.
ANÁLISIS
Los datos comunes seleccionados según se
muestra en la figura 3. Estos fueron
procesados en programas diseñados en los
Software MatLab (TW) e IDL (TF); los cuales
fueron estructurados para la utilización en
otros tipos de datos.
La figura 4. muestra los datos de
Precipitación, con sus respectivas equivalentes
en distintas escalas (abajo), es decir ha sido
reescalado a una resolución espacial (escala)
3
menor (menos píxeles), aplicando filtro de
baja y alta frecuencias (forma de la señal y
ruido respectivamente). A continuación se
aplica la Transformada de Fourier, para el
análisis periódico de la señal. El espectro de
frecuencia (figura 5.), muestra las
principales frecuencias que componen a la
Los detalles en los datos de precipitación con
alta frecuencia son mas caóticos debido a la
alta variación en los datos diarios; mientras los
picos de las señales de baja frecuencia, indican
la forma general de esta señal en distintas
escalas.
Con la técnica de Fourier fue posible obtener:
Dominio de Frecuencias que contiene
la señal de Precipitación y NDVI. Siendo muy
similares entre si.
La amplitud característica (n=0), nos
indica el valor medio de la señal y la
utilizamos para poder obtener el valor
proporcional de cambio de escala (relación de
escala), de las señales, en las amplitudes del
espectro real (ver (3)):
Flluvia (0)
1.52167
=
= 57.13068 (7)
FNDVI (0) 0.0266349
Como estamos trabajando con datos muy
diferentes en proporción es necesario
considerar
los
mayores
dígitos
no
significativos (5 dígitos). Con esta relación de
escala es posible luego multiplicarlos a los
valores de NDVI para obtener valores
correspondientes a los de Precipitación.
señal.
Figura 3. Interceptación en el dominio
temporal de los datos de Precipitación y
NDVI.
Igualmente los datos de NDVI, son
reescalados pero en una resolución espacial
mayor (mas píxeles), debido que los
detalles característicos (forma) de la señal
se perderían si bajamos de escala (figura 6
y 7).
“Se aprecia el comportamiento periódico
de las señales, en desfase, pero
proporcionales”.
La deducción del valor de longitud de
onda de cada señal, a partir de la señal
aplicada TF, es factible, con un índice de error
de 1pixel por periodo (píxel≡1día). Este error
es notado cuando se aprecian los datos
reconstruidos.
Los periodos son respectivamente:
tlluvia=364 datos (dato≡1día) y tNDVI=36 datos
(dato≡10dias). Se deduce así el desfase
existente en las señales (Ver figuras 4 y 6).
Con un valor aproximado de 175 píxeles (o
datos) de Precipitación. Estos datos son
obtenidos en un procesador de datos común;
con la teoría de ondas convencionales, gracias
al comportamiento periódico de la Señal de
4
Precipitación que se muestran en las figuras
5 y 7.
Figura 4. Datos de Precipitación de lluvia,
en la Estación Mazo Cruz, diarios. Las dos
Imágenes inferiores son aplicaciones de la
TW a una menor resolución espacial a
distintas frecuencias (baja y alta).
Figura 5. Datos de precipitación. Señal
reconstruida con la TF, con los 6 primeras
frecuencias (mas representativas)[2]. Abajo
el espectro de Frecuencias de la señal
respecto al valor de frecuencias (n; ver (3)).
Figura 6. Datos diarios de NDVI, ubicados en
las coordenadas de la Estación Mazo Cruz,
diarios. Las dos Imágenes inferiores son
aplicaciones de la TW a una menor resolución
espacial a distintas frecuencias (baja y alta).
Figura 7. Datos de NDVI. Señal reconstruida
con la TF, con los 6 primeras frecuencias
(mas representativas)[2]. Abajo el espectro de
Frecuencias (n; ver (3)).
V.
5
RESULTADOS
RECONSTRUCCIÓN.
La figura 8. muestra una reconstrucción sin
aplicar ninguna técnica de filtrado. Se
aprecia como los datos pierden detalle y no
nos indican los valores internos de la señal
(entre cada 10 píxeles). Además, hay que
tomar en cuenta el error existente en la
exoneración de datos (días) que no se
incluyen en los datos de NDVI (197
imágenes), debido que entre los datos,
existen, 10, 11, ó 9 días de diferencia en la
recopilación de los mismos, esto debido a
que los meses y años (regulares y
Bisiestos), contienen diferencias de días
entre ellos.
espacial, desde los datos mostrados en la
figura 8, y se obtuvo una señal en baja
frecuencia de la misma forma de la figura 6) y
de alta frecuencia (señal mostrada en la parte
inferior-derecha de la figura 4, baja resolución
espacial) obtenido al aplicar TW a la señal de
NDVI y Precipitación. Adicionalmente, al
realizar la aplicación de TW (tanto en
reescalamiento a mayor escala o menor
escala), la técnica nos permite aplicar un filtro
(Wavelet) por cada reescalamiento; así que se
utilizaron dos tipos de Wavelets: ‘Haar’ (baja
frecuencia) y ‘Symlets’ (alta frecuencia).
Se aplica a estos datos el procesamiento,
debido a que estos nuevos datos de NDVI,
tiene casi la misma resolución espacial
(píxeles) que los de Precipitación. Y se
procede a obtener los mismos datos y
señales, como se muestra en la figura 6.
Figura 9. Reconstrucción de la Señal de
NDVI, con características de Precipitación.
Aplicando Wavelet y Fourier, se caracteriza la
señal de NDVI, hasta llegar a una del tipo
Precipitación.
Figura 8. Simple reconstrucción espacial.
Cada 10 datos corresponde a un valor
respectivo de cada uno de los datos de
NDVI. (Reconstrucción no detallada).
La imagen a continuación muestra la
reconstrucción completa de la señal de
NDVI, con las técnicas de TW, y TF. Se
utilizó el ruido (alta frecuencia) de la señal
de Precipitación (figura 4) a una resolución
distinta
para
poder
producir
el
reescalamiento, filtrando datos con el ruido
(alta frecuencia) que contiene la
Precipitación. Es decir, se juntaron dos
funciones (señales): de baja frecuencia
(señal de NDVI, en baja resolución
Es posible la realización del mismo, mediante
el uso solamente de ruido de NDVI; ya que en
la figura 6 (en alta frecuencia), se observa la
gran similitud de ruido existente en la señal de
NDVI, con la señal mostrada en la figura 4 (en
alta frecuencia) de los datos de Precipitación.
COMPARACIÓN.
Aquí se muestran el resultado final. aplicando
TW y TF; los datos superiores son los de
NDVI reconstruidos con la señal frecuencia
baja (figura 6) y reescalada de la señal de
NDVI.
Para referencia alguna se indica el valor
cuadrático de la media residual, s2=19.9; y su
6
respectivo valor de media residual: s=4.46.
Por lo dicho en la introducción de la parte
de RECONSTRUCCIÓN, El error es
perturbado principalmente por la no
concordancia de datos al realizar el
reescalamiento (reconstrucción) de los
datos de NDVI.
La reproducción de las señales de
Precipitación, es realizada mediante un
proceso de reescalamiento y filtrado que nos
ayuda a seleccionar y clasificar los datos mas
útiles en el análisis agronómico a realizar; así
es útil este procesamiento para analizar la
variedad de eventos que suceden en el tiempo
que muestran los datos (mayor lluvia,
frecuencia de eventos, etc.. ).
Las señales tienen un comportamiento
periódico demostrado con la TF, es posible
realizar una extrapolación de los datos para
hallar valores posteriores a los tiempos
tomados, y tratar de producir una equivalente
señal de Precipitación de lluvia; es decir,
realizar pronósticos de lluvia en tiempos
indeterminados (futuro o pasado).
Figura 10. Comparación de data en fase.
La señal superior son los datos
segmentados de la Reconstrucción de la
señal de NDVI al tipo Precipitación. La
figura inferior datos reales (en fase con los
de NDVI), de la precipitación de lluvia.
Cabe recalcar que en el presente trabajo
únicamente se realiza el procesamiento de
datos sin tener en cuenta detalles técnicos del
sector agrícola (acotamiento de datos de
NDVI por saturación de lluvias ó sequías,
proporcionalidad de los datos de NDVI con
los de Precipitación, etc...), se vera mejorado
indiscutiblemente añadiendo detalles que
solamente pueden ser indicados por expertos
en dicho sector.
VI.
CONCLUSIONES
El análisis geoespacial confirma la
proporcionalidad de los datos tomados de
NDVI con los de Precipitación de lluvia en
una misma ubicación geográfica.
[1] RICHARDS JOHN A., (1995), Remote Sensing
Digital Image Análisis, Springer-Verlag, 2° edition.
Campbell, Australia.
[2]W. W. IMMERZEL, R. A. QUIROZ, S. M. DE
JONG; “Understanding complex spatiotemporal
weather patterns and land use interaction in the Tibetan
Autonomous Region using harmonic análisis of SPOT
VGT-S10 NDVI time series”. (2004)
[3] MARIAN PRUTSCHER (1998), Series de Fourier,
http://www.e-technik.uniulm.de/world/lehre/basic_mathematics/fourier/node2.ph
p3
[4] GONZÁLES RAFAEL C. AND RICHARD E.
WOODS, (1992), Digital Image Processing, Editorial
Addison-Wesley Publishing Company INC,Firts
edition.
[5] POLIKAR ROBI, (1996), “The Wavelet Tutorial”,
329 Durham Computation Center Iowa State
University.
http://users.rowan.edu/~polikar/WAVELETS/WTtutori
al.html
[6] EFI FOUFOULA-GEORGIOU AND PRAVEEN
KUMAR, (1994), Wavelets in Geophysics, editorial
Academic Press inc, Primera edición.
[7] Jhon C. Russ, THE IMAGE PROCESSING North
Carolina Handbook, Biblioteca del CIP. Editorial CRC
Press LLC, third Edition, año 1999.
[8] Plant R. Munk D., (2001), Application of Remote
Sensing to Strategic Questions in Cotton Management
and Research, the Journal of Cotton Science 5:30-41.
7