Entregable_3
David Pinedo
2023-10-06
# Carga de las librerÃas
library(raster)## Loading required package: sp## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
## (status 2 uses the sf package in place of rgdal)library(rgdal)## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/david/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/david/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:2.0-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.library(sp)
library(here)## here() starts at C:/Users/david/OneDrive - PUJ Cali/Análisis Información Geográficalibrary(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(readxl)
require(leaflet)## Loading required package: leafletrequire(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.9-2 (built on 2022-08-09) is now loaded
## --------------------------------------------------------------require(RColorBrewer)## Loading required package: RColorBrewerrequire(rasterVis)## Loading required package: rasterVis## Loading required package: lattice# Ruta de archivos
ruta <- here()# Importación de datos
aguacates <- read_xlsx(paste(ruta,"Datos_Completos_Aguacate.xlsx", sep = "/"))
head(aguacates)## # A tibble: 6 × 21
## id_arbol Latitude Longitude FORMATTED_DATE_TIME Psychro_Wet_Bulb_Tempe…¹
## <chr> <dbl> <dbl> <chr> <dbl>
## 1 1 2.38 -76.6 21/08/2019Â 9:22:57 a, m, 14.8
## 2 2 2.38 -76.6 21/08/2019Â 9:27:13 a, m, 11.6
## 3 3 2.38 -76.6 21/08/2019Â 9:36:36 a, m, 12.9
## 4 4 2.38 -76.6 21/08/2019Â 9:38:02 a, m, 14.1
## 5 5 2.38 -76.6 21/08/2019Â 9:39:38 a, m, 14.3
## 6 6 2.38 -76.6 21/08/2019Â 9:42:02 a, m, 14.2
## # ℹ abbreviated name: ¹​Psychro_Wet_Bulb_Temperature
## # ℹ 16 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>,
## # Estado_Fenologico_Predominante <dbl>, Frutos_Afectados <dbl>aguacates_limpia <- aguacates %>%
filter(grepl("01/10/2020",
FORMATTED_DATE_TIME))leaflet() %>%
addTiles() %>%
addCircleMarkers(lng = aguacates_limpia$Longitude,
lat = aguacates_limpia$Latitude,
radius = 0.2,
color = "blue")Observamos que la finca se encuentra en zona rural de Timbio, Cauca.
GeostadÃstica
geo_temp <- as.geodata(aguacates_limpia,coords.col = 3:2,data.col = 9)
plot(geo_temp)
summary(dist(aguacates_limpia[,3:2]))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.712e-05 4.051e-04 6.408e-04 6.827e-04 9.178e-04 1.959e-03variograma=variog(geo_temp,option = "bin",uvec=seq(0,0.001,0.00002))## variog: computing omnidirectional variogramdatos_env <- 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 envelopsplot(variograma)
lines(datos_env)
Se observa en el semivariograma una autocorrelación espacial no muy fuerte indicando que las mediciones cercanas pero las distantes tampoco distan demasiado de las lejanas.
Ajuste del modelo teórico
ini_vals = expand.grid(seq(1.2,1.5,l=10), seq(0.0001,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 "1.5" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 296807.231725087model_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 "1.5" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 288205.86784414model_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 "1.5" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 283743.808438713model_mco_spe## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## fixed value for tausq = 0
## parameter estimates:
## sigmasq phi
## 3.0186 0.0000
## Practical Range with cor=0.05 for asymptotic range: 0
##
## variofit: minimised weighted sum of squares = 17731.76plot(variograma)
lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="purple")
Se observa en los gráficos y en el resultado de SCE que el modelo esférico se ajusta de mejor manera a los datos.
Ahora se realizará la interpolación
max(aguacates_limpia[,3])## [1] -76.71022list(min(aguacates_limpia[,3]),
max(aguacates_limpia[,3]),
min(aguacates_limpia[,2]),
max(aguacates_limpia[,2]))## [[1]]
## [1] -76.7118
##
## [[2]]
## [1] -76.71022
##
## [[3]]
## [1] 2.392101
##
## [[4]]
## [1] 2.393634geo_aguacates_grid <- expand.grid( lon=seq(-76.710215,-76.711799,l=100),
lat=seq(2.392101 ,2.393634 ,l=100))
plot(geo_aguacates_grid)
points(aguacates_limpia[,3:2],col="blue")
#Predicción espacial krigging
geo_aguacates_krigging <- krige.conv(geo_temp, locations = geo_aguacates_grid,
krige= krige.control(nugget=0,trend.d="cte",
trend.l="cte",
cov.pars=c(sigmasq=3.0186,
phi=0.0001 )))## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhoodpar(mfrow=c(1,2))
image(geo_aguacates_krigging, main="kriging Predict", xlab="East", ylab="North")
contour(geo_aguacates_krigging,main="kriging Predict", drawlabels=TRUE)
Transformación de imagen en raster
prediccion <- cbind(geo_aguacates_grid,geo_aguacates_krigging$predict)
temp_prediccion <- rasterFromXYZ(cbind(geo_aguacates_grid,geo_aguacates_krigging$predict))
levelplot(temp_prediccion,par.settings =BuRdTheme)
temp_error <- rasterFromXYZ(cbind(geo_aguacates_grid,sqrt(geo_aguacates_krigging$krige.var)))
levelplot(temp_error,par.settings =BuRdTheme)