1. Realice un mapa exploratorio con leaflet que permita identificar la ubicación de los arboles con un mapa base de fondo que permita ver las características de la zona.

A continuación, se puede visualizar en leaflet la ubicación de los arboles con el mapa base de Esri World Imagery, el cual permite observar la imagen satelital de la zona de estudio y su entorno.

require(ggplot2)
require(table1)
require(geoR)
require(raster)
require(leaflet)
require(readxl)

datos_frutos <- read_excel("D:/ESP-GEOMATICA/2. TRATAMIENTO DE DATOS ESPACIALES/INFORME FINAL/datos_frutos.xlsx")

datos=datos_frutos

leaflet()%>% addProviderTiles('Esri.WorldImagery')%>%
  addCircleMarkers(lng = datos$long, lat= datos$lat, radius = 0.1, color="#FF0000")
  1. Realice un mapa exploratorio que permita identificar la distribución espacial de los frutos afectados en la finca utilizando ggplot2 o leaflet y una paleta de colores apropiada (ejemplo: azulrojo).

Discutir si efectivamente se encuentran zonas en la finca con mayor o menor cantidad de frutos afectados (cluster).

A continuación se puede visualizar la distribución espacial de los frutos afectados con la libreria leaflet, en este se puede observar que dentro de la finca existen cluster donde el de menos frutos afectados es de 20 y el mayor de 58.

datos_frutos <- read_excel("D:/ESP-GEOMATICA/2. TRATAMIENTO DE DATOS ESPACIALES/INFORME FINAL/datos_frutos.xlsx")

datos=datos_frutos



leaflet()%>% addProviderTiles('Esri.WorldImagery')%>%addCircleMarkers(
lng = datos$long,
lat= datos$lat, 
radius = 0.1,
clusterOptions =datos$frutos_afectados)
  1. Convertir la variable frutos afectados en binaria (1: si tiene frutos afectados, 0: si no tiene frutos afectados) y con base en esta variable calcular con table 1 una tabla que permita comparar los indicadores descriptivos de las variables climáticas comparando si tiene o no frutos afectados table1(~clima|afectados,data). Interpretar los resultados teniendo en cuenta las variables que presentan mayor diferencia en los arboles con y sin frutos afectados.

A continuación se reemplazan los datos donde hay más de 1 fruto afectado sea reemplazado por 1, y con esta se puede realizar table 1 que permite observar las estadisticas teniendo en cuenta las variables climaticas.

En esta tabla en las columnas se encuentra el valor 0 que significa que no hay frutos afectados y 1 donde hay más de un fruto afectado, con base a esta se infiere que la temperatura aumenta, la humedad relativa baja y la velocidad del viento sube teniendo en cuenta el promedio general en los puntos donde se encuentra frutos afectados, sin embargo las diferencias no ascienden en el promedio al 0.4.

datos_frutos <- read_excel("D:/ESP-GEOMATICA/2. TRATAMIENTO DE DATOS ESPACIALES/INFORME FINAL/datos_frutos.xlsx")

datos=datos_frutos

datos$frutos_afectados[datos$frutos_afectados>=1]<-1
datos
## # A tibble: 394 × 8
##    arbol   lat  long humedad_relativa temperatura velocidad_vie…¹ altura fruto…²
##    <chr> <dbl> <dbl>            <dbl>       <dbl>           <dbl>  <dbl>   <dbl>
##  1 1      2.38 -76.6             33.6        25.7             0.8   1896       0
##  2 2      2.38 -76.6             36.8        20.8             5.1   1895       1
##  3 3      2.38 -76.6             31.5        23.7             0.8   1889       0
##  4 4      2.38 -76.6             33.2        25               0.9   1890       0
##  5 5      2.38 -76.6             34.3        25               0.6   1894       1
##  6 6      2.38 -76.6             33.8        25               0.5   1900       0
##  7 7      2.38 -76.6             34.9        24.9             0.6   1900       0
##  8 8      2.38 -76.6             34.2        22.9             3.7   1899       0
##  9 9      2.38 -76.6             33.6        26.2             0.5   1897       0
## 10 10     2.38 -76.6             34          24.6             1     1892       1
## # … with 384 more rows, and abbreviated variable names ¹​velocidad_viento,
## #   ²​frutos_afectados
require(table1)
table1(~temperatura+humedad_relativa+velocidad_viento|frutos_afectados,datos)
0
(N=295)		 1
(N=99)		 Overall
(N=394)
temperatura
Mean (SD) 26.5 (2.38) 26.9 (2.63) 26.6 (2.45)
Median [Min, Max] 26.4 [21.6, 33.9] 27.0 [20.8, 32.0] 26.6 [20.8, 33.9]
humedad_relativa
Mean (SD) 36.5 (6.94) 36.2 (6.77) 36.4 (6.89)
Median [Min, Max] 35.2 [25.1, 76.1] 34.4 [26.7, 67.2] 35.0 [25.1, 76.1]
velocidad_viento
Mean (SD) 0.387 (0.411) 0.389 (0.607) 0.387 (0.467)
Median [Min, Max] 0.400 [0, 3.70] 0.400 [0, 5.10] 0.400 [0, 5.10]
  1. Ajuste un modelo de regresión poisson para predecir la variable total de frutos afectados y valide cuales son las variables climáticas que presentan mayor asociación con la presencia de frutos afectados. mod=glm(frutos_afectados~humedad_relativa+temperatura+velocidad_viento+altura,fa mily = “poisson”,data = datos_frutos)

En el siguiente código se ajusto un modelo de regresión poisson para predecir la variable total de frutos afectados, en este se incluyeron las variables climaticas y de terreno, teniendo en cuenta la significancia que arrojó el resumen del modelo la de mayor es la temperatura, seguida por la altura, por lo contrario la humedad relativa y velocidad del viento no aportaron significancia al modelo.

datos_frutos <- read_excel("D:/ESP-GEOMATICA/2. TRATAMIENTO DE DATOS ESPACIALES/INFORME FINAL/datos_frutos.xlsx")

datos=datos_frutos

mod=glm(frutos_afectados~humedad_relativa+temperatura+velocidad_viento+altura,family = "poisson",data = datos)

summary(mod)
## 
## Call:
## glm(formula = frutos_afectados ~ humedad_relativa + temperatura + 
##     velocidad_viento + altura, family = "poisson", data = datos)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9786  -1.5001  -1.3769  -0.6701  11.1520  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)   
## (Intercept)      35.593554  17.258673   2.062  0.03917 * 
## humedad_relativa -0.009461   0.009357  -1.011  0.31201   
## temperatura       0.076346   0.025474   2.997  0.00273 **
## velocidad_viento  0.121149   0.108412   1.117  0.26379   
## altura           -0.019691   0.009206  -2.139  0.03244 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1594.7  on 393  degrees of freedom
## Residual deviance: 1576.1  on 389  degrees of freedom
## AIC: 1872.6
## 
## Number of Fisher Scoring iterations: 7
mod
## 
## Call:  glm(formula = frutos_afectados ~ humedad_relativa + temperatura + 
##     velocidad_viento + altura, family = "poisson", data = datos)
## 
## Coefficients:
##      (Intercept)  humedad_relativa       temperatura  velocidad_viento  
##        35.593554         -0.009461          0.076346          0.121149  
##           altura  
##        -0.019691  
## 
## Degrees of Freedom: 393 Total (i.e. Null);  389 Residual
## Null Deviance:       1595 
## Residual Deviance: 1576  AIC: 1873
  1. Utilice geoestadística para realizar una interpolación de las variables climáticas y presente todos los pasos para el ajuste del modelo.

Temperatura:

A continuación se utiliza geoestadistica para interpolar la temperatura con el modelo exponencial, que fue el que mejor rendimiento tuvo para describir el fenomeno por tener menor SCE, además se caluclo el poder predictivo del modelo con la metrica MAE la cual dio como resultado 0.92.

require(ggplot2)
require(table1)
require(geoR)
require(raster)
require(leaflet)
require(readxl)
require (sf)
require (dplyr)

datos_frutos <- read_excel("D:/ESP-GEOMATICA/2. TRATAMIENTO DE DATOS ESPACIALES/INFORME FINAL/datos_frutos.xlsx")

datos=datos_frutos


##Interpolación con Geoestadistica para temperatura

geo_temp=as.geodata(datos,coords.col = 3:2,data.col=5)
plot(geo_temp)

summary(dist(geo_temp$coords))
##      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
##Desde 0 hasta el tercer cuartir se puede poner cada cuanto (un poquito más del minimo)

h=seq(0, 9.703e-04, l=10)

semi=variog(geo_temp, uvec=h)
## variog: computing omnidirectional variogram
plot(semi)

## Ajuste modelo teorico

vario_gauss=variofit(semi, cov.model = "gaussian")
## 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 "4.88"  "0"   "0.65" "0.5"
## status        "est"   "est" "est"  "fix"
## loss value: 15703.5520815038
lines(vario_gauss, col="red")

vario_sphe=variofit(semi, cov.model = "spherical")
## 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 "4.88"  "0"   "0.65" "0.5"
## status        "est"   "est" "est"  "fix"
## loss value: 10635.4521420932
lines(vario_sphe, col="blue")

vario_exp=variofit(semi, cov.model = "exponential")
## 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 "4.88"  "0"   "1.63" "0.5"
## status        "est"   "est" "est"  "fix"
## loss value: 7243.9067366975
lines(vario_exp, col="purple")

vario_gauss
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##   tausq sigmasq     phi 
##  0.7847  4.4775  0.0003 
## Practical Range with cor=0.05 for asymptotic range: 0.0004086373
## 
## variofit: minimised weighted sum of squares = 12835.56
vario_sphe
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## parameter estimates:
##   tausq sigmasq     phi 
##  0.6807  4.6492  0.0006 
## Practical Range with cor=0.05 for asymptotic range: 0.0005918533
## 
## variofit: minimised weighted sum of squares = 9464.237
vario_exp ##mejor por tener menor SCE
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
##  1.6260  4.8781  0.0005 
## Practical Range with cor=0.05 for asymptotic range: 0.001510302
## 
## variofit: minimised weighted sum of squares = 7012.658
##Prediccion espacial

##Delimitar el espacio

#plot(datos[,3:2])

min(datos$lat)
## [1] 2.380799
max(datos$lat)
## [1] 2.382694
min(datos$long)
## [1] -76.61415
max(datos$long)
## [1] -76.61254
x=seq(2.380799,2.382694,l=100)
y=seq(-76.61415,-76.61254,l=100)

loc0=expand.grid(y,x)

#plot(loc0)
#points(datos[,3:2],col="red")


temp_pred=krige.conv(geodata = geo_temp,locations = loc0,
                     krige = krige.control(obj.model=vario_exp))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
res=data.frame(loc0,temp_pred$predict)


temp_map=rasterFromXYZ(res)
crs(temp_map)=CRS("+init=epsg:4326")

#temp_map

#plot(temp_map)
#points(datos[,3:2])


#leaflet()%>% addProviderTiles('Esri.WorldImagery')%>%addRasterImage(temp_map, opacity=0.5)


##poder predictivo del modelo

valida=xvalid(geo_temp, model=vario_exp)
## 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.9290965
##cortar la zona de estudio

require(ggplot2)

datos2=data.frame(datos[,3:2])

occr <- st_as_sf(x = datos2, coords = c('long', 'lat'), crs = st_crs(4326))
#plot(st_geometry(occr), pch = 16, col = 'red')

occr <- mutate(occr, GID = 1)


hull <- occr %>% 
  dplyr::group_by(GID) %>% 
  dplyr::summarise() %>%
  st_cast('POLYGON') %>% 
  st_convex_hull()
  
#plot(st_geometry(hull), add = T, col = 'blue', opacity=0.5)
#plot(st_geometry(occr), add = TRUE, col = 'red', pch = 16)

mask_temp_map <- mask(temp_map, mask = hull)
leaflet()%>% addProviderTiles('Esri.WorldImagery')%>%addRasterImage(mask_temp_map, opacity=0.5)

Humedad relativa:

A continuación se utiliza geoestadistica para interpolar la temperatura con el modelo exponencial, que fue el que mejor rendimiento tuvo para describir el fenomeno por tener menor SCE, además se calculo el poder predictivo del modelo con la metrica MAE la cual dio como resultado 3.42.

require(ggplot2)
require(table1)
require(geoR)
require(raster)
require(leaflet)
require(readxl)
require (sf)
require (dplyr)

datos_frutos <- read_excel("D:/ESP-GEOMATICA/2. TRATAMIENTO DE DATOS ESPACIALES/INFORME FINAL/datos_frutos.xlsx")

datos=datos_frutos
#head(datos)

##Interpolación con Geoestadistica para humedad relativa

geo_hum=as.geodata(datos,coords.col = 3:2,data.col=4)
plot(geo_hum)

summary(dist(geo_hum$coords))
##      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
##Desde 0 hasta el tercer cuartir se puede poner cada cuanto (un poquito más del minimo)

h=seq(0, 9.703e-04, l=10)

semi=variog(geo_hum, uvec=h)
## variog: computing omnidirectional variogram
plot(semi)

## Ajuste modelo teorico

vario_gauss=variofit(semi, cov.model = "gaussian")
## 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 "27.18" "0"   "27.18" "0.5"
## status        "est"   "est" "est"   "fix"
## loss value: 609454.449149483
lines(vario_gauss, col="red")

vario_sphe=variofit(semi, cov.model = "spherical")
## 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 "40.77" "0"   "5.44" "0.5"
## status        "est"   "est" "est"  "fix"
## loss value: 1071893.23042456
lines(vario_sphe, col="blue")

vario_exp=variofit(semi, cov.model = "exponential")
## 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 "27.18" "0"   "27.18" "0.5"
## status        "est"   "est" "est"   "fix"
## loss value: 144683.964293356
lines(vario_exp, col="purple")

vario_gauss
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##   tausq sigmasq     phi 
## 24.0933 24.0933  0.0003 
## Practical Range with cor=0.05 for asymptotic range: 0.0004699665
## 
## variofit: minimised weighted sum of squares = 463401.6
vario_sphe
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## parameter estimates:
##   tausq sigmasq     phi 
##  6.1989 40.3345  0.0003 
## Practical Range with cor=0.05 for asymptotic range: 0.0003063143
## 
## variofit: minimised weighted sum of squares = 1063368
vario_exp ##mejor por tener menor SCE
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
## 27.1777 27.1777  0.0005 
## Practical Range with cor=0.05 for asymptotic range: 0.001422768
## 
## variofit: minimised weighted sum of squares = 144684
##Prediccion espacial

##Delimitar el espacio

#plot(datos[,3:2])

min(datos$lat)
## [1] 2.380799
max(datos$lat)
## [1] 2.382694
min(datos$long)
## [1] -76.61415
max(datos$long)
## [1] -76.61254
x=seq(2.380799,2.382694,l=100)
y=seq(-76.61415,-76.61254,l=100)

loc0=expand.grid(y,x)

#plot(loc0)
#points(datos[,3:2],col="red")


hum_pred=krige.conv(geodata = geo_hum,locations = loc0,
                     krige = krige.control(obj.model=vario_exp))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
res=data.frame(loc0,hum_pred$predict)


hum_map=rasterFromXYZ(res)
crs(hum_map)=CRS("+init=epsg:4326")
#hum_map

#plot(hum_map)
#points(datos[,3:2])


#leaflet()%>% addProviderTiles('Esri.WorldImagery')%>%addRasterImage(hum_map, opacity=0.5)


##poder predictivo del modelo

valida=xvalid(geo_hum, model=vario_exp)
## 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] 3.429647
##cortar la zona de estudio

require(ggplot2)

datos2=data.frame(datos[,3:2])

occr <- st_as_sf(x = datos2, coords = c('long', 'lat'), crs = st_crs(4326))
#plot(st_geometry(occr), pch = 16, col = 'red')

occr <- mutate(occr, GID = 1)


