INTRODUCCION

En este tutorial se trabajan los principales fundamentos, paquetes y atributos de metadatos/raster necesarios para trabajar en R.

Los paquaquetes de R que se utilizan son dplyr, ggplot2, raster y rgdal

La imagen que se va a utilizar en este tutorial pertenece a un recorte ubicado en el municipio de Villavicencio, departamento del Meta - Colombia

library(raster)
## Warning: package 'raster' was built under R version 3.4.4
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
##  Path to GDAL shared files: C:/Users/Natalia/Documents/R/win-library/3.4/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/Natalia/Documents/R/win-library/3.4/rgdal/proj
##  Linking to sp version: 1.2-7

Visualizacion de Atributos del archivo Raster

Esta leccion trabaja con una serie de GeoTIFF. GeoTIFF es un formato que contiene etiquetas integradas con metadatos a cerca de un dato raster. Para obtener la informacion sobre el dato raster antes de ser leido en R se utiliza la funcion GDALinfo. Esta funcion es buena antes de importar los datos que se van a trabajar.

# Seleccion de directorio de trabajo
setwd("C:\\Maestria en Geomatica\\Percepcion remota avanzada\\Talleres")
getwd()
## [1] "C:/Maestria en Geomatica/Percepcion remota avanzada/Talleres"
# Visualizar atributos de la imagen
GDALinfo("Landsat_recorte_B2.tif")
## rows        2618 
## columns     3571 
## bands       1 
## lower left origin.x        616785 
## lower left origin.y        397515 
## res.x       30 
## res.y       30 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## file        Landsat_recorte_B2.tif 
## apparent band summary:
##    GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float32           TRUE        -999          1       3571
## apparent band statistics:
##   Bmin  Bmax    Bmean      Bsd
## 1    0 35668 8550.375 2247.965
## Metadata:
## AREA_OR_POINT=Area

Para guardar la informacion obtenida en R

info_B2 <- capture.output(GDALinfo("Landsat_recorte_B2.TIF"))

Abrir un raster en R

Para importar el raster en R y explorar los metadatos mas en detalle se utiliza la funcion raster()

  1. Descargar el archivo raster en R y ver la estructura de los datos.
info_B2 <- raster("Landsat_recorte_B2.TIF")
info_B2
## class       : RasterLayer 
## dimensions  : 2618, 3571, 9348878  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 616785, 723915, 397515, 476055  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\Maestria en Geomatica\Percepcion remota avanzada\Talleres\Landsat_recorte_B2.TIF 
## names       : Landsat_recorte_B2 
## values      : 0, 35668  (min, max)

La informacion obtenida incluye los valores de minimos maximos, numero de filas y columnas, resolucion, sistema de coordenadas; pero no otros datos estadisticos.

Al igual que otras estructuras de datos en R como vectores y columnas data frames, la estadistica descriptiva de datos raster puede ser recuperada como:

summary(info_B2)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (1.07% of all cells)
##         Landsat_recorte_B2
## Min.                     0
## 1st Qu.               8647
## Median                8920
## 3rd Qu.               9198
## Max.                 35248
## NA's                     0

Sin embargo, la nota de advertencia se debe a que R calcula estas estadisticas utilizando todas las celdas en el raster, esto puede tomar una muestra aleatorea de 10.000 celulas y calcular a partir de esto. Para forzar el calculo en mas valores, o todos los valores, se puede utilizar el parametro maxsamp

summary(info_B2, maxsamp = ncell(info_B2))
##         Landsat_recorte_B2
## Min.                     0
## 1st Qu.               8647
## Median                8922
## 3rd Qu.               9197
## Max.                 37616
## NA's               4336366

No se ven grandes diferencias en el resumen estadistico de maxsap, excepto con rasters muy grandes

Para visualizar estos datos en R se utiliza ggplot2, ya que se necesita convertir est a un dataframe. El paquete raster tiene una funcion de construccion para convertir a un dataframe ploteable.

Datos claves

Hay muchos paquetes aportados por el usuarios para R, cada uno con su propio conjunto de funciones. Las funciones en diferentes paquetes pueden tener el mismo nombre, aunque hagan cosas diferentes. La notacion package_name::function_name le dice a R que biblioteca de paquete desea que busque para una función en particular. En este caso, se requiere esta notacion, porque el conjunto de funciones basicas de R (la cual siempre se busca primero a menos que se especifique otra cosa) tambien incluye una funcion as.data.frame()

Cuando se vea la estructura de los datos, se vera un formato dataframe estandar

info_B2_df <- raster::as.data.frame(info_B2, xy = TRUE)
str(info_B2_df)
## 'data.frame':    9348878 obs. of  3 variables:
##  $ x                 : num  616800 616830 616860 616890 616920 ...
##  $ y                 : num  476040 476040 476040 476040 476040 ...
##  $ Landsat_recorte_B2: num  NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...

Se puede utilizar ggplot() para plotear estos datos. Se podria usar un conjunto de escala de colores amigables scale_fill_viridis_c

library(ggplot2)
ggplot() +
    geom_raster(data = info_B2_df , aes(x = x, y = y, fill = Landsat_recorte_B2)) +
    coord_equal() +
    scale_fill_gradientn(colours = terrain.colors(10))

Datos claves

El paquete complementario ggspatial puede tambien se puede usar. El raster es convertido en un dataframe antes de plotearse, pero el trabajo de convertir los datos espaciales en un formato amigable ggplot se hace en segundo plano.

Este mapa muestra la elevaion del sitio de estudio en el municipio de Villavicencio, sinembargo no se puede saber si son pies o metros porque la leyenda no muestra las unidades. Para saber las unidades es necesario ver los metadatos. Muchos de los metadatos que interesan son parte de CRS

Ver el sistema coordenadas de referencia (CRS) del raster en R

crs(info_B2)
## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0

Entender CRS en formato de proyeccion (Proj4)

El sistema de coordenadas de los datos esta dado por R en formato proj4. La cadena de proj4 contiene todos los elementos individuales del sistema de coordenadas que R u otro SIG necesitan. Cada elemento se especifica con un signo +, similar a como un archivo .csv es desglosado. Despues de cada + se define el elemento CRS. Ejemplo, poryeccion(porj =)y dato (datum =)

Projeccion UTM Proj4

La proyeccion UTM para Image_info_df es +proj=utm + zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84=0,0,0

  • proj=utm: Proyeccion es UTM la cual tiene varias zonas
  • zone=18: Zona 18
  • datum=WGS84: Se refiere a la referencia es 0,0 del sistema de coordenadas utilizado para la proyección
  • units=m: Las unidades para las coordenadas es en metros
  • ellps=WGS84: El elipsoide para los datos es WGS84 (como se calcula la redondez de la tierra)

Es importante tener en cuenta que la zona es exclusiva de la proyeccion UTM

Calcular los valores maximos y minimos del Raster

Es importante conoer los valores minimo y maximo del raster, que para el caso particular representa los max de elevacion en el sitio de estudio

minValue(info_B2)
## [1] 0
maxValue(info_B2)
## [1] 35668

Si los valores minimo maximo no han sido calculados, se pueden calcular usando la funcion setMinMax()

Image_info <- setMinMax(info_B2)

Bandas raster

La imagen trabajada es un raster de una sola banda. Esto quiere decir que unicamente se encuentra almacenado un datased en el raster

Un raster dataset puede contener una o mas bandas. Se puede utilizar la funcion raster() para importar rasters de una sola banda o multiples bandas.Se puede ver el numero de bandas de un raster utilizando la funcion nlayers()

nlayers(info_B2)
## [1] 1

Sin embargo, un dato raster tambien puede ser multibanda, es decir un archivo raster contiene datos de mas de una variable o periodo de tiempo en cada celda. Por defecto, la funcion raster() solo importa la primera banda en un raster independientemente de si tiene una o mas bandas.

Tratar datos perdidos

Los datos raster con frecuencia tienen NoDataValue, estos son valores asignados a los pixeles donde los datos estan perdidos o no fueron colectados.

Por defecto, la forma de un raster es rectangular. Si se tiene un conjunto de datos que no tienen una forma rectangular, algunos pixeles en el borde del raster tendran NoDataValues. Esto sucede cuando los datos fueron recolectados por un avion que solo sobrevoló una parte de una región definida

Si el raster ya tiene valores NA configurados, y se quiere estar seguros de donde estan, se puede utilizar ggspatial para trazarlos deliberadamente en un color particular. Esto puede ayudar a verificar la cobertura de un conjunto de datos.

Para resaltar los valores NA en ggplot, se debe modificar la capa scale_fil_*() para que contenga una instruccion de color para los valores de NA, como scale_fill_gradientn(na.value=‘deeppink’)

El valor que se usa de forma convencional para tomar nota de los datos faltantes (NoDataValue) varia segun el tipo de datos raster. Para los raster de coma flotante, la cifra -3.4e + 38 es un valor predeterminado comun y para los entereos, -9999 es comun.

Un archivo GeoTIFF tiene una etiqueta que dice que es NoDataValue o se puede encontrar esta informacion en los metadatos del raster. Si se almaceno un NoDataValue en la etiqueta GeoTIFF, cuando R abre el raster, asignara a cada instancia del valor a NA. Los valores NA seran ignorados por R

Valores malos de los datos Raster

Los valores malos son diferentes de los NoDataValue. Son valores que quedan fuera del rango aplicable del conjunto de datos

Ejemplos

  • El indice de vegetacion normalizado (NDVI), que es una medida de verdor, se encuentra en el rango de -1 a 1. Cualquier valor fuera de este rango se considera un valor malo.
  • Los datos de reflectancia de una imagen frecuentemente se encuentran en el rango de 0-1 o 0-10.000 dependiendo del escalamiento de los datos. POr lo tanto, es probable que un valor mayor a 1 o 10.000 sea causado por un error en la recopilacion o procesamiento de datos.

Encontrar los valores malos de los datos.

  • Algunas veces los metadatos de los rasters pueden dar el rango esperado de los valores del raster. Valores fuera de este rango pueden ser valores malos y se deben tener en cuenta en el momento de analizar los datos.
  • Alguans veces se necesita sentido comun e intuicion cientifica como examinar los datos - al igual que con los datos de campo para identificar valores cuastionables.
  • Plotear los datos con un resaltado apropiado puede ayudar a revelar patrones en valores incorrectos.

Crar un Histograma de los valores de un Raster

Para explorar la distribucion de los datos contenidos en el raster se utiliza la funcion geom_histogram() la cual genera el histograma, el cual ñuede ser util para identificar valores atipicos y datos erroneos de los datos raster.

ggplot() +
    geom_histogram(data = info_B2_df, aes(Landsat_recorte_B2))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 4336366 rows containing non-finite values (stat_bin).

El mensaje de advertencia es causado por ggplot debido a la configuracion que impone que haya 30 contenedores para los datos. Este valor de contenedores se puede definir utilizando la funcion geom_histogram()

ggplot() +
    geom_histogram(data = info_B2_df, aes(Landsat_recorte_B2), bins = 20)
## Warning: Removed 4336366 rows containing non-finite values (stat_bin).

La forma de los dos histogramas es similar.