Created
December 2, 2016 20:33
-
-
Save geofis/5e5a93556e9ffd473f8113f0571a9b69 to your computer and use it in GitHub Desktop.
Análisis geoespacial y geoestadística con R
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment