1 INTRODUCCIÓN

El cálculo de valores estadísticos sobre datos raster es un aspecto esencial del manejo de este formato de datos. La obtención de estos valores estadísticos permite:

2 PASOS INICIALES

Para llevar a término este ejercicio, se requiere el uso de los siguientes paquetes.

library(raster)
## Loading required package: sp
#library(sp)
library(rasterVis)
## Loading required package: lattice
library(sf)
## Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE
#install.packages("kableExtra")
library(kableExtra)

También deberá establecerse el directorio de trabajo.

setwd("D:/G174_OCW/Lab_6_operaciones_estadisticas_imagenes_satelite/")

Desde ese directorio, se importará un fichero multibanda en formato GeoTIFF, solicitándo a R que muestre sus atributos.

imagen <- brick("./datos/analisis_estadistico.tif")
imagen
## class      : RasterBrick 
## dimensions : 1395, 1848, 2577960, 4  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 423405, 478845, 4452735, 4494585  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : analisis_estadistico.tif 
## names      : analisis_estadistico_1, analisis_estadistico_2, analisis_estadistico_3, analisis_estadistico_4 
## min values :                   9287,                   8398,                   7472,                   7127 
## max values :                  54511,                  43380,                  54694,                  65535

Para facilitar el trabajo con este fichero, se cambiará el nombre de las bandas

names(imagen) <- c("blue", "green", "red", "NIR")
names(imagen)
## [1] "blue"  "green" "red"   "NIR"

En las escenas Landsat originales, alrededor de la imagen aparece unas regiones con valores ausentes. Estos valores ausentes se codifican como 0. Si no se avisa al programa de esta circunstancia, considerará ese 0 como un valor numérico real. Por ello es necesario especificar esta circunstancia. Además, dado que se realizarán diferentes cálculos estadísticos, conviene evitar la conversión de los valores originales a notación científica.

NAvalue(imagen) <- 0                                          
options(scipen = 999)

Como es frecuente, conviene representar gráficamente el objeto importado mediante una imagen en falso color:

plotRGB(imagen, 
        r=4, g=3, b=2, 
        stretch = "lin", 
        axes = TRUE)

3 ANÁLISIS ESTADÍSTICO

La aplicación de diversos procedimientos estadísticos para el procesamiento de imágenes satelitales permite conocer, cuantificar y mejorar las características de la imagen, bien sea para ajustar la visualización o para orientar cualquier proceso posterior.

3.1 Cada banda por separado

El paquete raster proporciona la mayoría de las operaciones aritméticas para uso directo, como +, - , * y /, que habitualmente se realizan pixel a pixel. Por ejemplo, agregar dos rásteres requiere que cada píxel en el ráster A corresponda a un píxel en el ráster B; ambos serán sumados y almacenados en un nuevo ráster de las mismas dimensiones:

borrame1 <- imagen[[1]]+ imagen[[2]]
borrame1
## class      : RasterLayer 
## dimensions : 1395, 1848, 2577960  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 423405, 478845, 4452735, 4494585  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 17766, 97811  (min, max)

Además, se puede aplicar directamente un conjunto completo de funciones a los objetos ráster. Por ejemplo, para crear una nueva imagen que contenga la suma de todos los píxeles de un RasterStack podría usarse la función sum() de esta manera:

sumatotal <- sum(imagen)
sumatotal
## class      : RasterLayer 
## dimensions : 1395, 1848, 2577960  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 423405, 478845, 4452735, 4494585  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 32586, 216282  (min, max)

Otros argumentos son: log(), sqrt(), cumsum() o sign(). Algunas funciones, como cumsum(), devuelven más de una capa.

Los estadísticos que pueden extraerse son tanto de tendencia central como de dispersión. Los de tendencia central muestran alrededor de qué valores se sitúan la mayoría de los valores de los píxeles de los objetos ráster. El fundamental es la Media o promedio, valor que expresa el promedio de los niveles digitales de todos los pixeles de la imagen (Ecuación (2.1)).

\[\bar{x} = \frac{1}{n} \sum_{i=1}^{n}x_{i}\]

donde, \(x_i\) son los niveles digitales de los píxeles, y \(n\) el número total de pixeles en la imagen.

También es posible caracterizar las imágenes mediante estadísticos de dispersión. Uno de los básicos es la *Varianza**, un estadístico que muestra la dispersión de los niveles digitales de los pixeles alrededor de la media. Si los valores se concentran alrededor de la media la varianza es baja, por lo que su histograma será estrecho; en caso contrario, la varianza será. La fórmula de la varianza es la siguiente:

\[\sigma^{2} = \frac{\sum\limits_{i=1}^{n} \left(x_{i} - \bar{x}\right)^{2}} {n-1}\] donde, \(x_i\) es el nivel digital del pixel, \(n\) el número total de pixeles en la imagen y \(\bar{x}\) la media aritmética de la imagen.

Otro estadístico de dispersión ampliamente utilizado es la Desviación estándar, equivalente a la raíz cuadrada de la varianza

\[\sigma = \sqrt{\frac{\sum\limits_{i=1}^{n} \left(x_{i} - \bar{x}\right)^{2}} {n-1}}\]

El análisis puede comenzar con la función raster::summary (atención, esta función calcula los resultados a partir de una muestra de píxeles, no de todos el conjunto).

summary(imagen)
## Warning in .local(object, ...): summary is an estimate based on a sample of 100000 cells (3.88% of all cells)
##          blue    green      red   NIR
## Min.     9389  8606.00  7741.00  7483
## 1st Qu. 11094 10995.00 11491.00 15842
## Median  11703 11839.00 12831.00 17219
## 3rd Qu. 12342 12784.75 14260.75 18949
## Max.    28056 29686.00 32035.00 35926
## NA's    31760 31760.00 31760.00 31760

Estos estadísticos pueden también obtenerse de manera individual

minValue(imagen)                                
## analisis_estadistico_1 analisis_estadistico_2 analisis_estadistico_3 
##                   9287                   8398                   7472 
## analisis_estadistico_4 
##                   7127
maxValue(imagen)                                
## analisis_estadistico_1 analisis_estadistico_2 analisis_estadistico_3 
##                  54511                  43380                  54694 
## analisis_estadistico_4 
##                  65535

El cálculo de estos valores estadísticos se puede también obtener mediante la función raster::cellStats del paquete raster.

cellStats(imagen, "mean")
##     blue    green      red      NIR 
## 11817.60 11968.05 12940.29 17474.32
cellStats(imagen, "var")
##    blue   green     red     NIR 
## 1275728 2076143 3968944 5105280
cellStats(imagen, "sd")
##     blue    green      red      NIR 
## 1129.481 1440.883 1992.221 2259.487
cellStats(imagen, "mad")                                                                  # Desviación mediana absoluta (Median-absolute deviation, MAD)
##      blue     green       red       NIR 
##  925.1424 1325.4444 2054.8836 2247.6216
cellStats(imagen, function(x,...) quantile(x, probs=c(0.05, 0.25, 0.5, 0.75, 0.95),...))  # Cuantiles: 5%, 25%, 50%, 75% and 95%
##      blue green   red   NIR
## 5%  10309  9911  9957 14380
## 25% 11093 10992 11485 15838
## 50% 11702 11841 12832 17212
## 75% 12343 12785 14260 18942
## 95% 13710 14323 16186 21390

La función cellStats usa una muestra para calcular los estadísticos fundamentales (sum, mean, min, max, sd, skew y rms).

A su vez, estos valores estadísticos pueden transformarse en una tabla con el consiguiente código.

imagen_promedio <- cellStats(imagen, "mean")
imagen_varianza <- cellStats(imagen, "var")
imagen_desviacion_tipica <- cellStats(imagen, "sd")

estadisticos <- rbind(imagen_promedio, imagen_varianza, imagen_desviacion_tipica)

estadisticos <- cbind(estadisticos, 
                      apply(estadisticos,1, FUN = min), apply(estadisticos, 1, FUN = max), apply(estadisticos,1, FUN = mean))

rownames(estadisticos) <- c("Media", "Varianza" , "Desviación estandar")
colnames(estadisticos) <- c(colnames(estadisticos)[1:4], "Mínimo", "Máximo", "Medio")

kable(round(estadisticos,3), 
      caption="Resumen estadísticas de las bandas")
Resumen estadísticas de las bandas
blue green red NIR Mínimo Máximo Medio
Media 11817.599 11968.047 12940.292 17474.324 11817.599 17474.324 13550.066
Varianza 1275728.032 2076142.997 3968943.740 5105280.315 1275728.032 5105280.315 3106523.771
Desviación estandar 1129.481 1440.883 1992.221 2259.487 1129.481 2259.487 1705.518

Las desviaciones estándar son más elevadas en las bandas NIR, SWIR1 y SWIR2. El hecho de que la desviación estándar sea mayor indica que los datos están muy dispersos alrededor de la media, lo anterior se puede relacionar con imágenes bien contrastadas.

