Introduccion a datos raster en R
Lorena Orellana
Introducción
En este tutorial se explicará las funciones basicas para trabajar datos raster. Se utilizará las librerias: raster, rgdal y ggplot2
Se utilizará una banda en el cual se calculo el ndvi para la ciudad de cuenca
Cargar librerias
library(raster)## Loading required package: splibrary(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: D:/Lore/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: D:/Lore/Documents/R/win-library/3.4/rgdal/proj
## Linking to sp version: 1.2-5library(ggplot2)## Warning: package 'ggplot2' was built under R version 3.4.4Mirar los atributos de un raster
Para obtener información de sobre un raster antes de que este sea importado utilizamos la funcion GDALinfo()
#Colocar el directorio de trabajo en donde se encuentra las bandas de nuestra imagen satelite
setwd("D:\\maestria\\percepcion remota avanzada\\ejerciciosR")
# Observamos las caracteristicas de la Banda 1 de la imagen
GDALinfo ("ndvi_cuenca.TIF")## rows 429
## columns 671
## bands 1
## lower left origin.x 714225
## lower left origin.y 9675115
## res.x 30
## res.y 30
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## file ndvi_cuenca.TIF
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float32 TRUE -3.4e+38 3 671
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 -0.4412046 0.8200465 0.2993605 0.1564616
## Metadata:
## AREA_OR_POINT=Area# para gaurdar la informacion en R
info_isoterma <- capture.output(GDALinfo("ndvi_cuenca.TIF"))
info_isoterma## [1] "rows 429 "
## [2] "columns 671 "
## [3] "bands 1 "
## [4] "lower left origin.x 714225 "
## [5] "lower left origin.y 9675115 "
## [6] "res.x 30 "
## [7] "res.y 30 "
## [8] "ysign -1 "
## [9] "oblique.x 0 "
## [10] "oblique.y 0 "
## [11] "driver GTiff "
## [12] "projection +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs "
## [13] "file ndvi_cuenca.TIF "
## [14] "apparent band summary:"
## [15] " GDType hasNoDataValue NoDataValue blockSize1 blockSize2"
## [16] "1 Float32 TRUE -3.4e+38 3 671"
## [17] "apparent band statistics:"
## [18] " Bmin Bmax Bmean Bsd"
## [19] "1 -0.4412046 0.8200465 0.2993605 0.1564616"
## [20] "Metadata:"
## [21] "AREA_OR_POINT=Area "Abrir un Raster en R
Con la funcion raster() podemos abrir un archivo .tif, y nos muestra su estructura, asignamos la banda a una variable
ndvi <- raster("ndvi_cuenca.TIF")Con la funcion summary nos reporta las estadisticas: MIN, MAX, MEDIANA, LOS CUARTILES Y LOS NA’s. EN este caso el resumen es una estimacion basada en una muestra de pixeles
summary(ndvi)## Warning in sampleInt(length(x), size): size changed to n because it cannot
## be larger than n when replace is FALSE## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (34.74% of all cells)## ndvi_cuenca
## Min. -0.4412046
## 1st Qu. 0.1685399
## Median 0.2825156
## 3rd Qu. 0.4188275
## Max. 0.8200465
## NA's 0.0000000Para obtener el resumen con todos los valores de la banda utilizamos la funcion maxsamp
summary(ndvi, maxsamp = ncell(ndvi))## ndvi_cuenca
## Min. -4.412046e-01
## 1st Qu. 1.685399e-01
## Median 2.825156e-01
## 3rd Qu. 4.188275e-01
## Max. 8.200465e-01
## NA's 2.064630e+05Para visualizar los datos raster podemos utilizar la funcion plot() o ggplot(). Cuando utilizamos la libreria ggplot2 necesitamos convertir los datos de tipo raster a data.frame. utilizamos las escala de color scale_fill_viridis_c
#utilizando la funcion plot
plot(ndvi)
#utilizando la libreria ggplot2.
df_ndvi <- as.data.frame(ndvi, xy=TRUE)
str(df_ndvi)## 'data.frame': 287859 obs. of 3 variables:
## $ x : num 714240 714270 714300 714330 714360 ...
## $ y : num 9687970 9687970 9687970 9687970 9687970 ...
## $ ndvi_cuenca: num NA NA NA NA NA NA NA NA NA NA ...ggplot() +
geom_raster(data = df_ndvi , aes(x = x, y = y, fill = ndvi_cuenca)) +
coord_equal()+
scale_fill_gradientn(colours = rev(terrain.colors(10)))
Mirar el sistema de referencia de coordenadas en R
Podemos conocer el sistema de referencia de un raster utilizando la funcion crs()
crs(ndvi)## CRS arguments:
## +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## +ellps=WGS84 +towgs84=0,0,0Calcular los valores maximo y minimo
Para calcular estos valores utilizamos las funciones minValue y maxValue
minValue(ndvi)## [1] -0.4412046maxValue(ndvi)## [1] 0.8200465Bandas
Para conocer cuantas bandas tiene el raster podemos utilizar la funcion nlayers(). En este caso estamos trabajando con una sola banda
nlayers(ndvi)## [1] 1Datos perdidos
NoDataValue esta asignado a pixeles donde faltan datos Aveces cuando los rasters no tienen una forma rectangular, los pixeles en el borde tendran “valores” NoDataVAlue
En la imagen a continuación, los píxeles que son negros tienen valores NoData. La cámara no recopiló datos en estas áreas.
Crear histogramas
Utilizar histrogramas pueden ser muy util para conocer la distribucion de nuestros datos, podemos utilizar dentro de la fincion ggplot el argumento bins= para determinar el numero de separaciones de nuestro histograma
ggplot() +
geom_histogram(data = df_ndvi, aes(ndvi_cuenca), bins = 40)## Warning: Removed 206463 rows containing non-finite values (stat_bin).