Informe 1
Johan Urquijo
21/4/2020
Clasificación de cobertura del suelo en el área del municipio de Palmira
Clasificación mediante algoritmos no supervisados
Introducción
Datos y métodos
La zona de interés es el Valle del rio Cauca, especialmente toda la zona de emplazamiento de los ingenios azucareros. Al realizar una busqueda de imagenes Landsat 8 para esta zona, que cumplieran el criterio de tener un procentaje de nubosidad en la escena no mayor al 10%, se encontró que la mayoría de resultados arrojaban imagenes dañadas, tanto para el nivel de procesamiento L1 como L2. Teniendo esto en cuenta se optó por tomar una imagen Landsat 5 de la misma zona.
La escena en cuestión fue adquirida en Agosto de 1999 y su úlitmo procesamiento se realizó el 17 de Diciembre de 2017, tiene como coordenadas en el sistema mundial de referencia de Landsat el path “9” y el “row” 58. El porcentaje de nubosidad total de la escena es 4%, el RMSE geométrico en ambas direcciónes es de 4.597 metros.
A continuación se puede observar el resumen del propiedades para el conjunto de bandas agrupadas en en stack:
## class : RasterStack
## dimensions : 7001, 7771, 54404771, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 235485, 468615, 214785, 424815 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : LT05_L1TP//7_01_T1_B1, LT05_L1TP//7_01_T1_B2, LT05_L1TP//7_01_T1_B3, LT05_L1TP//7_01_T1_B4, LT05_L1TP//7_01_T1_B5, LT05_L1TP//7_01_T1_B6, LT05_L1TP//7_01_T1_B7
## min values : 0, 0, 0, 0, 0, 0, 0
## max values : 255, 255, 255, 255, 255, 255, 255Para mejorar la velocidad de procesamiento se recortó la escena al área que rodea el munició de Palmira, Valle del Cauca, las propiedades del conjunto de datos unz vez realizado el recorte fueron:
## class : RasterBrick
## dimensions : 1000, 1000, 1e+06, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 334995, 364995, 375015, 405015 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : C:/Users/jsurq/AppData/Local/Temp/RtmpUpeURq/raster/r_tmp_2020-04-22_095547_2528_98583.grd
## names : LT05_L1TP//7_01_T1_B1, LT05_L1TP//7_01_T1_B2, LT05_L1TP//7_01_T1_B3, LT05_L1TP//7_01_T1_B4, LT05_L1TP//7_01_T1_B5, LT05_L1TP//7_01_T1_B6, LT05_L1TP//7_01_T1_B7
## min values : 47, 18, 12, 8, 4, 128, 2
## max values : 255, 138, 161, 154, 255, 171, 255Debido a que esta imagen es de nivel de procesamiento L1, fue necesario calcular mediante la información de los metadatos la reflectancia en la parte alta de la atmosfera (TOA), en este caso como no se va a realizar un estudio multitemporal, el procesamiento hasta este nivel de reflectancia es suficiente.
A partir de esta imagen con reflectancia TOA, se calcularon las estádisticas unibanda. multibanda, los indices normalizados para agua, vegetación y construcciónes y se realizó una clasificación no supervisada mediante el procedimiento K-means. Es importante anotar que la banda térmica no se corrige radiométricamente para expresar la reflectancia sino para calcular los grádos celsius.
Para todos los procedimiento se utilizo el software estadístico R, (R Core Team, 2020).
Resultados y discución
Estadísticas y gráficos unibanda
En esta sección exploraremos la distribución, las medidas de tendencia central, de simetría y de dispersión de cada una de las bandas para sus niveles de reflectancia TOA, a través de este proceso verificaremos la calidad de los datos y en caso de ser necesario se hará modificación de los datos atípicos.
Lo primero que visualizaremos serán las bandas indivuaes en escala de gris:
IRC es la sigla para infrarrojo cercano e IROC es la sigla para infrarrojo de onda Corta
Se observa en la leyenda de todas las gráficas que todos los valores se encuentran entre 0 y 1 (a excepción de la banda térmica), lo que en terminos físicos es correcto teniendo en cuenta que la reflectancia es un porcentaje de la radiación total que recibe un objeto. Para verificar la distribución de los datos se calcularon los histrogramas y los estadísticos de tendencia central y dispersión para cada banda.
Esta es una tabla resumen de las estadísticas unibanda para todas las bandas:
## Mínimo Máximo Rango Media Mediana Moda Desv Estandar Asimetría Curtosis
## B1 0.07 0.41 0.34 0.10 0.10 0.09 0.01 3.22 27.75
## B2 0.05 0.46 0.41 0.09 0.08 0.08 0.01 2.75 18.49
## B3 0.03 0.46 0.43 0.07 0.07 0.06 0.02 2.05 7.45
## B4 0.02 0.53 0.51 0.26 0.26 0.27 0.06 -0.40 -0.17
## B5 0.00 0.59 0.59 0.15 0.14 0.12 0.04 0.92 2.15
## B6 22.43 40.37 17.94 27.33 26.83 25.53 2.44 1.00 0.34
## B7 0.00 0.82 0.83 0.08 0.07 0.05 0.04 1.32 3.54se observa que en general todas las bandas tienen valores correctos, que no son negativos o exceden los límites máximos adecuados (1 para visibles e infrarrojos y 255 para la térmica). Según la desviación estandar las bandas con mayor variabilidad son las del infrarrojo, dado que su desviación es un porcentaje mayor al 10% en relación a la media, especificamente el valor del porcentaje de la desviación frente a la media de cada banda es: 23.0769231% para B4, 28.5714286% para B5 y 57.1428571% para B7.
Las medidas de tendencia central son homogeneas para todas las bandas indicando que las distribuciones son monomodales. Ninguna banda tiene una distribución estadísticamente simétrica, todas estan desplazdas hacia la izquierda con excepcion de B4 que esta desplazada ligeramente hacia la derecha, lo que se evidencia tambien en que la moda sea mayor a la mediana y media. En terminos de curtosis nuevamente todas las dsitribuciones son leptocúrticas a excepción de B4 y B6 que estan cerca de tener una distribución mesocúrtica. Todas estas carácteristicas pueden verificarse al observar los histogramas graficados a continuación:
Desde este punto podemos empezar a determinar que con cierto grado de seguridad las bandas mas informativas serán las pertenecientes a los distintos tipos de infrarrojo, porque son la que tienen mayor variación, lo que unido a las bandas del espectro visible permitira discriminar mejor las distintas clases.
Estadísticas y gráficos multibanda
Ahora procederemos a calcular las estadísticas multibanda, se construiran mátrices de correlación y varianza-covarianza para las 7 bandas, y se graficaran algunas composiciones a color, que permitan diferenciar distintas clases de cobertura, para de esta manera usar esto mapas como linea de guía en el calculo de umbrales para la construcción de mascaras en el siguiente apartado.
La matriz de correlación de Pearson de las bandas es la siguiente:
## 7 x 7 Matrix of class "dtrMatrix"
## B1 B2 B3 B4 B5 B6 B7
## B1 1.00 . . . . . .
## B2 0.94 1.00 . . . . .
## B3 0.92 0.94 1.00 . . . .
## B4 -0.41 -0.32 -0.54 1.00 . . .
## B5 0.64 0.68 0.77 -0.36 1.00 . .
## B6 0.64 0.60 0.74 -0.67 0.68 1.00 .
## B7 0.77 0.74 0.88 -0.64 0.89 0.83 1.00En la matriz de correlación se observa que las bandas del espectro visibile estan altamente correlacionadas entre ellas, la banda 4 es la unica que tienen una tendencia negativa en la correlacion, además de tener los menores valores. Es importante tener en cuenta esta ausencia de correlación de B4 a la hora de incluirla en los siguientes pasos de procesamiento, la banda térmica tambien presenta un nivel medio de correlación con el resto de bandas.
La matriz de covarianza es esta:
## 7 x 7 Matrix of class "dtrMatrix"
## B1 B2 B3 B4 B5 B6 B7
## B1 0.000131 . . . . . .
## B2 0.000157 0.000212 . . . . .
## B3 0.000235 0.000303 0.000495 . . . .
## B4 -0.000284 -0.000280 -0.000727 0.003634 . . .
## B5 0.000281 0.000380 0.000660 -0.000847 0.001489 . .
## B6 0.017895 0.021237 0.040164 -0.099181 0.064409 5.977352 .
## B7 0.000370 0.000455 0.000822 -0.001619 0.001452 0.085903 0.001776Aquí podemos observar que la covarianza entre la banda térmica y las otras es mayor en comparación al resto de datos de la matríz, pero no se puede dejar de lado que este efecto se da por la escala de la variable mas no por un fenomeno de correlación de las bandas.
Composiciones a color
Para elegir cuales composiciones a color utilizar, se implementó el criterio del optimum index factor (OIF), para todas las posibles combinaciones de 3 bandas, y se seleccionaron las que según este estadístico son las mas informativas además de la composición a color natural.
A continuación podemos observar los valores de OIF para las 35 posibles combinaciones:
## Ba Bb Bc OIF
## 1 1 2 3 0.014
## 2 1 2 4 0.048
## 3 1 2 5 0.027
## 4 1 2 6 1.128
## 5 1 2 7 0.024
## 6 1 3 4 0.048
## 7 1 3 5 0.030
## 8 1 3 6 1.074
## 9 1 3 7 0.027
## 10 1 4 5 0.078
## 11 1 4 6 1.459
## 12 1 4 7 0.060
## 13 1 5 6 1.270
## 14 1 5 7 0.039
## 15 1 6 7 1.112
## 16 2 3 4 0.050
## 17 2 3 5 0.029
## 18 2 3 6 1.083
## 19 2 3 7 0.027
## 20 2 4 5 0.081
## 21 2 4 6 1.579
## 22 2 4 7 0.065
## 23 2 5 6 1.270
## 24 2 5 7 0.039
## 25 2 6 7 1.147
## 26 3 4 5 0.072
## 27 3 4 6 1.292
## 28 3 4 7 0.058
## 29 3 5 6 1.142
## 30 3 5 7 0.039
## 31 3 6 7 1.020
## 32 4 5 6 1.485
## 33 4 5 7 0.074
## 34 4 6 7 1.187
## 35 5 6 7 1.050A partir de esta tabla concluimos que según este criterio las combinaciones de las bandas 642, 654, 641 son las mas informativas. Esto sucede porque tanto B4 cómo B6 son las bandas menos correlacionadas con las demás, pero no es común usar está banda en composiciones de color, por lo que se graficaran estás tres combinaciones, pero tambien 4 combinaciones estandar.
A partir de estas composiciones podemos observar que cada una aporta una información diferente para la discriminación de las clases de cobertura. La combinación 432 es muy buena diferenciando las zonas costruidas en colores claros, de la vegetanción en rojos que van desde el palido hasta el intenso y los suelos descubiertos en color verde; aunque la combinación 453, se asemeja bastante esta última no discrimina bien la diferencia entre áreas construidas y suelo despejado, pero si permite distinguir muy bien la zona de relieve montañoso en el costado oriental de la imagen, algo que tambien logra muy bien la composición 451.
Hasta este momento podemos observar que se distinguen con cierta facilidad cuatro tipo de coberturas: Construido, Vegetación oscura, vegetación clara y suelos descubiertos.Aunque no se vea a simple vista hay dos afluentes de agua que se encuentran en la imágen, hasta el momento no ha sido facil distinguirlos, tal vez esto se pueda conseguir con las combinaciones, dadas por el OIF.
Estas composiciones no aportan mucha información nueva, en general permiten diferenciar muy bien el área construida de la vegetación, pero a la hora de distinguir el suelo descubierto de las áreas construidas presentan el mismo tono pero con diferente intensidad. La composición 645 fue adecuada para separar la zona del pie de monte en colores azules.
Calculo de indices de Cobertura
Se calcularon los indices normalizados para vegetación (NDVI), agua (MNDWI, AWEI) y construcciones (NDBI).
Primero visualizaremos NDVI
