Análisis geoespacial y geoestadística con R

Una selección de códigos de R sobre manipulación de datos geoespaciales y geoestadística. Lo presenté en el Santo Domingo useR Group.

El guión del punto “0” no se desarrolla en sentencias, sino que forma parte del contenido teórico explicado en la sesión del grupo. Los demás se centran en sencillas operaciones que muestran las capacidades de R en la visualización y gestión de datos. Se trata de una pequeña muestra de las capacidades de R. Espero sea de utilidad.

#GUIÓN DE ANÁLISIS ESPACIAL/GEOESPACIAL CON R
##0) CONCEPTOS: ANÁLISIS ESPACIAL/GEOESPACIAL, GEOESTADÍSTICA (SEMIVARIOGRAMA, AUTOCORRELACIÓN ESPACIAL)
##1) PAQUETES ESPACIALES. OBJETOS ESPACIALES
##2) CARGAR MAPAS EN R. DATOS ASOCIADOS A PAQUETES Y DATOS PROPIOS
##2.1) MAPAS DENTRO DE R
##2.2) GENERAR MAPAS PARA GOOGLEEARTH
##2.3) CREAR OBJETOS ESPACIALES A PARTIR DE COORDENADAS Y GENERAR SALIDAS PARA GOOGLEMAPS
##3) UNIONES
##3.1) UNIÓN POR ATRIBUTOS. SALUD A NIVEL DE PROVINCIAS
##4) GEOESTADÍSTICA. EJEMPLO CON LA LLUVIA DE RD PROMEDIADA PARA EL PERIODO 1979-2014

#1) PAQUETES
##ESPACIALES
library(raster) #PARA GESTIONAR IMÁGENES RASTER. TAMBIÉN APORTA FUENTES A NIVEL DE PAÍS O DE TERRITORIOS SOBRE DIVISIÓN TERRITORIAL, CLIMA Y TOPOGRAFÍA
library(sp) #ANÁLISIS ESPACIAL
library(maps) #PARA GENERAR MAPAS, SOBRE TODO DE ESTADOS UNIDOS, AUNQUE TAMBIÉN MUNDIALES
library(mapdata) #BASE DE DATOS DE MAPAS EXTRA, QUE AÑADE ALGUNOS TEMAS
library(maptools) #PARA MANEJAR SHAPEFILES
library(rgdal) #PARA MANEJAR ARCHIVOS DE LA BIBLIOTECA GDAL (GEOSPATIAL DATA ABSTRACTION LIBRARY)
library(RgoogleMaps) #PARA GENERAR, EN VENTANA GRÁFICA DE R, MAPAS DE GOOGLE
library(plotGoogleMaps) #PARA GENERAR ARCHIVOS html CON LOS QUE SE VISUALIZAN OBJETOS ESPACIALES EN GOOGLEMAPS
library(plotKML) #VISUALIZACIÓN DE OBJETOS ESPACIALES Y ESPACIO-TEMPORAL DE OBJETOS EN GOOGLE EARTH
library(rworldmap) #PARA GENERAR MAPAS
library(gstat) #MODELIZACIÓN GEOESTADÍSTICA ESPACIAL Y ESPACIO-TEMPORAL, PREDICCIÓN Y SIMULACIÓN
library(geoR) #ANÁLISIS DE DATOS GEOESTADÍSTICOS
library(dismo) #PARA MODELIZAR DISTRIBUCIÓN DE ESPECIES. SE UTILIZARÁ SÓLO EL COMANDO gbif, PERO TIENE MUCHAS FUNCIONES INTERESANTES

##ÚTILES PARA PREPARACIÓN Y REPRESENTACIÓN DE DATOS. NO USADOS EN ESTE SCRIPT ESPECÍFICO, PERO MUCHOS GRÁFICOS Y TABLAS FUERON HECHOS CON ESTOS
#library(dplyr) #MANIPULACIÓN DE DATOS EN TABLAS, TABLAS RESUMEN
#library(ggplot2) #GRÁFICOS ESTILIZADOS

#2) CREAR/CARGAR OBJETOS ESPACIALES EN R

##2.2) DENTRO DE R
###DEL PAQUETE MAPS
dev.new()
map()	#MAPAMUNDI DE BAJA RESOLUCIÓN
dev.new()
map('worldHires','Dominican Republic',col='black',fill=T)
box() #AÑADE BORDE AL MAPA ACTUAL

