Análisis geoespacial y geoestadística con R

geost.R

#GUION DE ANALISIS ESPACIAL/GEOESPACIAL CON R
##0) CONCEPTOS: ANALISIS ESPACIAL/GEOESPACIAL, GEOESTADISTICA (SEMIVARIOGRAMA, AUTOCORRELACION 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) UNION POR ATRIBUTOS. SALUD A NIVEL DE PROVINCIAS
##4) GEOESTADISTICA. EJEMPLO CON LA LLUVIA DE RD PROMEDIADA PARA EL PERIODO 1979-2014


#1) PAQUETES
##ESPACIALES
library(raster) #PARA GESTIONAR IMAGENES RASTER. TAMBIEN APORTA FUENTES A NIVEL DE PAIS O DE TERRITORIOS SOBRE DIVISION TERRITORIAL, CLIMA Y TOPOGRAFIA
library(sp) #ANALISIS ESPACIAL
library(maps) #PARA GENERAR MAPAS, SOBRE TODO DE ESTADOS UNIDOS, AUNQUE TAMBIEN MUNDIALES
library(mapdata) #BASE DE DATOS DE MAPAS EXTRA, QUE ANADE 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 GRAFICA DE R, MAPAS DE GOOGLE
library(plotGoogleMaps) #PARA GENERAR ARCHIVOS html CON LOS QUE SE VISUALIZAN OBJETOS ESPACIALES EN GOOGLEMAPS
library(plotKML) #VISUALIZACION DE OBJETOS ESPACIALES Y ESPACIO-TEMPORAL DE OBJETOS EN GOOGLE EARTH
library(rworldmap) #PARA GENERAR MAPAS
library(gstat) #MODELIZACION GEOESTADISTICA ESPACIAL Y ESPACIO-TEMPORAL, PREDICCION Y SIMULACION
library(geoR) #ANALISIS DE DATOS GEOESTADISTICOS
library(dismo) #PARA MODELIZAR DISTRIBUCION DE ESPECIES. SE UTILIZARA SOLO EL COMANDO gbif, PERO TIENE MUCHAS FUNCIONES INTERESANTES

##UTILES PARA PREPARACION Y REPRESENTACION DE DATOS. NO USADOS EN ESTE SCRIPT ESPECIFICO, PERO MUCHOS GRAFICOS Y TABLAS FUERON HECHOS CON ESTOS
#library(dplyr) #MANIPULACION DE DATOS EN TABLAS, TABLAS RESUMEN
#library(ggplot2) #GRAFICOS ESTILIZADOS

setwd(tempdir()) #FIJAR EL DIRECTORIO DE TRABAJO

#2) CREAR/CARGAR OBJETOS ESPACIALES EN R


##2.2) DENTRO DE R
###DEL PAQUETE MAPS
dev.new()
dev.new()
map()	#MAPAMUNDI DE BAJA RESOLUCION
dev.new()
map('worldHires','Dominican Republic',col='black',fill=T)
map('worldHires',c('Dominican Republic', 'Haiti'),col='black',fill=T)
box() #ANADE BORDE AL MAPA ACTUAL

###DEL PAQUETE Rgooglemaps. MOSTRANDO MAPA DE GOOGLE DE LA ISLA DENTRO DE R
dev.new(); 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 ESTAN POR PAIS, 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)
dev.new()
sp::spplot(dop['NAME_1'],col.regions = rainbow(32, start = 0, end = 1))
dev.new()
sp::spplot(dop['NAME_1'],col.regions = sample(rainbow(32),32))
dev.new()
sp::spplot(dom['ID_1'])

####DATOS DESCARGADOS VIA http Y TRABAJADOS EN LOCAL
getwd() #MUESTRA DONDE ESTA 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 ESPANOLA CON LA FUNCION readGDAL DEL PAQUETE rgdal, EL CUAL SE ENCONTRABA EN EL COMPRIMIDO ZIP, Y LO ASIGNA AL OBJETO dem
dem #MUESTRA LAS CARACTERISTICAS DEL raster DENOMINADO dem
dev.new() #CREA UNA VENTANA GRAFICA NUEVA
plot(dem) #VISUALIZA EL RASTER dem EN UN MAPA
raster::zoom(dem) #PARA HACER ZOOM INTERACTIVO MEDIANTE DOS PUNTOS DE UN RECTANGULO IMAGINARIO. HACER CLIC PRIMERO EN UNO DE LOS PUNTOS DEL MAPA Y LUEGO HACER CLIC EN EL SEGUNDO; EL AREA CUBIERTA POR EL RECTANGULO IMAGINARIO SERA LA MOSTRADA EN EL PLOT. COMENTADA PARA EVITAR DETENCION DE SCRIPT
sneyba<-readOGR(dsn='sneyba.kml',layer='sneyba') #LEE UN ARCHIVO KML CON LIMITES ARBITRARIOS DE LA SIERRA DE NEYBA Y LO ASIGNA AL OBJETO sneyba, EL CUAL SE ENCONTRABA EN EL COMPRIMIDO ZIP
sneyba #MUESTRA LAS CARACTERISTICAS 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, ANADIENDOLO AL MAPA ANTERIOR
rsneyba<-raster::mask(crop(dem,extent(sneyba)),sneyba)
dev.new()
plot(rsneyba)
valsneyba <- getValues(rsneyba)
head(valsneyba) #RECUPERA LOS PRIMEROS 6 VALORES
mean(valsneyba, na.rm=T)
dev.new()
hist(valsneyba) #HISTOGRAMA DE LOS VALORES DE ALTURA
dev.new()
hist(rsneyba[valsneyba>100]) #HISTOGRAMA DE LOS VALORES DE ALTURA DE MAS DE 100 METROS


###2.2) GENERAR MAPAS PARA GOOGLEEARTH
plotKML(dop['ID_1']) #SE CREA UN KML EN EL DIRECTORIO DE TRABAJO (PARA CONOCER EL DIRECTORIO DE TRABAJO: "getwd()" SIN COMILLAS). GOOGLE EARTH SE ABRE DE FORMA AUTOMATICA. EN CASO CONTRARIO, LOCALIZAR EL KML EN EL DIRECTORIO DE TRABAJO Y ABRIRLO MANUALMENTE. SI SE ESPECIFICA EL ARGUMENTO open.kml=F SE EVITA LA APERTURA AUTOMATICA DE GOOGLEEARTH
plotKML(dop['NAME_1'],shape="", open.kml = F)
plotKML(dop['NAME_1'],shape='http://maps.google.com/mapfiles/kml/shapes/placemark_circle.png', open.kml = F)


###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 RIO
cd.env<-read.csv('http://geografiafisica.org/geo112/practica_04/integrados_env.csv',encoding = 'latin1', head=T) #DESCARGA Y CONVIERTE A UNA TABLA (data.frame) EL ARCHIVO CONTENIENDO LA INFORMACION AMBIENTAL DE LOS MUESTREOS
cd.env #MUESTRA EL OBJETO cd.env
cd.env.ok<-cd.env #CREA UN OBJETO CON AQUELLOS REGISTROS DE d.env CUYAS COORDENADAS ESTAN 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 CONTENIAN TEXTO
cd.env.ok$x_inicial<-as.numeric(as.character(cd.env.ok$x_inicial)) #CONVIERTE A NUMERICO 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 NUMERICO 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 NUMERICOS
coordinates(cd.env.ok)<- ~x_inicial+y_inicial #INDICA AL PROGRAMA QUE LAS COORDENADAS DEL OBJETO t.env.ok ESTAN EN LOS CAMPOS x_inicial E y_inicial. ESTO CONVERTIRA 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 LANZARA EL NAVEGADOR DE MANERA AUTOMATICA, Y CARGARA LOS PUNTOS. NOTA: SI R NO TIENE PERMISOS DE ESCRITURA EN EL DIRECTORIO DE TRABAJO, NO PODRA CREAR EL ARCHIVO Y DARA 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")

##3) UNIONES
##3.1) UNION 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_salud2.csv', check.names = T, encoding = 'latin1')
dopds<-merge(dop,ds)
dev.new()
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), main = 'Porcentaje de poblacion urbana, 2012')
sp::spplot(dopds['pct.poblacion.urbana.2012'],col.regions=rev(heat.colors(100)), main = 'Porcentaje de poblacion urbana, 2012')
sp::spplot(dopds['pct.poblacion.urbana.2012'],col.regions=colorRampPalette(c("white", "black"))(20), main = 'Porcentaje de poblacion urbana, 2012')
correl <- cor(ds[,sapply(ds,is.numeric)])
correl
View(correl)

##4) GEOESTADISTICA. 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)

###CONVERSION A SP Y EXPLORACION
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'], open.kml = F)