hull <- occr %>% 
  dplyr::group_by(GID) %>% 
  dplyr::summarise() %>%
  st_cast('POLYGON') %>% 
  st_convex_hull()
  
#plot(st_geometry(hull), add = T, col = 'blue', opacity=0.5)
#plot(st_geometry(occr), add = TRUE, col = 'red', pch = 16)

mask_hum_map <- mask(hum_map, mask = hull)
leaflet()%>% addProviderTiles('Esri.WorldImagery')%>%addRasterImage(mask_hum_map, opacity=0.5)

Velocidad del viento:

A continuación se utiliza geoestadistica para interpolar la temperatura con el modelo exponencial, que fue el que mejor rendimiento tuvo para describir el fenomeno por tener menor SCE, además se calculo el poder predictivo del modelo con la metrica MAE la cual dio como resultado 0.31.

require(ggplot2)
require(table1)
require(geoR)
require(raster)
require(leaflet)
require(readxl)
require (sf)
require (dplyr)

datos_frutos <- read_excel("D:/ESP-GEOMATICA/2. TRATAMIENTO DE DATOS ESPACIALES/INFORME FINAL/datos_frutos.xlsx")

datos=datos_frutos
#head(datos)

##Interpolación con Geoestadistica para humedad relativa

geo_vien=as.geodata(datos,coords.col = 3:2,data.col=6)
plot(geo_vien)

summary(dist(geo_vien$coords))
##      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
##Desde 0 hasta el tercer cuartir se puede poner cada cuanto (un poquito más del minimo)

h=seq(0, 9.703e-04, l=10)

semi=variog(geo_vien, uvec=h)
## variog: computing omnidirectional variogram
plot(semi)

## Ajuste modelo teorico

vario_gauss=variofit(semi, cov.model = "gaussian")
## 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 "0.11"  "0"   "0.11" "0.5"
## status        "est"   "est" "est"  "fix"
## loss value: 16.8032903945883
lines(vario_gauss, col="red")

vario_sphe=variofit(semi, cov.model = "spherical")
## 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 "0.11"  "0"   "0.11" "0.5"
## status        "est"   "est" "est"  "fix"
## loss value: 15.5805927193741
lines(vario_sphe, col="blue")

vario_exp=variofit(semi, cov.model = "exponential")
## 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 "0.17"  "0"   "0.11" "0.5"
## status        "est"   "est" "est"  "fix"
## loss value: 7.03955430369931
lines(vario_exp, col="purple")

vario_gauss
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##   tausq sigmasq     phi 
##  0.1243  0.0888  0.0004 
## Practical Range with cor=0.05 for asymptotic range: 0.0006491261
## 
## variofit: minimised weighted sum of squares = 10.9637
vario_sphe
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## parameter estimates:
##   tausq sigmasq     phi 
##  0.1121  0.1121  0.0010 
## Practical Range with cor=0.05 for asymptotic range: 0.00100461
## 
## variofit: minimised weighted sum of squares = 7.679
vario_exp ##mejor por tener menor SCE
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
##  0.1121  0.1681  0.0008 
## Practical Range with cor=0.05 for asymptotic range: 0.002448113
## 
## variofit: minimised weighted sum of squares = 6.3613
##Prediccion espacial

##Delimitar el espacio

#plot(datos[,3:2])

min(datos$lat)
## [1] 2.380799
max(datos$lat)
## [1] 2.382694
min(datos$long)
## [1] -76.61415
max(datos$long)
## [1] -76.61254
x=seq(2.380799,2.382694,l=100)
y=seq(-76.61415,-76.61254,l=100)

loc0=expand.grid(y,x)

##plot(loc0)
#points(datos[,3:2],col="red")


vien_pred=krige.conv(geodata = geo_vien,locations = loc0,
                     krige = krige.control(obj.model=vario_exp))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
res=data.frame(loc0,vien_pred$predict)


vien_map=rasterFromXYZ(res)
crs(vien_map)=CRS("+init=epsg:4326")


##plot(vien_map)
##points(datos[,3:2])


#leaflet()%>% addProviderTiles('Esri.WorldImagery')%>%addRasterImage(vien_map, opacity=0.5)


##poder predictivo del modelo

valida=xvalid(geo_vien, model=vario_exp)
## 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.3105212
##cortar la zona de estudio

require(ggplot2)

datos2=data.frame(datos[,3:2])

occr <- st_as_sf(x = datos2, coords = c('long', 'lat'), crs = st_crs(4326))
##plot(st_geometry(occr), pch = 16, col = 'red')

occr <- mutate(occr, GID = 1)


hull <- occr %>% 
  dplyr::group_by(GID) %>% 
  dplyr::summarise() %>%
  st_cast('POLYGON') %>% 
  st_convex_hull()
  
##plot(st_geometry(hull), add = T, col = 'blue', opacity=0.5)
##plot(st_geometry(occr), add = TRUE, col = 'red', pch = 16)
mask_vien_map <- mask(vien_map, mask = hull)
leaflet()%>% addProviderTiles('Esri.WorldImagery')%>%addRasterImage(mask_vien_map, opacity=0.5)

TOPOGRAFIA:

A continuación se utiliza geoestadistica para interpolar la temperatura con el modelo gaussiano, que fue el que mejor rendimiento tuvo para describir el fenomeno por tener menor SCE, además se calculo el poder predictivo del modelo con la metrica MAE la cual dio como resultado 1.455503.

require(ggplot2)
require(table1)
require(geoR)
require(raster)
require(leaflet)
require(readxl)
require (sf)
require (dplyr)

datos_frutos <- read_excel("D:/ESP-GEOMATICA/2. TRATAMIENTO DE DATOS ESPACIALES/INFORME FINAL/datos_frutos.xlsx")

datos=datos_frutos
#head(datos)

##Interpolación con Geoestadistica para humedad relativa

geo_top=as.geodata(datos,coords.col = 3:2,data.col=7)
plot(geo_top)

summary(dist(geo_top$coords))
##      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
##Desde 0 hasta el tercer cuartir se puede poner cada cuanto (un poquito más del minimo)

h=seq(0, 9.703e-04, l=10)

semi=variog(geo_top, uvec=h)
## variog: computing omnidirectional variogram
plot(semi)

## Ajuste modelo teorico

vario_gauss=variofit(semi, cov.model = "gaussian")
## 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 "49.14" "0"   "4.91" "0.5"
## status        "est"   "est" "est"  "fix"
## loss value: 306976.047044923
lines(vario_gauss, col="red")

vario_sphe=variofit(semi, cov.model = "spherical")
## 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 "36.85" "0"   "4.91" "0.5"
## status        "est"   "est" "est"  "fix"
## loss value: 2954749.14806853
lines(vario_sphe, col="blue")

vario_exp=variofit(semi, cov.model = "exponential")
## 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 "49.14" "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 2232861.08256103
lines(vario_exp, col="purple")

vario_gauss ##mejor por tener menor SCE
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##   tausq sigmasq     phi 
##  4.9140 49.1396  0.0006 
## Practical Range with cor=0.05 for asymptotic range: 0.001101614
## 
## variofit: minimised weighted sum of squares = 306976
vario_sphe
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## parameter estimates:
##   tausq sigmasq     phi 
##  4.9140 36.8547  0.0010 
## Practical Range with cor=0.05 for asymptotic range: 0.0009554351
## 
## variofit: minimised weighted sum of squares = 1989554
vario_exp 
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
##  0.0000 49.1396  0.0005 
## Practical Range with cor=0.05 for asymptotic range: 0.001451016
## 
## variofit: minimised weighted sum of squares = 2226952
##Prediccion espacial

##Delimitar el espacio

#plot(datos[,3:2])

min(datos$lat)
## [1] 2.380799
max(datos$lat)
## [1] 2.382694
min(datos$long)
## [1] -76.61415
max(datos$long)
## [1] -76.61254
x=seq(2.380799,2.382694,l=100)
y=seq(-76.61415,-76.61254,l=100)

loc0=expand.grid(y,x)

#plot(loc0)
#points(datos[,3:2],col="red")


top_pred=krige.conv(geodata = geo_top,locations = loc0,
                     krige = krige.control(obj.model=vario_gauss))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
res=data.frame(loc0,top_pred$predict)


top_map=rasterFromXYZ(res)
crs(top_map)=CRS("+init=epsg:4326")


##plot(top_map)
##points(datos[,3:2])


##leaflet()%>% addProviderTiles('Esri.WorldImagery')%>%addRasterImage(top_map, opacity=0.5)


##poder predictivo del modelo

valida=xvalid(geo_top, model=vario_gauss)
## 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] 1.455503
##cortar la zona de estudio

require(ggplot2)

datos2=data.frame(datos[,3:2])

occr <- st_as_sf(x = datos2, coords = c('long', 'lat'), crs = st_crs(4326))
##plot(st_geometry(occr), pch = 16, col = 'red')

occr <- mutate(occr, GID = 1)


hull <- occr %>% 
  dplyr::group_by(GID) %>% 
  dplyr::summarise() %>%
  st_cast('POLYGON') %>% 
  st_convex_hull()
  
##plot(st_geometry(hull), add = T, col = 'blue', opacity=0.5)
##plot(st_geometry(occr), add = TRUE, col = 'red', pch = 16)

mask_topo <- mask(top_map, mask = hull)
leaflet()%>% addProviderTiles('Esri.WorldImagery')%>%addRasterImage(mask_topo, opacity=0.5)

6. Una vez tenga la interpolación aplicar el modelo a las imágenes para construir una predicción de frutos afectados en la finca y crear una imagen que muestre las zonas con daño mayor a cero. Compare esta imagen con la de los puntos del numeral 2.¿Considera que el modelo representa bien las zonas con daño?, como podría utilizarseeste modelo para hacer control en la finca a futuro?

FrutosAfectados = function(temp,humedad,velocidad,topo){
  return (35.893554+temp*0.076346+velocidad*0.121149-humedad*0.0094961-topo*0.019691)
}

map_pred = FrutosAfectados(mask_vien_map,mask_temp_map,mask_hum_map,mask_topo)





leaflet(datos_frutos) %>% 
  addProviderTiles('Esri.WorldImagery')%>%
  addCircleMarkers(lng = datos_frutos$long,lat = datos_frutos$lat,radius = 0.1, popup = paste0("Altura: ",datos_frutos$altura," ,Frutos Afectados: ",datos_frutos$frutos_afectados),group = "Arboles" ) %>%
  addCircleMarkers(lng = datos$long,lat= datos$lat, radius = 0.1,clusterOptions =datos$frutos_afectados, group ='Cluster')%>%
  addRasterImage(mask_temp_map, opacity=0.8, group ='Temperatura')%>%
  addRasterImage(mask_hum_map, opacity=0.8,group ='Humedad')%>%
  addRasterImage(mask_vien_map, opacity=0.8,group = 'Velocidad del Viento')%>%
  addRasterImage(mask_topo, opacity=0.8,group = 'Altura')%>%
  addRasterImage(map_pred, opacity=0.8,group = 'Predicción frutos afectados')%>%
  addLayersControl(
    overlayGroups = c("Arboles", "Temperatura",'Humedad','Velocidad del Viento', 'Altura','Cluster','Predicción frutos afectados'),
    options = layersControlOptions(collapsed = FALSE)
  )%>%
  addMiniMap(toggleDisplay = TRUE)

Con base al modelo generado previamente se realiza el calculo con las variables climaticas y topográficas, sin embargo según la significancia y basado en el mapa de los cluster hecho previamente que se puede de igual modo ver en este mapa que el modelo no describe del todo bien el fenomeno pues solo dos de cuatro variables son representativas, y solo una es de caracter climatico (temperatura), si este modelo pudiera realizarse con mayor cantidad de datos y con más precisión sería posible utilizarlo de tal modo que se pueda controlar las variables que activan directamente el daño de los frutos, como por ejemplo controlando la temperatura, el riego, entre otros cuidados en la hacienda.

7. Utilice ahora otro enfoque realizando un análisis de clúster con el método kmeans sobre las imágenes de clima interpoladas y los pasos necesarios para identificar y caracterizar las zonas. Genere un mapa de los cluster y compare este con las zonas de mayor daño del punto 2. ¿Considera que existe algún cluster con mayor presencia de daños en los frutos? Que características tiene ese cluster. ¿Como podríamos utilizar este resultado para hacer control en la finca a futuro?

require(raster)
require(rasterVis)
require(ggplot2)

#Preparación de los datos stack

datos=stack(vien_map,temp_map,hum_map,top_map)
plot(datos)

#kmeans

datos_tabla=datos[]


datos_escala=scale(datos_tabla)


datos2=datos

datos2[]=datos_escala


cl=kmeans(datos_escala, centers=10)

#cl$cluster

cluster_map=datos2[[1]]
#plot(cluster_map)
cluster_map[]=cl$cluster
#plot(cluster_map)

##PASO 1- IDENTIFICAR VALOR DE K MAS APROPIADO

calc_desempeño=function(k){
cl=kmeans(datos_escala, centers=k)
desempeño=cl$betweenss/cl$totss*100
return(desempeño)
}

calc_desempeño(k=2)
## [1] 36.76562
ks=2:10
desempeño=sapply(ks, calc_desempeño)

plot(ks, desempeño, type="b")

#nos quedamos con k=5

#Paso2- cuales son los grupos (mapa de cluster)

cl=kmeans(datos_escala, centers=5)

#cl$cluster

cluster_map=datos2[[1]]
plot(cluster_map)

cluster_map[]=cl$cluster
plot(cluster_map)

levelplot(cluster_map)

##PASO 3 - CARACTERIZACIÓN DE LOS GRUPOS

require(table1)

res=cbind.data.frame(datos_tabla, cl$cluster)
table1(~vien_pred.predict
+temp_pred.predict
+hum_pred.predict
+top_pred.predict
|cl$cluster, data=res)
1
(N=943)		 2
(N=2034)		 3
(N=2980)		 4
(N=2255)		 5
(N=1788)		 Overall
(N=10000)
vien_pred.predict
Mean (SD) 1.02 (0.195) 0.432 (0.0605) 0.349 (0.0820) 0.435 (0.135) 0.263 (0.0564) 0.433 (0.224)
Median [Min, Max] 0.998 [0.696, 1.65] 0.425 [0.293, 0.717] 0.336 [0.189, 0.761] 0.415 [0.259, 0.851] 0.258 [0.156, 0.444] 0.366 [0.156, 1.65]
temp_pred.predict
Mean (SD) 25.7 (0.885) 24.4 (1.39) 26.0 (0.890) 27.5 (0.620) 29.4 (0.720) 26.6 (1.91)
Median [Min, Max] 25.7 [23.5, 27.5] 25.0 [21.9, 26.5] 26.1 [23.8, 27.7] 27.6 [25.6, 29.0] 29.3 [27.9, 31.3] 26.6 [21.9, 31.3]
hum_pred.predict
Mean (SD) 34.0 (1.07) 40.4 (2.77) 38.7 (2.08) 33.9 (1.56) 32.4 (1.24) 36.4 (3.66)
Median [Min, Max] 34.1 [31.5, 40.1] 39.8 [33.8, 47.9] 38.3 [35.1, 47.7] 34.1 [30.7, 37.1] 32.2 [29.9, 35.1] 36.1 [29.9, 47.9]
top_pred.predict
Mean (SD) 1900 (4.47) 1880 (3.36) 1890 (4.08) 1890 (3.96) 1890 (2.09) 1890 (6.76)
Median [Min, Max] 1900 [1880, 1910] 1880 [1870, 1890] 1890 [1890, 1900] 1890 [1880, 1900] 1890 [1890, 1900] 1890 [1870, 1910]

Según los resultados se puede evidenciar que si existe un cluster de frutos afectados, el cual se describe de mejor manera en la tabla de resumen de estadisticas el cual resume las variables climaticas y topográficas por la zona del cluster, sin embargo el nivel de información debería detallarse de mejor manera, para que así los dueños de la hacienda tomen decisiones para acabar con los focos de los frutos afectados.