Análisis geoespacial y geoestadística con 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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment