R/geoR-Map-BASE.R

In carlosperezoft/hipervizr: Dashboard para Hipercartas y Cartas de Control de Red Rio.

# [24.may.22 4pm]
# Uso del paquete geoR para manejo de informacion geografica o mapas.
#
# IMPORTANTE: Tiene como dependencias los paquetes: RandomFields_3.3.14 y geoR_1.8-1
#
library(geoR, quietly=TRUE)    # v1.8-1
library(plotly, quietly=TRUE)  # v4.10.0
library(ggplot2, quietly=TRUE) # v3.3.5
library(scatterplot3d, quietly=TRUE)  #v0.3-41
library(plot3D, quietly=TRUE)  #v1.4
library(sm, quietly=TRUE)  #v2.2-5.7
library(npsp, quietly=TRUE)  #v0.5-3
# data(package = "geoR")  # lista los conjuntos de datos en el paquete geoR

data(wolfcamp) # carga el archivo de datos wolfcamp
head(wolfcamp)
summary(wolfcamp)
#
# genera por defecto gráficos de los valores en las posiciones espaciales
plot(wolfcamp)
# También, en lugar del histograma, nos puede interesar un gráfico de dispersión 3D
plot(wolfcamp, lowess = TRUE, scatter3d = TRUE)
# Si se asume que hay una tendencia puede interesar eliminarla:
plot(wolfcamp, trend=~coords)
#
# genera un gráfico con las posiciones de los datos
points(wolfcamp)
# Se pueden establecer los tamaños de los puntos, simbolos y colores a partir de los valores de los datos.
points(wolfcamp, col = "gray", pt.divide = "equal")
##********************
##*  Modelado de la dependencia
##********************
data(s100) # Cargar datos estacionarios
head(s100)
summary(s100)
#
plot(s100)
#*****
#* Variogramas empíricos
#******
oldpar <- par(mfrow=c(1,2))
plot(variog(s100))
#
vario <- variog(s100, max.dist = 0.6)
plot(vario)
#
vario.b <- variog(s100, max.dist = 0.6) #discretizado
# variog: computing omnidirectional variogram
vario.c <- variog(s100, max.dist=0.6, op="cloud")  #nube
# variog: computing omnidirectional variogram
vario.bc <- variog(s100, max.dist=0.6, bin.cloud=TRUE)  #discretizado+nube
# variog: computing omnidirectional variogram
vario.s <- variog(s100, max.dist=0.6, op="sm", band=0.2)  #suavizado
# variog: computing omnidirectional variogram
# Representación gráfica
oldpar<-par(mfrow=c(2,2)) # Preparar para 4 gráficos por ventana
plot(vario.b, main="Variograma empírico")
plot(vario.c, main="Nube de puntos variograma")
plot(vario.bc, bin.cloud=TRUE, main="Graficos de cajas")
title("Gráficos de cajas") # Corregir fallo del comando anterior
head(vario.s)
plot(vario.s, main="Variograma suavizado")
#
#*************
#*Predicción espacial (kriging)
#*************
vario.ols <- variofit(vario.b, ini = c(1, 0.5), weights = "equal")  #ordinarios
vario.wls <- variofit(vario.b, ini = c(1, 0.5), weights = "cressie")  #ponderados
##
xx <- seq(0, 1, l = 51)
yy <- seq(0, 1, l = 51)
pred.grid <- expand.grid(x = xx, y = yy)
plot(s100$coords, pch = 20)
points(pred.grid, pch = 3, cex = 0.2)
#
ko.wls <- krige.conv(s100, loc = pred.grid, krige = krige.control(obj.m = vario.wls))
#
oldpar <- par(mfrow = c(1, 2))
image(ko.wls) #superficie de predicción
title("Predicciones")
points(s100$coords, pch=20) #añadir posiciones datos
contour(ko.wls,add=T) #añadir gráfico de contorno
#
image(ko.wls, val = ko.wls$krige.var) #superficie de varianzas
title("Superficie de varianzas")
points(s100$coords, pch=20)
contour(ko.wls,val=sqrt(ko.wls$krige.var),add=T)
#
par(oldpar)
contour(ko.wls,filled = TRUE)
#
fcol <- topo.colors(10)[cut(matrix(ko.wls$pred,nrow=51,ncol=51)[-1,-1],10,include.lowest=TRUE)]
persp(ko.wls, theta=-60, phi=40, col=fcol)
head(ko.wls)
xx
predictMtx <- matrix(ko.wls$predict, nrow = length(xx))
# TIP: # El uso de la "~" le indica al plot_ly que use el nombre de la variable
# para los datos y el "Label" del respectivo eje
#
fig <- plot_ly(z=~predictMtx, type = "heatmap")
fig
fig <- plot_ly(z=~predictMtx, type = "contour")
fig
#fig <- plot_ly(x=~xx, y=~yy, z=~predictMtx)
fig <- plot_ly(z=~predictMtx)
fig <- fig %>% add_surface()
fig
#
# Surface Plot With Contours, FULL OK!!
fig <- plot_ly(z=~predictMtx)
fig <- fig %>% add_surface(
   contours = list(
      z = list(
         show=TRUE,
         usecolormap=TRUE,
         highlightcolor="#ff0000",
         project=list(z=TRUE)
      )
   )
)
fig <- fig %>% layout(
   scene = list(
      camera=list(
         eye = list(x=1.87, y=0.88, z=-0.64)
      )
   )
)
fig
#
persp3D(xx, yy, matrix(ko.wls$predict, nrow = length(xx)), theta=-60, phi=40)
#
spersp(xx, yy, ko.wls$predict, theta=-60, phi=40)
contour(spersp(xx, yy, ko.wls$predict, theta=-60, phi=40))
#