###DEL PAQUETE Rgooglemaps. MOSTRANDO MAPA DE GOOGLE DE LA ISLA DENTRO DE R
PlotOnStaticMap(GetMap(center=c(lat=19,lon=-71.5), zoom=7,maptype="satellite"))

###DEL PAQUETE RASTER
####DATOS ASOCIADOS AL PAQUETE
dop<-getData('GADM', country='DOM', level=1)
dom<-getData('GADM', country='DOM', level=2)
#getData('SRTM',lat=19,lon=-71) #DEM POR GRANDES TERRITORIOS. TIENE UN INCONVENIENTE: OBLIGA A DESCARGAR ARCHIVOS MUY GRANDES. NO ESTÁN POR PAÍS, POR TILES. RD SE CONSTRUYE CONCATENANDO DOS TILES. ARCHIVOS GRANDES
dev.new()
sp::spplot(dop)
dev.new()
sp::spplot(dom)
dev.new()
sp::spplot(dop['ID_1'],col.regions = rainbow(32, start = 0, end = 1))
dop$NAME_1=factor(dop$NAME_1)
sp::spplot(dop['NAME_1'],col.regions = rainbow(32, start = 0, end = 1))
dev.new()
sp::spplot(dom['ID_1'])

####DATOS DESCARGADOS VÍA http Y TRABAJADOS EN LOCAL
setwd(tempdir()) #FIJAR EL DIRECTORIO DE TRABAJO
getwd() #MUESTRA DÓNDE ESTÁ EL DIRECTORIO DE TRABAJO
temp <- tempfile() #ASIGNAR EL NOMBRE DE UN ARCHIVO TEMPORAL AL OBJETO temp
download.file("http://geografiafisica.org/r/maestria_geom/practica_01.zip",temp) #DESCARGA EL ARCHIVO ZIP Y LO NOMBRA COMO TEMPORAL
unzip(temp) #DESCOMPRIME EN EL DIRECTORIO DE TRABAJO
unlink(temp) #DESVINCULA DEL ARCHIVO TEMPORAL Y LO BORRA
dem<-raster(readGDAL("demesp.tif")) #LEE EL TIF DEL DEM DE LA ESPAÑOLA CON LA FUNCIÓN readGDAL DEL PAQUETE rgdal, EL CUAL SE ENCONTRABA EN EL COMPRIMIDO ZIP, Y LO ASIGNA AL OBJETO dem
dem #MUESTRA LAS CARACTERÍSTICAS DEL raster DENOMINADO dem
dev.new() #CREA UNA VENTANA GRÁFICA NUEVA
plot(dem) #VISUALIZA EL RASTER dem EN UN MAPA
raster::zoom(dem) #PARA HACER ZOOM INTERACTIVO MEDIANTE DOS PUNTOS DE UN RECTÁNGULO IMAGINARIO. HACER CLIC PRIMERO EN UNO DE LOS PUNTOS DEL MAPA Y LUEGO HACER CLIC EN EL SEGUNDO; EL ÁREA CUBIERTA POR EL RECTÁNGULO IMAGINARIO SERÁ LA MOSTRADA EN EL PLOT. COMENTADA PARA EVITAR DETENCIÓN DE SCRIPT
sneyba<-readOGR('sneyba.kml','sneyba.kml') #LEE UN ARCHIVO KML CON LÍMITES ARBITRARIOS DE LA SIERRA DE NEYBA Y LO ASIGNA AL OBJETO sneyba, EL CUAL SE ENCONTRABA EN EL COMPRIMIDO ZIP
sneyba #MUESTRA LAS CARACTERÍSTICAS DEL OBJETO ESPACIAL DENOMINADO sneyba
sneyba<-spTransform(sneyba,CRS("+init=epsg:32619")) #TRANSFORMA EL OBJETO ESPACIAL sp DEL SISTEMA DE REFERENCIA ESPACIAL LATLONG/WGS84 AL UTM/WGS84
dev.new()
plot(dem) #VISUALIZA EL RASTER dem EN UN MAPA
plot(sneyba,add=T) #VISUALIZAR EL OBJETO ESPACIAL sneyba, AÑADIÉNDOLO AL MAPA ANTERIOR
rsneyba<-raster::mask(crop(dem,extent(sneyba)),sneyba)
dev.new()
plot(rsneyba)
head(getValues(rsneyba)) #RECUPERA LOS PRIMEROS 6 VALORES
dev.new()
hist(getValues(rsneyba)) #HISTOGRAMA DE LOS VALORES DE ALTURA
dev.new()
hist(rsneyba[getValues(rsneyba)>100]) #HISTOGRAMA DE LOS VALORES DE ALTURA DE MÁS DE 100 METROS

