#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