Geo Datos Clima
David Arango
4/19/2021
Importar Datos y Librerias
library(readxl)
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
## --------------------------------------------------------------library(raster)## Loading required package: splibrary(ggplot2)
library(leaflet)
geo_datos =read_excel("~/Desktop/geo_datos.xlsx")
geo_datos## # A tibble: 394 x 7
## id_arbol Latitude Longitude `Relative Humidity` Temperature `Wind Speed`
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 2.38 -76.6 61.2 22.9 0
## 2 2 2.38 -76.6 62.3 23.8 0
## 3 3 2.38 -76.6 62.6 22.6 0.4
## 4 4 2.38 -76.6 61.3 22.3 0.6
## 5 5 2.38 -76.6 61.9 22.7 0.8
## 6 6 2.38 -76.6 65.7 23.1 0.4
## 7 7 2.38 -76.6 67 21 0.5
## 8 8 2.38 -76.6 65.8 21.6 0.4
## 9 9 2.38 -76.6 66.6 21.4 0
## 10 10 2.38 -76.6 66.3 20.8 0.7
## # … with 384 more rows, and 1 more variable: Altitude <dbl>Se observa que los datos contienen información sobre las localizaciones de los sitios donde se tomaron datos de clima como temperatura, humedad relativa entre otras.
leaflet() %>% addTiles() %>% addCircleMarkers(lng = geo_datos$Longitude,lat = geo_datos$Latitude,radius = 0.2,color = "black")# leaflet() %>% addTiles() %>% addCircleMarkers(lng = geo_datos$Longitude,lat = geo_datos$Latitude,radius = 0.2,color = "black",clusterOptions = markerClusterOptions())Se observa que la finca se encuentra en la zona rural del municipio de Popayan y su distribucción espacial es bastante regular.
A. Geoestadística
Como primer paso realizamos la creación de la variable regionalizada:
geo_temp=as.geodata(geo_datos,coords.col = 3:2,data.col = 5)
plot(geo_temp)
summary(dist(geo_datos[,3:2]))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.302e-05 4.017e-04 6.544e-04 7.038e-04 9.703e-04 1.912e-03variograma=variog(geo_temp,option = "bin",uvec=seq(0,0.0009703,9.703e-05))## variog: computing omnidirectional variogramplot(variograma)
variograma_mc=variog.mc.env(geo_temp,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 envelopslines(variograma_mc)
Se observa en el semivariograma una autocorrelación espacial muy fuerte indicando que las mediciones cercanas, presentan temperaturas similares en contraste con las distantes.
Ahora vamos a realizar el ajuste del modelo teorico:
##Ajuste del Modelo de Semivarianza
ini.vals = expand.grid(seq(8,12,l=10), seq(6e-04,8e-04,l=10))
plot(variograma)
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: 174964.991544732model_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: 130506.318494236model_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 "10.22" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 107137.155200497lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="purple")
Podemos observar tanto graficamente como por el valor de SCE que el modelo gaussiano es el que mejor se ajusta a los puntos del semivariograma muestral.
Ahora vamos a realizar la interpolación (predicción espacial)
max(geo_datos[,3])## [1] -76.61254##Predicción Espacial Kriging
geodatos_grid=expand.grid( lon=seq(-76.61415,-76.61254,l=100),lat=seq(2.380799 ,2.382694 ,l=100))
plot(geodatos_grid)
points(geo_datos[,3:2],col="red")
geodatos_ko=krige.conv(geo_temp, loc=geodatos_grid,
krige= krige.control(nugget=0,trend.d="cte",
trend.l="cte",cov.pars=c(sigmasq=12.7080, phi=0.0007 )))## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhoodpar(mfrow=c(1,2))
image(geodatos_ko, main="kriging Predict", xlab="East", ylab="North")
contour(geodatos_ko,main="kriging Predict", add=TRUE, drawlabels=TRUE)
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)
Ahora vamos a tranformar la imagen en raster:
require(rasterVis)## Loading required package: rasterVis## Loading required package: terra## terra version 1.1.4## Loading required package: lattice## Loading required package: latticeExtra##
## Attaching package: 'latticeExtra'## The following object is masked from 'package:ggplot2':
##
## layerpred=cbind(geodatos_grid,geodatos_ko$predict)
temp_predict=rasterFromXYZ(cbind(geodatos_grid,geodatos_ko$predict))
plot(temp_predict)
levelplot(temp_predict,par.settings =BuRdTheme)
temp_error=rasterFromXYZ(cbind(geodatos_grid,sqrt(geodatos_ko$krige.var)))
levelplot(temp_error,par.settings =BuRdTheme)
por ultimo vamos a cortar los raster con el shape de la finca:
finca=shapefile("~/Desktop/Fresno_area.shp")
plot(temp_predict)
plot(finca,add=T)
temp_predict_finca=mask(temp_predict,finca)
plot(temp_predict_finca)
plot(finca,add=T)
crs(temp_predict_finca)=crs("+proj=longlat +datum=WGS84")
leaflet() %>% addTiles() %>% addRasterImage(temp_predict_finca,opacity = 0.6)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 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 +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 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 +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definitionpor ultimo realizamos la validación cruzada:
valida=xvalid(geodata = geo_temp,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))
MAE## [1] 1.105438El modelo presenta un error de predicción de acuerdo con la validación cruzada de 1.1 grados.