###2.2) GENERAR MAPAS PARA GOOGLEEARTH
plotKML(rsneyba)
####LO MISMO, PERO GUARDANDO ARCHIVO KML, USANDO EL PAQUETE RASTER
####rsneyball<-projectRaster(rsneyba, crs="+proj=longlat +datum=WGS84", method='ngb')
####KML(rsneyba,file='rsneyba.kml')
plotKML(dop['ID_1'])
plotKML(dop['NAME_1'],shape="")
plotKML(dop['NAME_1'],shape='http://maps.google.com/mapfiles/kml/shapes/placemark_circle.png')

###2.3) CREAR OBJETOS ESPACIALES A PARTIR DE COORDENADAS Y GENERAR UNA SALIDA EN FORMA DE MAPA PARA GOOGLEMAPS

####2.3.1 MUESTRAS DE "CALLAOS" DE RÍO
cd.env<-read.csv('http://geografiafisica.org/geo112/practica_04/integrados_env.csv',head=T) #DESCARGA Y CONVIERTE A UNA TABLA (data.frame) EL ARCHIVO CONTENIENDO LA INFORMACIÓN AMBIENTAL DE LOS MUESTREOS
cd.env #MUESTRA EL OBJETO cd.env
cd.env.ok<-subset(cd.env, cd.env$coord_validadas=='sí') #CREA UN OBJETO CON AQUELLOS REGISTROS DE d.env CUYAS COORDENADAS ESTÁN VALIDADAS
str(cd.env.ok) #MUESTRA LA ESTRUCTURA DEL OBJETO d.env.ok, DONDE LOS CAMPOS x_inicial Y y_inicial APARECEN COMO FACTORIALES, DADO QUE AL SER IMPORTADOS CONTENÍAN TEXTO
cd.env.ok$x_inicial<-as.numeric(as.character(cd.env.ok$x_inicial)) #CONVIERTE A NUMÉRICO EL CAMPO x_inicial, PARA QUE SE ASUMAN ADECUADAMENTE LAS COORDENADAS
cd.env.ok$y_inicial<-as.numeric(as.character(cd.env.ok$y_inicial)) #CONVIERTE A NUMÉRICO EL CAMPO y_inicial, PARA QUE SE ASUMAN ADECUADAMENTE LAS COORDENADAS
str(cd.env.ok) #MUESTRA LA ESTRUCTURA DEL OBJETO d.env.ok, DONDE LOS CAMPOS x_inicial Y y_inicial APARECEN AHORA COMO NUMÉRICOS
coordinates(cd.env.ok)<-~x_inicial+y_inicial #INDICA AL PROGRAMA QUE LAS COORDENADAS DEL OBJETO t.env.ok ESTÁN EN LOS CAMPOS x_inicial E y_inicial. ESTO CONVERTIRÁ A d.env.ok EN UN OBJETO ESPACIAL
proj4string(cd.env.ok) <- CRS("+init=epsg:32619") #LE ASIGNA EL SISTEMA DE REFERENCIA DE COORDENADAS CORRESPONDIENTE AL OBJETO
plotGoogleMaps(cd.env.ok,iconMarker='',zcol="codigo","cantometria.html") #GENERA UN ARCHIVO .html EN LA CARPETA DE TRABAJO, LA CUAL SE SABE ESCRIBIENDO getwd(). EL PROGRAMA LANZARÁ EL NAVEGADOR DE MANERA AUTOMÁTICA, Y CARGARÁ LOS PUNTOS. NOTA: SI R NO TIENE PERMISOS DE ESCRITURA EN EL DIRECTORIO DE TRABAJO, NO PODRÁ CREAR EL ARCHIVO Y DARÁ ERROR. EN ESE CASO, SE RECOMIENDA FIJAR UNA CARPETA DE TRABAJO EN UN DIRECTORIO CON TODOS LOS PRIVILEGIOS, USANDO EL COMANDO setwd(CARPETA). POR EJEMPLO, setwd("C:/"), O setwd("D:/"), O setwd("C:/Users/Public")

####2.3.2 OCURRENCIAS DE UN GÉNERO DE PLANTAS: Lantana
lantana<-gbif('Lantana','*',args='originisocountrycode=DO',geo=T) #GENERA EL OBJETO lantana, QUE CONTIENE TODAS LAS OCURRENCIAS, DENTRO DE LA BASE DE DATOS GBIF, DE ESPECIES CON COORDENADAS DEL GÉNERO Lantana, REFERIDAS A REPÚBLICA DOMINICANA
lantana<-lantana[,c('species','catalogNumber','lat','lon','alt')] #EXTRAE SÓLO LAS COLUMNAS QUE INTERESA CONSERVAR
str(lantana) #MUESTRA LA ESTRUCTURA DEL OBJETO lantana
lantana<-subset(lantana,lon<(-68)&lat>17) #FILTRA (CONSERVA) LAS FILAS QUE TIENEN LONGITUDES INFERIORES A -68°W Y SUPERIORES A 17°N (RD)
coordinates(lantana) <- c("lon", "lat")  #ASIGNA COORDENADAS AL OBJETO DESDE LAS COLUMNAS "lat" "lon", LO CUAL CONVIERTE A LANTANA EN OBJETO ESPACIAL
proj4string(lantana) <- CRS("+init=epsg:4312") #ASIGNA LA PROYECCIÓN QUE CORRESPONDE AL OBJETO
plotGoogleMaps(lantana,iconMarker='',zcol="catalogNumber","lantana.html") #GENERA EL html, LO CUAL ABRIRÁ UNA VENTANA DEL NAVEGADOR POR DEFECTO

##3) UNIONES
##3.1) UNIÓN POR ATRIBUTOS. SALUD A NIVEL DE PROVINCIAS. FUENTE DE LOS DATOS DE SALUD: http://www.bvs.org.do/bvs/htdocs//local/File/indicadores--basicos-de-salud-2013.pdf
ds<-read.csv('http://geografiafisica.org/r/datos_salud.csv')
dopds<-merge(dop,ds)
sp::spplot(dopds['NAME_1'],col.regions = rainbow(32, start = 0, end = 1))
sp::spplot(dopds['pct.poblacion.urbana.2012'],col.regions=heat.colors(100))
sp::spplot(dopds['pct.poblacion.urbana.2012'],col.regions=rev(heat.colors(100)))
sp::spplot(dopds['pct.poblacion.urbana.2012'],col.regions=colorRampPalette(c("white", "black"))(20))
cor(ds[,sapply(ds,is.numeric)])

##4) GEOESTADÍSTICA. EJEMPLO CON LA LLUVIA DE RD PROMEDIADA PARA EL PERIODO 1970-2014
d.env.t<-read.csv('http://geografiafisica.org/r/pm1979_2015.csv')

###PRUEBAS DE NORMALIDAD, HISTOGRAMAS
dev.new()
hist(d.env.t$Pm)
qqnorm(d.env.t$Pm)
shapiro.test(d.env.t$Pm)

###CONVERSIÓN A SP Y EXPLORACIÓN
coordinates(d.env.t)<-~x+y
proj4string(d.env.t)<-CRS("+init=epsg:32619")
dev.new()
sp::spplot(d.env.t['Pm'])
plotKML(d.env.t['Pm'])

###VARIOGRAMA
plot(variogram(Pm~x+y, d.env.t, alpha=seq(0,360,by=10),cloud=T))
dev.new();plot(variogram(Pm~x+y, d.env.t, alpha=c(110)))
dev.new();plot(variogram(Pm~x+y, d.env.t, alpha=c(280)))
rango <- sqrt(diff(d.env.t@bbox['x',])^2 + diff(d.env.t@bbox['y',])^2)/4 
rango
Pm.sph <- vgm(nugget=0,model='Sph',range=rango)
Pm.sv <- variogram(Pm~x+y, d.env.t, alpha=c(280))
Pm.sph.f <- fit.variogram(Pm.sv,model=Pm.sph)
plot(Pm.sv,Pm.sph.f,pch="+", pl=TRUE,col="black", main="Pm")

###MALLA
grd <- expand.grid(x=seq(from=min(d.env.t@bbox['x',]), to=max(d.env.t@bbox['x',]), by=1000),y=seq(from=min(d.env.t@bbox['y',]), to=max(d.env.t@bbox['y',]), by=1000)) 
str(grd) #MUESTRA LA ESTRUCTURA DEL OBJETO grd. DOS COLUMNAS, x E y, CONTIENEN LAS COORDENADAS. EN ESTE CASO, SE TRATA DE UN data.frame
coordinates(grd) <- ~x+y #CONVIERTE A grd EN UN OBJETO ESPACIAL (SpatialPoints) INDICÁNDOLE QUÉ COLUMNAS CONTIENEN LAS COORDENADAS x E y
spplot(grd) #VISUALIZA EL OBJETO grd. COMO NO TIENE NINGUNA VARIABLE DE DATOS, REPRESENTA LOS VALORES DE LAS COORDENADAS
proj4string(grd) <- CRS("+init=epsg:32619") #ASIGNANDO EL MISMO SISTEMA DE REFERENCIA DE COORDENADAS DEL OBJETO d

###INTERPOLACIÓN
OKPm <- krige(id="Pm", formula=Pm~1,d.env.t, newdata=grd, model = Pm.sph.f)
str(OKPm) #ESTRUCTURA DE LA INTERPOLACIÓN
gridded(OKPm)<-TRUE #CUADRICULANDO LA INTERPOLACIÓN
str(OKPm) #NUEVA ESTRUCTURA DE LA INTERPOLACIÓN

###VISTA EN GE
plotKML(d.env.t, colour_scale=rep("#FFFF00", 2), points_names=d.env.t$estacion) #DESPLAZADAS. NO SE CORRESPONDEN CON LOS LUGARES PRECISOS DE LAS ESTACIONES. LA ESTACIÓN DE BAYAGUANA COINCIDÍA CON LA DE RANCHO ARRIBA
plotKML(OKPm["Pm.pred"], colour_scale = SAGA_pal[[1]]) 
writeGDAL(OKPm["Pm.pred"],"Pm.pred.tif", drivername="GTiff")

#ALGUNOS RECURSOS

##BIBLIOGRAFÍA
###Hengl, T. (2009): "A Practical Guide to Geostatistical Mapping"
###Bivand, R.; Pebesma, E.; Gómez-Rubio, Virgilio (2013). Applied Spatial Data Analysis with R". Springer

##DATOS DE SALUD
###Ministerio de Salud Pública (MSP); Organización Panamericana de la Salud (OPS); Organización Mundial de la Salud (2013): "Indicadores básicos de salud 2013". URL: http://www.bvs.org.do/bvs/htdocs//local/File/indicadores--basicos-de-salud-2013.pdf

##DATOS ONE (ÚTILES PARA UNIONES POR ATRIBUTOS)
###SHAPEFILES (DESCOMPRIMIR, FIJAR RUTA DE TRABAJO): http://www.one.gov.do/cartografia/276/informaciones-cartograficas
###REDATAM (DESCARGAR TABLAS EN FORMATO EXCEL. ESTOS ARCHIVOS, CUANDO SE ABREN EN EXCEL, EL CAMPO "codigo" PIERDE EL PRIMER DIGITO CUANDO ÉSTE ES "0"): http://redatam.one.gob.do/cgibin/RpWebEngine.exe/PortalAction?&MODE=MAIN&BASE=CPV2010&MAIN=WebServerMain.inl

##DIRECCIONES WEB:
###http://geostat-course.org/node
###https://sites.google.com/site/rodriguezsanchezf/news/usingrasagis
###https://www.youtube.com/watch?v=Ha06LxeI4Z8
###TASK VIEW DE ANÁLISIS ESPACIAL EN R (AUTOR: BIVAND): http://cran.r-project.org/web/views/Spatial.html
###http://www.molecularecologist.com/2012/09/making-maps-with-r/
###http://www.r-bloggers.com/create-maps-with-maptools-r-package/
###https://www.nceas.ucsb.edu/scicomp/usecases/CreateMapsWithRGraphics
###https://procomun.wordpress.com/2012/02/18/maps_with_r_1/
###https://www.youtube.com/watch?v=mMaMmaTfsQE
###https://www.youtube.com/watch?v=YYh6D8x8qng
###https://www.youtube.com/watch?v=e2XkA-_OovM #Mobilize, privado
###http://r-video-tutorial.blogspot.com/2014/02/displaying-spatial-sensor-data-from.html

Created by Pretty R at inside-R.org

Adicionalmente, recomiendo este trabajo de Tomislav Hengl. Es una buena síntesis sobre procedimientos básicos en geoestadística.

Dr. José Ramón Martínez Batlle (Ph.D)

Comments

Un pensamiento en “Análisis geoespacial y geoestadística con R

  1. Pingback: Análisis geoespacial y geoestadística con R ¡actualizado! | Geografía Física – República Dominicana – Dr. José Ramón Martínez Batlle

Deja un comentario

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *