TareaClase11
AndersonCT
22/11/2020
Instalar paquetes
Instalamos los paquetes readxl para el cargue de los datos en formato de excel y ggplot2 para visualizar datos.
## Intalar paquetes
library(readxl)
library(ggplot2)Cargue de los datos
Cargamos los datos del sector de Fresno, el cual tiene informaciĆ³n meteorolĆ³gica espacializada en terreno.
## Importar datos
Fresno_1_Da <- read_excel("D:/Posgrado/1. Primer semestre/1. Clase tratamiento de datos espaciales/Clase 11/Fresno_1_Da.xlsx", sheet = "Sheet1")
Fresno_1_Da## # A tibble: 2,363 x 9
## id_arbol Latitude Longitude Date2 `Psychro Wet Bu~ `Station Pressu~
## <chr> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 1 2.38 -76.6 21/0~ 14.8 805.
## 2 2 2.38 -76.6 21/0~ 11.6 805.
## 3 3 2.38 -76.6 21/0~ 12.9 806.
## 4 4 2.38 -76.6 21/0~ 14.1 806.
## 5 5 2.38 -76.6 21/0~ 14.3 805.
## 6 6 2.38 -76.6 21/0~ 14.2 805.
## 7 7 2.38 -76.6 21/0~ 14.4 805.
## 8 8 2.38 -76.6 21/0~ 12.8 805.
## 9 9 2.38 -76.6 21/0~ 15 805
## 10 10 2.38 -76.6 21/0~ 14 805.
## # ... with 2,353 more rows, and 3 more variables: Relative_Humidity <dbl>,
## # Crosswind <dbl>, Temperature <dbl>names(Fresno_1_Da)## [1] "id_arbol" "Latitude"
## [3] "Longitude" "Date2"
## [5] "Psychro Wet Bulb Temperature" "Station Pressure"
## [7] "Relative_Humidity" "Crosswind"
## [9] "Temperature"Graficos exploratorios del terreno
Planteamos algunos grĆ”ficos para observar los datos en terreno. Primero Analizamos la distribuciĆ³n de los puntos en este mediante las variables de longitud y latitud.
plot(Fresno_1_Da[,c(3,2)])
Ahora, le asociamos a los puntos la variable humedad relativa, la cual fue necesario ajustarle el nombre en la base de datos a Relative_Humidity. El grƔfico se realiza por cada muestreo temporal realizado en el aƱo 2019.
ggplot(Fresno_1_Da,aes(x=Longitude,y=Latitude,color=Relative_Humidity))+geom_point()+theme_bw()+facet_wrap(~Date2)+
scale_color_continuous(type = "viridis")
Filtros de la base de datos
Filtramos el instante de tiempo 1, equivalente al 01/10/2019.
##Filtramos el instante 1 de tiempo
fechas <- unique(Fresno_1_Da$Date2)
pos <- which(Fresno_1_Da$Date2==fechas[4])
datos_m1 <- Fresno_1_Da[pos,]
datos_m1## # A tibble: 394 x 9
## id_arbol Latitude Longitude Date2 `Psychro Wet Bu~ `Station Pressu~
## <chr> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 1 2.38 -76.6 1/10~ 17.4 799.
## 2 2 2.38 -76.6 1/10~ 18.4 800.
## 3 3 2.38 -76.6 1/10~ 17.4 800
## 4 4 2.38 -76.6 1/10~ 16.9 800.
## 5 5 2.38 -76.6 1/10~ 17.3 800.
## 6 6 2.38 -76.6 1/10~ 18.4 799.
## 7 7 2.38 -76.6 1/10~ 16.6 799.
## 8 8 2.38 -76.6 1/10~ 17.1 799.
## 9 9 2.38 -76.6 1/10~ 17 800.
## 10 10 2.38 -76.6 1/10~ 16.4 800.
## # ... with 384 more rows, and 3 more variables: Relative_Humidity <dbl>,
## # Crosswind <dbl>, Temperature <dbl>Uso de librerĆa leaflet
Para visualizar el mapa, con un mapa base de fondo, como el de open street maps, utilizamos la librerĆa leaflet.
library(leaflet) ##Para hacer mapas iteractivos como open street
leaflet() %>% addCircleMarkers(lng = datos_m1$Longitude,lat = datos_m1$Latitude, radius = 0.05,color="black",opacity = 0.9,label=datos_m1$id_arbol) %>% addTiles()GeoestadĆstica, evaluaciĆ³n de los datos
Mediante los siguientes grƔficos evaluamos la correlacion espacial y problemas de estacionaridad.
library(geoR)## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.8-1 (built on 2020-02-08) is now loaded
## --------------------------------------------------------------geo_dato <- as.geodata(datos_m1,coords.col = c(3,2), data.col = 7)
plot(geo_dato)
Calculo del semivariograma
Se realiza el calculo del semivariograma y su regiĆ³n. Observando que los datos tienen correlaciĆ³n alta a cortas distancias y pierde este factor a medida que las distancias aumentan. El Primer valor de donde inicia, podria ser el minimo, segundo valor es el 3 cuartil equivalente a donde se estabiliza, y el tercer valor es el mĆnimo, todo esto lo tomamos del summary.
plot(variograma <- variog(geo_dato,option = "bin"))## variog: computing omnidirectional variogram
variograma <- variog(geo_dato,option = "bin",uvec=seq(0,0.0009703,0.00003)) ## variog: computing omnidirectional variogramplot(variograma)
La regiĆ³n no atrapa casi datos.
variograma.env <- variog.mc.env(geo_dato,obj=variograma)## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelopsplot(variograma,envelope=variograma.env)
Ajuste del Modelo de Semivarianza
Modelos para el ajuste de la semivarainza, tales como el exponencia, gaussiano y el esferico.
ini.vals <- expand.grid(seq(8,12,l=10), seq(0.0006,0.0008,l=10))
model_mco_exp <- variofit(variograma, ini=ini.vals, cov.model="exponential", wei="npair", min="optim")## variofit: covariance model used is exponential
## variofit: weights used: npairs
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
## sigmasq phi tausq kappa
## initial.value "12" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 1860668904.52791## Warning in variofit(variograma, ini = ini.vals, cov.model = "exponential", :
## unreasonable initial value for sigmasq + nugget (too low)model_mco_gaus <- variofit(variograma, ini=ini.vals, cov.model="gaussian", wei="npair", min="optim",nugget = 0)## variofit: covariance model used is gaussian
## variofit: weights used: npairs
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
## sigmasq phi tausq kappa
## initial.value "12" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 1866115182.83403## Warning in variofit(variograma, ini = ini.vals, cov.model = "gaussian", :
## unreasonable initial value for sigmasq + nugget (too low)model_mco_spe <- variofit(variograma, ini=ini.vals, cov.model="spheric",fix.nug=TRUE, wei="npair", min="optim")## variofit: covariance model used is spherical
## variofit: weights used: npairs
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
## sigmasq phi tausq kappa
## initial.value "12" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 1786002872.60648## Warning in variofit(variograma, ini = ini.vals, cov.model = "spheric", fix.nug =
## TRUE, : unreasonable initial value for sigmasq + nugget (too low)plot(variograma,envelope=variograma.env)
lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="purple")
Aparentemente ningun modelo tiene un buen ajuste.
Escoger el modelo y krigin
Auqnue ninguno presenta ajuste, continuaremos el ejercicio con el modelo gaussiano, para la toma de parƔmetros para el Krigin.
model_mco_gaus## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
## tausq sigmasq phi
## 165.1571 105.5938 0.0014
## Practical Range with cor=0.05 for asymptotic range: 0.002501026
##
## variofit: minimised weighted sum of squares = 107701229Lienzo de coordenadas.
min(geo_dato$coords[,1])## [1] -76.61415max(geo_dato$coords[,1])## [1] -76.61254min(geo_dato$coords[,2])## [1] 2.380799max(geo_dato$coords[,2])## [1] 2.382694Krige.
geodatos_grid=expand.grid(lat=seq(-76.61415,-76.61254,l=100),long=seq(2.380799,2.382694,l=100))
plot(geodatos_grid)
points(geo_dato$coords,col="red",pch=16)
geodatos_ko=krige.conv(geo_dato, loc=geodatos_grid,
krige= krige.control(nugget=0,trend.d="cte",
trend.l="cte",cov.pars=c(sigmasq=105.5938, phi=0.0014)))## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhoodPredicciĆ³n del kriggin - grafica de contorno
par(mfrow=c(1,1))
image(geodatos_ko, main="kriging Predict", xlab="Lat", ylab="Long")
contour(geodatos_ko,main="kriging Predict", add=TRUE, drawlabels=TRUE)
DesviaciĆ³n estandar de la predicciĆ³n por kriging
image(geodatos_ko, main="kriging StDv Predicted",val=sqrt(geodatos_ko$krige.var), xlab="East", ylab="North")
contour(geodatos_ko,main="kriging StDv Predict",val=sqrt(geodatos_ko$krige.var), add=TRUE, drawlabels=TRUE)
MAE (Media Absoluta del Error)
valida <- xvalid(geo_dato,model = model_mco_gaus)## xvalid: number of data locations = 394
## xvalid: number of validation locations = 394
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394,
## xvalid: end of cross-validationMAE <- mean(abs(valida$error)) #Mean Absolute Error
MAE## [1] 8.567895Raster de la humedad relativa en la finca
library(raster)## Loading required package: spmat_imagen <- data.frame(geodatos_grid,geodatos_ko$predict)
temp_m1_raster <- rasterFromXYZ(mat_imagen)
plot(temp_m1_raster)
finca <- shapefile("D:/Posgrado/1. Primer semestre/1. Clase tratamiento de datos espaciales/Clase 11/Shape/Fresno_area.shp")
plot(finca,add=TRUE)
Raster recortado a la zona de la finca
temp_m1_raster_finca <- mask(temp_m1_raster,finca)
plot(temp_m1_raster_finca)
plot(finca,add=TRUE)
VisualizaciĆ³n final con open street maps
crs(temp_m1_raster_finca) <- CRS("+init=epsg:4326")
library(rasterVis)## Loading required package: lattice## Loading required package: latticeExtra##
## Attaching package: 'latticeExtra'## The following object is masked from 'package:ggplot2':
##
## layerlevelplot(temp_m1_raster_finca,par.settings=BuRdTheme)
library(leaflet)
leaflet() %>% addTiles() %>% addRasterImage(temp_m1_raster_finca, opacity = 0.9, colors="Spectral")## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition