Omitir los comandos de cinta
Saltar al contenido principal
Inicio de sesión

Saber más

Administrar permisosAdministrar permisos
Caso 2: Imágenes Sentinel2 - Seguimiento del ciclo de cultivos en una serie de parcelas (post 1/2)
07/04/2020 - Producto, Formación, Noticia

​Primer post en el que evaluaremos imágenes NDVI del año agrícola 2018-2019 (1-septiembre-2018 a 31-agosto-2019); para hacer el seguimiento de cultivos en parcelas del comunal del municipio de Cortes. Sirva de ejemplo para casos similares a cargo de Cooperativas Agrarias, propietarios de múltiples parcelas dispersas, etc.

Destacar que el sistema agrario de este municipio está principalmente basado en el regadío de los canales de Tauste e Imperial de Aragón, con agua del río Ebro; teniendo, por tanto, gran importancia la agricultura en general, como la hortofruticultura en particular.

Empezamos con una sesión nueva de QGIS. Para tener una referencia geográfica, incorporaremos tanto el Mapa Base de IDENA como la Ortofotografía de 2018 (más reciente de verano, si más adelante ya está la de 2020, mejor) por servicio teselado (WMTS) dando de alta el servicio:

 

http://idena.navarra.es/ogc/wmts/1.0.0/WMTSCapabilities.xml

De esta forma, fijamos el SRC del proyecto en EPSG:25830, que es el oficial. QGIS hará una transformación al vuelo excelente con información que tenga otros SRC.

Descargaremos las imágenes desde:

https://filescartografia.navarra.es/3_ORTOFOTOGRAFIA/3_00_SATELITES/3_00_1_SENTINEL2/

En el ejemplo, descargaremos e incorporaremos a QGIS las imágenes NDVI disponibles de ese periodo:

Como el algoritmo principal que utilizaremos, presenta las imágenes NDVI conforme se encuentran en el Panel de capas, es conveniente ordenarlas de forma cronológica.

Ahora descargamos e incorporamos un Excel con las parcelas del comunal:

Cortes_ParcelaRustica_Comunal

Contiene un atributo REFCAT con la agregación de municipio, polígono y parcela, otros datos catastrales y un DNI imaginario del adjudicatario de la parcela.

Y desde la página de descargas del Servicio de Riqueza Territorial el catastro gráfico de Cortes:

 

De todos los archivos que contiene, solo necesitamos descompactar los 4 que se llaman CATAST_Pol_ParcelaRusti (prj, shp, shx y dbf). Incorporamos el .shp al proyecto y lo situamos el primero en el Panel de capas. Podemos eliminar del disco, grafico.zip.

Queda simbolizada como Símbolo único – Relleno sencillo, con color al azar. Lo dejamos así, porque necesitamos ver cuáles son parcelas del comunal.

Ahora hay que seleccionar las parcelas de nuestro interés y que están en el archivo Excel. Para ello, sobre CATAST_Pol_ParcelaRusti, botón derecho, Propiedades, vamos a la pestaña Uniones, adjuntando solo los datos que no vienen ya en el shapefile:

 

Al abrir la Tabla de atributos de CATAST_Pol_ParcelaRusti, podemos ver que se han incluido los tres atributos, pero dónde no hay unión, porque no son parcelas de comunal (no estaban en el Excel), ambos atributos se quedan vacíos (valor NULL):

 

Con “Seleccionar objetos por expresión”, nos quedamos con los que venían en el Excel 

 

De las 1120 parcelas, se seleccionan 147.

Sobre la capa, Exportar – Guardar objetos seleccionados como… creamos un nuevo shapefile "Cortes_ParcelaRustica_Adjudic" y lo añadimos al proyecto.

Esta nueva capa se situará por encima de la anterior y debe tapar exclusivamente aquellas parcelas que estaban seleccionadas en amarillo.

Podemos desconectar la visibilidad de la capa CATAST_Pol_ParcelaRusti y deseleccionar todos los objetos

Modificar la simbología de la nueva capa, dentro de símbolo único a Línea exterior sencilla, con un color que no se encuentre en las imágenes, azul, por ejemplo, grosor 0,4. También podemos etiquetarlas con el DNI del adjudicatario "ADJUDIC", pero con cuidado en la visibilidad con el control de escalas: Texto en tamaño 8, color azul oscuro, con buffer blanco, ubicación alrededor del centroide y representación en escalas mayores a 1:10.000:

 

Conectando de abajo a arriba la visibilidad de las imágenes NDVI, puede apreciarse la dinámica de los cultivos de nuestras parcelas a lo largo del año agrícola: verde en pleno vigor del cultivo hasta roja en suelo desnudo. Incluso las partes, dentro de la parcela, que no hay cultivo.

 

Pero "visualizar" no es suficiente. Podemos procesar las imágenes a nivel de parcela y asignarles valores cuantitativos (complejos de interpretar) o cualitativos (con una parametrización previa). Elegimos el 2º método, más interesante y de mejor capacidad de comunicación con los agricultores.

La cuestión técnica:

El rango 0 a 255 de las imágenes es un valor re-normalizado en su tratamiento, ya que las bandas originales recogen valores de reflectancia más amplios y que permiten un trabajo más rico utilizando técnicas propias de la teledetección. La re-normalización supone la reducción de información de cada pixel a un byte (28 = 256 valores posibles).

Una "inspección" por el método "de arriba abajo" de estos ráster permite conocer los valores de ambas bandas en distintas situaciones de cultivo.

La banda 1 de la imagen (B2 del satélite, rojo visible) es muy sensible al suelo desnudo, con valores mayores de 200 en terreno sin vegetación y menores de 100 con cultivo. La banda 2 (el índice) responde a la inversa.

Como se trata de determinar el ciclo de cambios en los cultivos del municipio, eminentemente regadíos de alta productividad, existen varios métodos para alcanzar ese propósito.

1. Desde las imágenes, trabajando con ráster:

Consiste en procesar las imágenes NDVI, mes a mes, para obtener una clasificación con 4 niveles*:

  • Suelo desnudo: clase 1
  • Cultivo en fase 1 (en floración): clase 2
  • Cultivo en fase 2 (en crecimiento): clase 3
  • Cultivo en pleno desarrollo: clase 4

    * Los rangos de las bandas de la imagen y su correspondencia en estados vegetativos deben ser establecidos con parcelas de "verdad-campo", incluso ajustarlo a los cultivos específicos del área y su ciclo agrícola. En este caso se trata de valores globales sin contrastar.

    La operación consiste en reclasificar las imágenes aprovechando estas relaciones:
  • Suelo desnudo:   banda 1 >=200 y banda 2 <=100
  • Cultivo en fase 1:            banda 1 >=200 y banda 2 >100 y <=200
  • Cultivo en fase 2:            banda 1 >=100 y banda 2 >200
  • Cultivo pleno:      banda 1 <100 y banda 2 >200

    Además, las imágenes originales NDVI tienen una extensión más allá de Navarra completa (son francamente grades), con lo que aprovecharemos la reclasificación para recortarlas al área de interés. Y pasaremos imágenes de 3 bandas a otras más sencilla, de una sola banda, con el resultado de la ecuación.

    2. Almacenando el promedio de valores de ambas bandas para cada parcela en el propio shapefile, creando posteriormente nuevos atributos como la asignación de clases anterior.

    Este método es bastante más farragoso que el anterior y se basa en el uso del mismo algoritmo, pero hay que ejecutarlo para cada imagen y banda, luego se trabaja el doble. Optamos por el primero y dejamos éste segundo como idóneo para pocas imágenes.

     
    1. Método 1, proceso:

1. Necesitamos un ráster que sirva de marco para recortar las imágenes. Como una vez concluido el proceso podremos eliminarlo, sirve cualquiera de las imágenes, pero ajustada al área de interés.

