Clase11_DatosFinca
Gabriela Chamorro
23/11/2020
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Mapa de localizacion de ofertas
library(readxl)
fresmo_1_Da <- read_excel("D:/ESPECIALIZACION/SEMESTRE_1/1. Tratamiento de datos/Clase 11/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))+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 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(leaflet)## Loading required package: leafletleaflet() %>% 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-03variograma=variog(geo_dato,option = "bin",uvec = seq(0,0.0009703,0.00003))## variog: computing omnidirectional variogramplot(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 envelopsplot(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 "12" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 161534.404990903model_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 "12" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 128786.888699576model_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 "10.22" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 99030.0417375569lines(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.61415max(geo_dato$coords[,1])## [1] -76.61254min(geo_dato$coords[,2])## [1] 2.380799max(geo_dato$coords[,2])## [1] 2.382694geodatos_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 neighbourhoodpar(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-validationMAE=mean(abs(valida$error))
MAE## [1] 1.10579## Cargamos el ratser
require(raster)## Loading required package: raster## Loading required package: spmat_imagen=data.frame(geodatos_grid,geodatos_ko$predict)
temp_m1_raster=rasterFromXYZ(mat_imagen)
plot(temp_m1_raster)
finca=shapefile("D:/ESPECIALIZACION/SEMESTRE_1/1. Tratamiento de datos/Clase 11/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)
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':
##
## layerlevelplot(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