Método de interpolación Kriging - finca

library(readxl)
fresmo_1_Da <- read_excel("~/Especializacion_geomatica/Tratamiento_de_datos_espaciales/geoestadistica4/Fresno_1_Da.xlsx")
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in F2365 / R2365C6: got '***'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in G2365 / R2365C7: got '***'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in H2365 / R2365C8: got '***'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in I2365 / R2365C9: got '***'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in J2365 / R2365C10: got '***'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in K2365 / R2365C11: got '***'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in L2365 / R2365C12: got '***'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in M2365 / R2365C13: got '***'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in N2365 / R2365C14: got '***'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in O2365 / R2365C15: got '***'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in P2365 / R2365C16: got '***'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in Q2365 / R2365C17: got '***'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in R2365 / R2365C18: got '***'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in S2365 / R2365C19: got '***'
## Warning in read_fun(path = enc2native(normalizePath(path)), sheet_i = sheet, :
## Expecting numeric in T2365 / R2365C20: got '***'
## mapa de localizacion de ofertas 
require(ggplot2)
## Loading required package: ggplot2
ggplot(fresmo_1_Da,aes(x=Longitude,y=Latitude))+geom_point()+theme_bw()

ggplot(fresmo_1_Da,aes(x=Longitude,y=Latitude))+geom_point()+theme_bw()+facet_wrap(~Date2)

ggplot(fresmo_1_Da,aes(x=Longitude,y=Latitude,color=Temperature))+geom_point()+theme_bw()+
  facet_wrap(~Date2)+ scale_color_continuous(type = "viridis")

## filtramos la base de datos 
fechas=unique(fresmo_1_Da$Date2)
pos=which(fresmo_1_Da$Date2==fechas[4])
datos_m1=fresmo_1_Da[pos,]
datos_m1
## # A tibble: 394 x 20
##    id_arbol Latitude Longitude Date2 `FORMATTED DATE~ `Psychro Wet Bu~
##    <chr>       <dbl>     <dbl> <chr> <chr>                       <dbl>
##  1 1            2.38     -76.6 1/10~ 1/10/2019  3:19~             17.4
##  2 2            2.38     -76.6 1/10~ 1/10/2019  3:19~             18.4
##  3 3            2.38     -76.6 1/10~ 1/10/2019  3:18~             17.4
##  4 4            2.38     -76.6 1/10~ 1/10/2019  3:18~             16.9
##  5 5            2.38     -76.6 1/10~ 1/10/2019  3:17~             17.3
##  6 6            2.38     -76.6 1/10~ 1/10/2019  3:17~             18.4
##  7 7            2.38     -76.6 1/10~ 1/10/2019  3:09~             16.6
##  8 8            2.38     -76.6 1/10~ 1/10/2019  3:09~             17.1
##  9 9            2.38     -76.6 1/10~ 1/10/2019  3:08~             17  
## 10 10           2.38     -76.6 1/10~ 1/10/2019  3:07~             16.4
## # ... with 384 more rows, and 14 more variables: `Station Pressure` <dbl>,
## #   `Relative Humidity` <dbl>, Crosswind <dbl>, Temperature <dbl>, `Barometric
## #   Pressure` <dbl>, Headwind <dbl>, `Direction - True` <dbl>, `Direction -
## #   Mag` <dbl>, `Wind Speed` <dbl>, `Heat Stress Index` <dbl>, Altitude <dbl>,
## #   `Dew Point` <dbl>, `Density Altitude` <dbl>, `Wind Chill` <dbl>
require(leaflet)
## Loading required package: leaflet
leaflet() %>% addCircleMarkers(lng = datos_m1$Longitude,lat = datos_m1$Latitude,
 radius = 0.05,color = "black",opacity = 0.9)%>% addTiles()
require(geoR)
## Loading required package: 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=3:2,data.col=9)
plot(geo_dato)

## calcular el semivariograma 
summary(dist(geo_dato[[1]]))
##      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-03
variograma=variog(geo_dato,option = "bin",uvec = seq(0,0.0009703,0.00003))
## variog: computing omnidirectional variogram
plot(variograma)

## bandas de confianza, si hay o no autocorrelacion 
variograma.emv=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 envelops
plot(variograma,envelope=variograma.emv)

##Ajuste del Modelo de Semivarianza
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 "8"     "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 860010.762926661
## Warning in variofit(variograma, ini = ini.vals, cov.model = "exponential", :
## unreasonable initial value for sigmasq (too high)
## Warning in variofit(variograma, ini = ini.vals, cov.model = "exponential", :
## unreasonable initial value for sigmasq + nugget (too high)
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 "8"     "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 664428.506291517
## Warning in variofit(variograma, ini = ini.vals, cov.model = "gaussian", :
## unreasonable initial value for sigmasq (too high)
## Warning in variofit(variograma, ini = ini.vals, cov.model = "gaussian", :
## unreasonable initial value for sigmasq + nugget (too high)
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 "8"     "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 2348861.41957811
## Warning in variofit(variograma, ini = ini.vals, cov.model = "spheric", fix.nug =
## TRUE, : unreasonable initial value for sigmasq (too high)
## Warning in variofit(variograma, ini = ini.vals, cov.model = "spheric", fix.nug =
## TRUE, : unreasonable initial value for sigmasq + nugget (too high)
lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="purple")

#---- el gaussiono presenta mayor ajuste 

##Predicción Espacial Kriging

min(geo_dato$coords[,1])
## [1] -76.61415
max(geo_dato$coords[,1])
## [1] -76.61254
min(geo_dato$coords[,2])
## [1] 2.380799
max(geo_dato$coords[,2])
## [1] 2.382694
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=12.6072, phi=0.0006)))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
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)

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)

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-validation
MAE=mean(abs(valida$error))
MAE
## [1] 0.1948358
## Cargamos el ratser
require(raster)
## Loading required package: raster
## Loading required package: sp
mat_imagen=data.frame(geodatos_grid,geodatos_ko$predict)
temp_m1_raster=rasterFromXYZ(mat_imagen)
plot(temp_m1_raster)

finca=shapefile("~/Especializacion_geomatica/Tratamiento_de_datos_espaciales/geoestadistica4/Fresno_area.shp")


temp_m1_raster_finca=mask(temp_m1_raster,finca)
plot(temp_m1_raster_finca)
plot(finca,add=TRUE)

crs(temp_m1_raster_finca)=crs(finca)
crs(temp_m1_raster_finca)=CRS("+init=epsg:4326")

require(rasterVis)
## Loading required package: rasterVis
## Loading required package: lattice
## Loading required package: latticeExtra
## 
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
## 
##     layer
levelplot(temp_m1_raster_finca,par.settings=BuRdTheme)

require(leaflet)
leaflet() %>% addTiles() %>% addRasterImage(temp_m1_raster_finca,opacity = 0.7, 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