Desde la Caja de herramientas de procesos, con el algoritmo "Geometría mínima delimitadora", utilizando como capa de entrada CATAST_Pol_ParcelaRusti, sin campo de clases, generamos un tipo de geometría "Envolver (recuadro delimitador)" en archivo temporal (le asigna el nombre de Geometría delimitadora).

Desde el menú Ráster – Conversión – Rasterizar, desde Geometría delimitadora, con valor fijo 1, unidades del tamaño del ráster georreferenciadas: 10 x 10 m, utilizando la extensión de la propia capa, en archivo temporal (le asigna el nombre de OUTPUT).

Es conveniente, crear un directorio específico para las nuevas imágenes.

La nueva imagen, OUTPUT, podemos situarla en la parte final del Panel de capas, ya que no es más que una extracción de una ya existente.

2. En el menú: Ráster – Calculadora ráster, se muestran los archivos ráster del proyecto seguidos de @ y el número de banda. Construir la ecuación de reclasificación para la primera imagen: NDVI_20180916 y guardar el nuevo ráster como 20180916_NDVI.

La calculadora ráster evalúa condiciones igual que hace la vectorial, devolviendo un valor "1" cuando la condición es cierta y "0" cuando es falsa, resultado que habrá que multiplicar por el valor asignado a la clase correspondiente. Así, la clase 1 se construye con:

("NDVI_20180916@1">=200 and " NDVI_20180916@2"<=100) * 1

a la que seguirán las clases 2 a 4, que como son mutuamente excluyentes, su resultado cuando sea verdadero vale 1, se suma a las demás, que devolverán 0: clase 1 + clase 2 + clase 3 + clase 4.

NOTA: es recomendable auxiliarse de la calculadora ráster para construir las ecuaciones, pero resulta básico ayudarse de un documento dónde poder almacenar y modificar para no escribir otras 11 veces la ecuación. Con el bloc de notas documentar clases y ecuación en nuevo archivo que se guarda como Ecuacion_Calculadora_raster_NDVI.txt.

La ecuación completa es:

("NDVI_20180916@1" >= 200 and "NDVI_20180916@2" <= 100) * 1 +

("NDVI_20180916@2" >100 and "NDVI_20180916@2" <=200 and "NDVI_20180916@1">200) * 2 +

("NDVI_20180916@2" > 200 and "NDVI_20180916@1" >= 100) * 3 +

("NDVI_20180916@2" > 200 and "NDVI_20180916@1" < 100) * 4

Ver imagen:

 

3. Analizar el nuevo ráster 20180916_NDVI. Revisar el histograma para comprobar que el rango de valores está dentro de lo esperado. Estilo: unibanda pseudocolor, interpolación exacta (son valores discretos), paleta RdYlGn, modo de clasificación por intervalos iguales con 4 clases: 1 - rojo, 2 - amarillo, 3 – verde claro, 4 – verde oscuro. Completar la leyenda y situar el ráster en grupo Cultivos, en la parte superior del panel de capas (bajo las vectoriales) y comparar con NDVI_20180916. Se hacen apreciables estados intermedios de cultivo que en la imagen NDVI original pasaban inadvertidos.

 

4. Repetir el proceso con el resto de imágenes NDVI, con ayuda del archivo de texto: Edición – Reemplazar – 20180916 por 20181021, copiar y pegar en la calculadora ráster y crear el respectivo archivo 20181021_NDVI. Desplazar en el panel de capas a su grupo y copiar-pegar estilo desde 20180916_NDVI. Lo mismo para todas las demás imágenes con sus nombres correspondientes:

20181021, 20181230, 20190213, 20190305, 20190330, 20190429, 20190514, 20190628, 20190718 y 20190817

La colección de imágenes procesadas es impresionante. Si se quiere, ya podemos eliminar los originales NDVI del proyecto (y del disco, con lo que pesan).

En el siguiente post, continuamos el proceso.

SITNA

Datos adjuntos