library(readxl)
Fresno_1_Da <- read_excel("C:/Users/Andres/Desktop/ESPECIALIZACION/1ER SEMESTRE/TRATAMIENTO DE DATOS ESPACIALES/CLASE10/Fresno_1_Da.xlsx",
sheet = "Sheet1")
View(Fresno_1_Da)
plot (Fresno_1_Da [,3:2])
require (ggplot2)
## Loading required package: ggplot2
ggplot (Fresno_1_Da, aes(x= Longitude, y= Latitude))+geom_point()+theme_bw()
ggplot (Fresno_1_Da, aes(x= Longitude, y= Latitude))+geom_point()+theme_bw()+facet_wrap(~Date2)
ggplot (Fresno_1_Da, aes(x= Longitude, y= Latitude, color= Temperature))+geom_point()+theme_bw()+facet_wrap(~Date2) + scale_color_continuous(type= "viridis")
### Filtrar un instante del tiempo
fechas=unique(Fresno_1_Da$Date2)
pos=which(Fresno_1_Da$Date2==fechas[1])
datos_m1=Fresno_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 21/0~ 14.8 805.
## 2 2 2.38 -76.6 21/0~ 11.6 805.
## 3 3 2.38 -76.6 21/0~ 12.9 806.
## 4 4 2.38 -76.6 21/0~ 14.1 806.
## 5 5 2.38 -76.6 21/0~ 14.3 805.
## 6 6 2.38 -76.6 21/0~ 14.2 805.
## 7 7 2.38 -76.6 21/0~ 14.4 805.
## 8 8 2.38 -76.6 21/0~ 12.8 805.
## 9 9 2.38 -76.6 21/0~ 15 805
## 10 10 2.38 -76.6 21/0~ 14 805.
## # ... with 384 more rows, and 3 more variables: `Relative Humidity` <dbl>,
## # Crosswind <dbl>, Temperature <dbl>
require(leaflet)
## Loading required package: leaflet
require(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)
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)
variograma.env=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.env)
##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: 24954.5072627246
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: 135574.864752597
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: 209253.508136788
lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="purple")
##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=1.5665, phi=0.0004)))
## 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="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)
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] 1.348928
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("C:/Users/Andres/Desktop/ESPECIALIZACION/1ER SEMESTRE/TRATAMIENTO DE DATOS ESPACIALES/CLASE10/Fresno_area.shp")
plot(finca, add=TRUE)
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)
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