Método de interpolacion espacial Kriging

## Importar datos

library(readxl)
fresmo_1_Da <- read_excel("Fresno_1_Da.xlsx",sheet ="Sheet1")
plot (fresmo_1_Da[,3:2])


## Mapa de localizacion de ofertas 

require(ggplot2)
## Loading required package: ggplot2

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 por fecha del 21 
 
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 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>
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)

##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")

image(geodatos_ko, main="kriging StDv Predicted",val=sqrt(geodatos_ko$krige.var), xlab="East", ylab="North")

## 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("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)
## Loading required package: 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