###VARIOGRAMA
dev.new();plot(variogram(Pm~x+y, d.env.t, alpha=seq(0,360,by=10))) #PANEL DE GRAFICOS DE SEMIVARIOGRAMA EN DISTINTAS DIRECCIONES. APARENTEMENTE, LA DIRECCION 110° ES LA QUE MUESTRA UN ALCANCE REAL DE MESETA
dev.new();plot(variogram(Pm~x+y, d.env.t, alpha=c(110), cloud = T)) #cloud = T HACE QUE SE CALCULE EL SEMIVARIOGRAMA EN NUBE
dev.new();plot(variogram(Pm~x+y, d.env.t, alpha=c(110), cloud = F)) #SIN NUBE
dev.new();plot(variogram(Pm~x+y, d.env.t, alpha=c(280), cloud = T)) #UNA DIRECCION SIMILAR EN ANGULO COMPLEMENTARIO
rango <- sqrt(diff(d.env.t@bbox['x',])^2 + diff(d.env.t@bbox['y',])^2)/4  #DIAGONAL DEL ENCUADRE DIVIDIDO ENTRE 4. ESTA RELACION SUELE USARSE COMO REGLA DE ORO PARA ESTABLECER UN RANGO INICIAL DEL SEMIVARIOGRAMA MODELO
rango
Pm.sph <- vgm(nugget=0,model='Sph',range=rango) #SEMIVARIOGRAMA MODELO, EN ESTE CASO EL ELEGIDO ES ESFERICO. SE UTILIZA EL RANGO CALCULADO (DIAGONAL ENTRE CUATRO)
Pm.sv <- variogram(Pm~x+y, d.env.t, alpha=c(110)) #SEMIVARIOGRAMA MUESTRAL, PERO ASIGNADO A UN OBJETO
Pm.sph.f <- fit.variogram(Pm.sv,model=Pm.sph) #SEMIVARIOGRAMA AJUSTADO: CONSISTE EN EL MEJOR AJUSTE POSIBLE ENTRE EL SEMIVARIOGRAMA MUESTRAL Y EL MODELO
dev.new();plot(Pm.sv,Pm.sph.f,pch="+", pl=TRUE,col="black", main="Pm") #GRAFICO DEL SEMIVARIOGRAMA AJUSTADO CON LOS PUNTOS DEL MUESTRAL

##EL SEMIVARIOGRAMA AJUSTADO CONTIENE LA INFORMACIÓN NECESARIA PARA REALIZAR LA INTERPOLACION DEL KRIGGING SOBRE UNA MALLA VACIA
###GENERAR UNA MALLA VACIA
grd <- expand.grid(x=seq(from=min(d.env.t@bbox['x',]), to=max(d.env.t@bbox['x',]), by=5000),y=seq(from=min(d.env.t@bbox['y',]), to=max(d.env.t@bbox['y',]), by=5000)) 
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) INDICANDOLE QUE COLUMNAS CONTIENEN LAS COORDENADAS x E y
dev.new();spplot(grd) #VISUALIZA EL OBJETO grd. DADO QUE 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

###INTERPOLACION
OKPm <- krige(id="Pm", formula=Pm~1,d.env.t, newdata=grd, model = Pm.sph.f) #KRIGGING ORDINARIO. NOTESE QUE EL SEMIVARIOGRAMA AJUSTADO ENTRA COMO MODELO EN EL KRIGGING. LA INTERPOLACION SE REALIZA SOBRE LA MALLA VACIA
str(OKPm) #ESTRUCTURA DE LA INTERPOLACION
gridded(OKPm)<-TRUE #CUADRICULANDO LA INTERPOLACION
dev.new();plot(OKPm)
str(OKPm) #NUEVA ESTRUCTURA DE LA INTERPOLACION

###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 ESTACION DE BAYAGUANA COINCIDIA CON LA DE RANCHO ARRIBA
plotKML(OKPm["Pm.pred"], colour_scale = SAGA_pal[[9]]) #ABRE GOOGLEEARTH Y MUESTRA LA MALLA INTERPOLADA
OKPm.p <- grid2poly(OKPm["Pm.pred"]) #CONVIERTE LA MALLA A POLIGONOS
OKPm.p@data$Pm.pred <- round(OKPm.p@data$Pm.pred,0) #REDONDEA VALORES DE PRECIPITACION
plotKML(OKPm.p["Pm.pred"], colour_scale = SAGA_pal[[9]]) #ABRE GOOGLEEARTH Y MUESTRA LA MALLA INTERPOLADA CONVERTIDA A POLIGONOS
writeGDAL(OKPm["Pm.pred"],"Pm.pred.tif", drivername="GTiff") #CONVIERTE LA MALLA INTERPOLADA A GTIFF


#ALGUNOS RECURSOS

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

##DATOS DE SALUD
###Ministerio de Salud Publica (MSP); Organizacion Panamericana de la Salud (OPS); Organizacion Mundial de la Salud (2013): "Indicadores basicos de salud 2013". URL: http://www.bvs.org.do/bvs/htdocs//local/File/indicadores--basicos-de-salud-2013.pdf

##DATOS ONE (UTILES 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 ESTE 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 ANALISIS 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