El análisis numérico debe ser complementado con un análisis gráfico, usando para ello las diferentes herramientas disponibles en R. Por ejemplo, los histogramas muestran la distribución de frecuencias de los valores de los píxeles de una imagen raster. Con su ayuda se puede cuantificar la frecuencia de niveles digitales, ya que el histograma indica el número de pixeles presentes en cada rango. Además, con este gráfico se puede extraer una primera aproximación de las posibles clases temáticas y definir la calidad de la imagen por medio del contraste: una imagen que genere un histograma muy estrecho corresponde a una imagen con un pobre contraste visual.

El rango de los valores que puede tener un pixel se muestra en el eje X, mientras que su frecuencia aparece en el eje Y. En las siguientes líneas de código se representará el histograma correspondiente a una única banda.

hist(imagen[[2]], 
     maxpixels=1000000, 
     ylim = c(0, 0.00030),
     plot=TRUE, 
     prob = TRUE, 
     ylab="Frecuencia", 
     col="darkblue", 
     xlab="Banda blue")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 39%
## of the raster cells were used. 1000000 values used.
# Generación de una curva de densidad superpuesta al histograma
ds <- density(imagen[[2]], 
              plot = FALSE)
lines(ds, 
      col = "red", 
      lwd = 2)

Normalmente, es más usual representar simultáneamente todas las bandas disponibles en el objeto ráster.

par(mfrow=c(2,2))
hist(imagen[[1]], maxpixels=1000000, plot=TRUE, ylab="Frecuencia", col="darkblue", xlab=" Banda blue")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 39%
## of the raster cells were used. 1000000 values used.
hist(imagen[[2]], maxpixels=1000000, plot=TRUE, ylab="Frecuencia", col="green",xlab=" Banda green")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 39%
## of the raster cells were used. 1000000 values used.
hist(imagen[[3]], maxpixels=1000000, plot=TRUE, ylab="Frecuencia", col="red", xlab="Banda red")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 39%
## of the raster cells were used. 1000000 values used.
hist(imagen[[4]], maxpixels=1000000, plot=TRUE, ylab="Frecuencia", col="darkred", xlab="Banda NIR")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 39%
## of the raster cells were used. 1000000 values used.

par(mfrow=c(1,1))

Una interesante alternativa para la visualización de objetos ráster es el paquete rasterVis. La función rasterVis::histogram() dibuja tanto el histograma de una única banda, como los histogramas de todas al bandas simultáneamente.

  • Histograma de una única banda
histogram(imagen[[4]], 
          main = "Histograma de la banda NIR", 
          xlab = "Valores de los píxeles", 
          ylab = "Frecuencia", col ="wheat")

  • Histograma de todas las bandas simultáneamente
histogram(imagen, 
          main = "Histograma de todas las bandas", 
          xlab = "Valores de los píxeles", 
          ylab = "Frecuencia", 
          col = "wheat")

Un gráfico de densidad muestra la distribución de los valores de una imagen raster con intervalos continuos (a diferencia del histograma), y es muy útil para determinar la forma de la distribución de esa imagen. Rl pico del gráfico de densidad muestra el intervalo de los datos con una mayor frecuencia (la moda). Para ello se usa la función rastervis::densityplotdel paquete rasterVis

densityplot(imagen,
            xlab = "Valores",
            main = "Gráfico de densidad de cada banda")

El gráfico de violín es una combinación del gráfico de caja y bigotes y del gráfico de densidad rotado. El primero muestra la distribución de los datos a partir de los cuartiles. La caja incluye todos los datos entre los cuartiles superior e inferior; los bigotes representan la variabilidad de los datos más alllá de los cuartiles superior e inferior. Además, el gráfico de caja y bigotes incluye otros estadísticos como la media, el rango intercuartílico, etc..

El gráfico de densidad muestra la distribución de los datos en términos de densidad de píxeles con respecto a sus valores. Por lo tanto, un gráfico de violín muestra a la vez un resumen de los datos y la distribución completa de las variables. Para realizarlo sobre datos raster acudiremos a la funcion bwplotdel paquete rasterVis.

bwplot(imagen, 
       main = "Gráfico de violín")

3.2 Todas las celdas o píxeles

También pueden calcularse esos mismos estadísticos a nivel de píxel de todas las capas de una imagen multibanda, creando una nueva capa raster que almacena esos nuevos datos. Esto se realiza mediante la función raster::calcdel paquete raster. El estadístico deseado se incorpora con el argumento fun. Por ejemplo:

promedio_imagen <- calc(imagen, 
                        fun = mean, 
                        na.rm = TRUE,
                        filename = "promedio_imagen.tif", 
                        overwrite = TRUE)
promedio_imagen
## class      : RasterLayer 
## dimensions : 1395, 1848, 2577960  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 423405, 478845, 4452735, 4494585  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : promedio_imagen.tif 
## names      : promedio_imagen 
## values     : 8146.5, 54070.5  (min, max)
plot(promedio_imagen, main = "promedio de todas las bandas, axes = TRUE")

La distribución de los valores medios puede ser verificada recurriendo a un histograma

histogram(promedio_imagen, 
          col = "wheat", 
          main = "Histograma de valores medios de todas las bandas")

Es recomendable usar la función raster::calc para cálculos complejos, porque una llamada a esta función no implica la generación de archivos temporales. Una función personalizada podría escribirse así

borrame2 <- calc(imagen, 
                 fun = function(x) {
                   2*x + 0.5*x^2 - 10})

Como en casos anteriores, se puede crear una función y aplicarla a todos los píxeles de un RasterStack o RasterBrick. La sintaxis siempre debe enunciarse de la siguiente manera: funcion(variable) {formula matemática}. A continuación se define una función multiplicando todos los píxeles por 100, que es aplicada posteriormente al objeto ráster.

func <- function(x) {x * 100}

imagen_multiplicada <- calc(imagen, func)
imagen_multiplicada
## class      : RasterBrick 
## dimensions : 1395, 1848, 2577960, 4  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 423405, 478845, 4452735, 4494585  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      :    blue,   green,     red,     NIR 
## min values :  928700,  839800,  747200,  712700 
## max values : 5451100, 4338000, 5469400, 6553500

3.3 Medidas de relación entre bandas

La aplicación de otros estadísticos permite conocer la relación entre las bandas de una imagen. Los estadísticos multibanda más habituales son la covarianza y la correlación, normalmente acompañados de un diagrama de dispersión por pares.

Covarianza mide la variación conjunta de los niveles digitales de los pixeles de dos bandas alrededor de una media común. Es esencialmente una medida de la varianza entre dos variables. Toma valores positivos cuando ambas variables tienden a aumentar o disminuir a la vez y valores negativos cuando una variable tiene a incrementar mientras la otra disminuye o viceversa. Sin embargo, no evalúa la dependencia entre variables y se mide en las mismas unidades que las que aparecen medidas dichas variables: la varianza puede tomar cualquier valor positivo o negativo.

  • Covarianza positiva: Indica que dos variables tienden a moverse en la misma dirección.
  • Covarianza negativa: Revela que dos variables tienden a moverse en direcciones inversas.

La fórmula de cálculo es la siguiente

\[cov_{x,y} = \frac{\sum\limits_{i=1}^{n}{(x_i-\overline{x}) \cdot (y_i-\overline{y})} }{n-1}\]

donde \(X_i\) es el nivel digital del pixel en la banda \(X\), \(Y_i\): Nivel digital del pixel en la banda \(Y\), \(\bar{x}\) la media aritmética de la banda \(X\), \(\bar{y}\)y la media aritmética de la banda \(Y\) y \(m\) el número total de píxeles en la imagen

Para el cálculo de la matriz de covarianza entre las bandas de una misma imagen es necesario transformar cada una de ellas en una matriz.

cov <- cov(as.matrix(imagen), use = "complete.obs")
round(cov*1000,2)
##             blue      green        red        NIR
## blue  1275728032 1555661556 1941097249 1398264024
## green 1555661556 2076142997 2753069621 2236621479
## red   1941097249 2753069621 3968943740 3294359200
## NIR   1398264024 2236621479 3294359200 5105280315

Coeficiente de correlación también mide la relación entre los niveles digitales de dos bandas, pero toma valores entre entre -1 y 1, siendo rxy=1 una relación perfecta entre las dos bandas, rxy=-1 una relación inversa entre las dos bandas y rxy=0 no existe relación entre las dos bandas.

\[\rho = \frac{\text{cov}(X,Y)}{\sigma_x \sigma_y}\]

siendo \({cov}(X,Y)\) la covarianza de las bandas X y Y, \(sigma_x\) la desviación estándar de la banda \(X\) y \(sigma_y\) la desviación estándar de la banda \(Y\).

El coeficiente de correlación puede calcularse entre dos o más bandas, aunque el fichero ráster debe transformarse en un dataframe.

x <- as.data.frame(na.omit(imagen))
cor(x$blue, x$NIR, method="spearman", use='complete.obs')
## [1] 0.5585544
cor(x, method="spearman", use='complete.obs')
##            blue     green       red       NIR
## blue  1.0000000 0.9551156 0.8714753 0.5585544
## green 0.9551156 1.0000000 0.9639095 0.6937604
## red   0.8714753 0.9639095 1.0000000 0.7254449
## NIR   0.5585544 0.6937604 0.7254449 1.0000000

Un diagrama de dispersión por pares constituye una buena alternativa a los valores numéricos. La comparacion entre los valores espectrales de dos bandas de una misma imagen raster puede realizarse mediante la función pairsdel paquete raster. Esta función devuelve un gráfico de dispersión entre las dos bandas especificadas (que muestra el valor de un determinado pixel en ambas bandas y el grado de relacion entre ambas), sus histogramas y el valor del coeficiente de correlación entre ambas.

pairs(imagen[[1:4]], 
      hist=TRUE,                                # Muestra el histograma
      cor=TRUE,                                 # Muestra el coeficiente de correlación
      main = "Bandas Red versus NIR",
      use="pairwise.complete.obs", 
      maxpixels=100000)
## Warning in graphics::par(usr): argument 1 does not name a graphical parameter

## Warning in graphics::par(usr): argument 1 does not name a graphical parameter

## Warning in graphics::par(usr): argument 1 does not name a graphical parameter

## Warning in graphics::par(usr): argument 1 does not name a graphical parameter

## Warning in graphics::par(usr): argument 1 does not name a graphical parameter

## Warning in graphics::par(usr): argument 1 does not name a graphical parameter

## Warning in graphics::par(usr): argument 1 does not name a graphical parameter

## Warning in graphics::par(usr): argument 1 does not name a graphical parameter

## Warning in graphics::par(usr): argument 1 does not name a graphical parameter

## Warning in graphics::par(usr): argument 1 does not name a graphical parameter

La funcion rasterVis::splomdel paquete rasterVis también ofrece el gráfico de dispersión y la curva de densidad de todas las bandas al mismo tiempo.

splom(imagen, 
      main = "Gráfico de dispersion con todas las bandas")

Existen otras opciones para la representación gráfica de la relación entre capas de una imagen

library(corrplot)
## corrplot 0.92 loaded
matriz.corr <- round(cor(x, method="spearman", use='complete.obs'), 2)
corrplot(matriz.corr, 
         type = "upper", 
         order = "hclust", 
         tl.col = "black", 
         tl.srt = 45)

tmp1 <- as.matrix(sample(imagen,500000))
library(PerformanceAnalytics)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
chart.Correlation(tmp1, 
                  histogram=TRUE, 
                  pch=19)
## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

4 CONSULTAS ESPACIALES: EXTRACCIÓN DE VALORES DE LAS IMÁGENES

En el caso de las imágenes de satélite, se pueden obtenr los valores correspondientes a una determinada función, por ejemplo, la media, superponiendo un fichero vectorial con polígonos a la imagen en formato raster, también mediante la función raster::extractdel paquete raster. Esta función permite que un objeto, bien un matrix o un data.frame, con dos columnas (correspondientes a las coordenadas x e y), o bien cualquier otro objeto espacial creado con el paquete sp (por ejemplo, un SpatialPoints*, SpatialPolygons*, SpatialLines o Extent) como entrada.

Las funciones disponibles deben incluirse con el argumento fun. Como en casos anteriores, tanto la imagen raster como el fichero vectorial deberían tener el mismo CRS. A continuación se van a plantear varios ejemplos.

4.1 Extracción de los valores correspondientes a puntos.

En caso de no disponer de esos puntos, se pueden crear esos puntos utilizando diferentes procedimientos.

imagen

  • Un único punto
xy <- data.frame(x = 453405, y = 4472735)                     # Creación de un dataframe con un único punto 
xy <- st_as_sf(xy, coords = c("x","y"))

Dado que xy es un dataframe, y no un objeto espacial, es necesaria su transformación, adquiriendo además la misma proyección del objeto ráster con el que se está trabajando, y en caso necesario, se convertirán en un fichero vectorial con formato shp.

crs(imagen)
## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["unknown",
##     BASEGEOGCRS["unknown",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ID["EPSG",6326]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8901]]],
##     CONVERSION["UTM zone 30N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]],
##         ID["EPSG",16030]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]
xy <- st_as_sf(xy, coords = c("x","y"))
st_crs(xy) <- st_crs(imagen)                                  # Transformación del dataframe en un objeto espacial con sistema de proyección)
xy
## Simple feature collection with 1 feature and 0 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 453405 ymin: 4472735 xmax: 453405 ymax: 4472735
## Projected CRS: +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
##                 geometry
## 1 POINT (453405 4472735)
st_write(xy, "D:/G174_OCW/Lab_6_operaciones_estadisticas_imagenes_satelite/datos/puntos.shp", append = FALSE)
## Deleting layer `puntos' using driver `ESRI Shapefile'
## Writing layer `puntos' to data source 
##   `D:/G174_OCW/Lab_6_operaciones_estadisticas_imagenes_satelite/datos/puntos.shp' using driver `ESRI Shapefile'
## Writing 1 features with 0 fields and geometry type Point.

Otra posibilidad para que este objeto adquiera el sistema de proyección del objeto ráster es el siguiente.

#compareCRS(xy, imagen)
#xy <- st_set_crs(xy, "EPSG:32618")

Tras ello, se deben superponer los datos vectoriales a la imagen para visualizar la ubicación del punto.

plotRGB(imagen, r=4, g=3, b=2, stretch = "lin", axes = TRUE)
plot(xy, 
     col = "yellow",                # Color de los puntos
     cex = 3,                       # Tamaño del símbolo
     pch = 19,                      # Tipo de puntos
     add = TRUE)

extract(imagen,                                               # Objeto ráster
        xy)                                                   # Puntos de los que se va a extraer los valores de un pixel
##       blue green   red   NIR
## [1,] 12923 13080 13872 16965
  • 20 puntos distribuidos aleatoriamente por la zona de análisis. Estos puntos se situarán dentro de las coordenadas del objeto ráster (argumentos xmin, xmax, ymin, ymin`).
puntos <- data.frame(x = runif(20, xmin(imagen), xmax(imagen)), 
                 y = runif(20, ymin(imagen), ymax(imagen)))
puntos <- st_as_sf(puntos, coords = c("x","y"))
st_crs(puntos) <- st_crs(imagen)                                  

plotRGB(imagen, r=4, g=3, b=2, stretch = "lin", axes = TRUE)
plot(puntos, 
     col = "yellow",                # Color de los puntos
     cex = 3,                       # Tamaño del símbolo
     pch = 19,                      # Tipo de puntos
     add = TRUE)

A partir de este momento ya podemos extraer los valores de todas las bandas de la imagen en las localizaciones especificadas mediante la función raster::extractdel paquete raster, creando con ello un dataframe

puntos_valores <- extract(imagen, 
                          puntos)          
puntos_valores
##        blue green   red   NIR
##  [1,] 11492 11421 12059 16215
##  [2,] 11755 11578 13156 15606
##  [3,] 12189 12643 13693 16979
##  [4,] 12128 12342 14086 17321
##  [5,] 11479 11826 13172 17468
##  [6,] 13182 14567 16652 20829
##  [7,] 10345  9890 10084 15336
##  [8,] 11728 12518 14700 19368
##  [9,] 10181  9609  9593 14632
## [10,] 10963 10801 11598 16139
## [11,] 10190  9723  9832 15096
## [12,] 11317 11427 12470 16288
## [13,] 13965 14226 15449 17770
## [14,] 10850 10828 11152 16643
## [15,] 13033 13072 13948 17068
## [16,] 10820 10836 11458 16430
## [17,] 11770 12382 14187 18929
## [18,] 11071 10756 10788 15568
## [19,] 12366 12680 13875 17787
## [20,] 13094 12867 13432 16654

Estos valores extraidos de la imagen raster pueden ser añadidos al fichero vectorial original (puntos.shp) para crear un objeto espacial y posteriormente pueden ser salvados otra vez como fichero vectorial o como fichero de datos (*.csv).

puntos <- cbind(puntos, 
                puntos_valores)
puntos
## Simple feature collection with 20 features and 4 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 423648.3 ymin: 4456799 xmax: 478219.8 ymax: 4494466
## Projected CRS: +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## First 10 features:
##     blue green   red   NIR                 geometry
## 1  11492 11421 12059 16215 POINT (449432.8 4472243)
## 2  11755 11578 13156 15606 POINT (445743.2 4475348)
## 3  12189 12643 13693 16979 POINT (438897.2 4460932)
## 4  12128 12342 14086 17321 POINT (426648.3 4457008)
## 5  11479 11826 13172 17468 POINT (477553.6 4486426)
## 6  13182 14567 16652 20829 POINT (447533.9 4481984)
## 7  10345  9890 10084 15336 POINT (443314.9 4494466)
## 8  11728 12518 14700 19368 POINT (478219.8 4476422)
## 9  10181  9609  9593 14632 POINT (435388.2 4483862)
## 10 10963 10801 11598 16139 POINT (425756.7 4485633)

En caso de querer exportarlos como fichero de datos, el código sería el siguiente:

write.table(puntos, file = "./datos/fichero_puntos.csv", 
            append = FALSE, 
            sep = ",", 
            row.names = FALSE, 
            col.names = TRUE)

Los datos extraídos pueden ser representados gráficamente para visualizar la ditribución de los datos en un píxel particular de todas las bandas.

plot(puntos_valores[1, ],            # Indicación para la primera columa
     type = "o", 
     axes = TRUE,
     main = "Representacion espectral de los puntos extraidos")

Se puede mejorar esta representacion gráfica incluyendo más detalles.

plot(puntos_valores[1, ],            # Indicación para la primera columna
     type = "o",                    # Tipo de símbolo
     axes = FALSE,
     ylim = c(7000,16000),              # Valores límite del eje Y
     ylab = "Valores de reflectancia", xlab = "Nombre de las bandas", col = "blue",
     main = "Representacion espectral de los puntos extraidos")

# Adición de los ejes al dibujo existente, pero con nuevos valores
axis(side = 1, at = 1:4, labels = colnames(puntos_valores), tick =  TRUE)
axis(side = 2, at = NULL, labels = TRUE, tick =  TRUE)

# Adición de otras filas de datos al gráfico existente
lines(puntos_valores[2, ], type = "o", col = "red")
lines(puntos_valores[3, ], type = "o", col = "purple")
lines(puntos_valores[4, ], type = "o", col = "green")

# Leyenda
legend("topleft", 
       legend = c(1,2,3,4), 
       lty = 1, 
       bty = "n", 
       col = c("blue", "red", "purple", "green"), 
       y.intersp = 0.5, 
       cex = 0.7,
       title = "ID de los puntos espaciales")

4.2 Extracción de valores correspondientes a polígonos: estadísticas zonales.

Por estadísticas zonales se entiende el cálculo de diversos estadísticos para una o varias zonas previamente definidas, a partir de una imagen ráster existente. Por ejemplo, se puede conocer la pendiente media de varias cuencas hidrográficas combinando un modelo digital de elevaciones (DEM) en formato ráster y un fichero vectorial con los límites de esas cuencas. El resultado es una tabla que contiene valores estadísticos para cada una de las cuencas hidrográficas.

Lo más habitual es que estemos interesados en extraer los valores correspondientes a una o varias zonas de interés (ROI), que habitualmente conforman objetos espaciales en formato vectorial (polígonos). Estos polígonos pueden ser digitalizados u obtenidos de una fuente de información externa.

En caso de tener que digitalizar los polígonos sobre la imagen deberemos utilizar este código

#plotRGB(imagen[[4:2]], stretch = "lin", axes = TRUE)

# Digitalización de cada polígono. Esto deberá repetirse tantas veces como polígonos deseemos dibujar. Una vez se llama a la función, comienza la digitalización de los vértices de los polígonos. Para finalizar, se hace click con el botón de la derecha del ratón y se selecciona stop

#poligono1 <- drawPoly(sp = TRUE, col='yellow', lwd=2)
#st_crs(poligono1) <- st_crs(imagen)

# Creación de un único fichero con todos los polígonos y grabación del mismo
#zonas <- rbind(poligono1,poligono2, poligono3, poligono4, poligono5, poligono6)
#st_write(zonas, "./datos/zonas.shp")

# Verificación de la posición de los polígonos
# plotRGB(imagen[[4:2]], stretch = "lin", axes = TRUE)
# plot(zonas, col = "yellow", add =TRUE)

En las siguientes lineas se importará un fichero vectorial que incluye varios poligonos para calcular posteriormente la moda de todos los píxeles situados dentro de cada polígono, en cada una de las bandas.

zonas <- st_read("./datos/zonas.shp")
## Reading layer `zonas' from data source 
##   `D:\G174_OCW\Lab_6_operaciones_estadisticas_imagenes_satelite\datos\zonas.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 4 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 427981.7 ymin: 4456098 xmax: 476274.3 ymax: 4492225
## Projected CRS: Undefined Cartesian SRS with unknown unit
class(zonas)
## [1] "sf"         "data.frame"
compareCRS(zonas, imagen)
## [1] FALSE

Dado que el nuevo fichero vectorial no comparte el mismo CRS que la imagen, es requisito ineludible asignarle el mismo CRS.

st_crs(zonas) <- 32630
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that
st_crs(zonas)
## Coordinate Reference System:
##   User input: EPSG:32630 
##   wkt:
## PROJCRS["WGS 84 / UTM zone 30N",
##     BASEGEOGCRS["WGS 84",
##         ENSEMBLE["World Geodetic System 1984 ensemble",
##             MEMBER["World Geodetic System 1984 (Transit)"],
##             MEMBER["World Geodetic System 1984 (G730)"],
##             MEMBER["World Geodetic System 1984 (G873)"],
##             MEMBER["World Geodetic System 1984 (G1150)"],
##             MEMBER["World Geodetic System 1984 (G1674)"],
##             MEMBER["World Geodetic System 1984 (G1762)"],
##             MEMBER["World Geodetic System 1984 (G2139)"],
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ENSEMBLEACCURACY[2.0]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 30N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Between 6°W and 0°W, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Burkina Faso. Côte' Ivoire (Ivory Coast). Faroe Islands - offshore. France. Ghana. Gibraltar. Ireland - offshore Irish Sea. Mali. Mauritania. Morocco. Spain. United Kingdom (UK)."],
##         BBOX[0,-6,84,0]],
##     ID["EPSG",32630]]

A continuación, se superpone la imagen raster y el fichero vectorial con los polígonos.

plotRGB(imagen[[4:2]], stretch = "lin", axes = TRUE)
plot(zonas, 
     col = "yellow", 
     add =TRUE)

Para obtener todos los valores de los píxeles ubicados dentro de los polígonos, se usa la función raster::extract(). Esta función requiere:

  • El ráster del que deseamos extraer valores
  • La capa vectorial que contiene los polígonos que deseamos utilizar como límite o límites

NOTA: Podemos decirle que almacene los valores de salida en un data.frame usando df=TRUE (opcional, el valor predeterminado es NO devolver un marco de datos).

borrame1 <- raster::extract(imagen, 
                           zonas, 
                           df = TRUE)

Cuando usamos la función extract(), R extrae el valor de cada píxel ubicado dentro de los límites de los polígonos, en este caso el objeto zonas; la función extrajo valores de 540314 píxeles. Para ver los primeros datos del dataframe

head(borrame1)
##   ID  blue green   red   NIR
## 1  1 10288  9924 10125 16011
## 2  1 10538 10314 10848 15986
## 3  1 10347 10006 10387 15088
## 4  1  9982  9539  9548 15250
## 5  1 10074  9706  9843 15219
## 6  1 10335 10002 10438 15113

Podemos crear un histograma de los valores de altura de los árboles dentro del límite para comprender mejor la estructura o la distribución de la altura de los árboles.

hist(borrame1$blue, 
     main="Histograma",
     col="springgreen",
     xlab="Etiqueta eje X", 
     ylab="Etiqueta eje y")
abline(v = mean(borrame1$blue), 
       lwd = 2) ## Mean line

También podemos usar la función summary() para ver estadísticas descriptivas que incluyen valores mínimo, 1er cuartil, mediana, promedio, 3er cuartil y máximo de la banda analizada.

summary(borrame1$blue)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    9404   10949   11555   11671   12172   28058

También pueden calcularse los 5 números de Tukey (valor minimo, percentil 25, mediana (percentil 50), percentil 75 y valor máximo)

fivenum(borrame1$red)
## [1]  7728 11253 12866 14357 34429

4.3 Resumir valores ráster extraídos

A menudo queremos extraer valores de resumen de un ráster. Podemos decirle a R el tipo de resumen estadístico que nos interesa usando el método fun=. Extraigamos un valor de altura media para nuestro AOI.

Extracción de los valores promedios de los píxeles localizados dentro de los polígonos

borrame2 <- raster::extract(imagen, 
                            zonas,
                            fun = mean, 
                            df=TRUE)
borrame2
##   ID     blue    green      red      NIR
## 1  1 10646.36 10350.68 10758.91 15819.12
## 2  2 11858.66 11697.60 12214.97 16219.51
## 3  3 11857.08 12422.25 14107.29 18270.96
## 4  4 11925.88 12245.39 13494.40 18048.09

Si se qusiera calcular los 5 números de Tukey para todas las bandas habría que recurrir a la siguiente sintaxis.

sapply(borrame2[c('blue', 'green', 'red', 'NIR')], fivenum)
##          blue    green      red      NIR
## [1,] 10646.36 10350.68 10758.91 15819.12
## [2,] 11251.72 11024.14 11486.94 16019.32
## [3,] 11857.87 11971.49 12854.68 17133.80
## [4,] 11892.27 12333.82 13800.84 18159.52
## [5,] 11925.88 12422.25 14107.29 18270.96