Geo Datos Clima
David Arango
4/19/2021
Importar Datos y Librerias
library(readxl)
library(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
## --------------------------------------------------------------library(raster)## Loading required package: splibrary(ggplot2)
library(leaflet)
geo_datos =read_excel("~/Desktop/geo_datos.xlsx")
geo_datos## # A tibble: 394 x 7
## id_arbol Latitude Longitude `Relative Humidity` Temperature `Wind Speed`
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 2.38 -76.6 61.2 22.9 0
## 2 2 2.38 -76.6 62.3 23.8 0
## 3 3 2.38 -76.6 62.6 22.6 0.4
## 4 4 2.38 -76.6 61.3 22.3 0.6
## 5 5 2.38 -76.6 61.9 22.7 0.8
## 6 6 2.38 -76.6 65.7 23.1 0.4
## 7 7 2.38 -76.6 67 21 0.5
## 8 8 2.38 -76.6 65.8 21.6 0.4
## 9 9 2.38 -76.6 66.6 21.4 0
## 10 10 2.38 -76.6 66.3 20.8 0.7
## # … with 384 more rows, and 1 more variable: Altitude <dbl>Se observa que los datos contienen información sobre las localizaciones de los sitios donde se tomaron datos de clima como temperatura, humedad relativa entre otras.
leaflet() %>% addTiles() %>% addCircleMarkers(lng = geo_datos$Longitude,lat = geo_datos$Latitude,radius = 0.2,color = "black")# leaflet() %>% addTiles() %>% addCircleMarkers(lng = geo_datos$Longitude,lat = geo_datos$Latitude,radius = 0.2,color = "black",clusterOptions = markerClusterOptions())Se observa que la finca se encuentra en la zona rural del municipio de Popayan y su distribucción espacial es bastante regular.
A. GeoestadÃstica
Como primer paso realizamos la creación de la variable regionalizada:
geo_temp=as.geodata(geo_datos,coords.col = 3:2,data.col = 5)
plot(geo_temp)
summary(dist(geo_datos[,3:2]))## 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_temp,option = "bin",uvec=seq(0,0.0009703,9.703e-05))## variog: computing omnidirectional variogramplot(variograma)
variograma_mc=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 envelopslines(variograma_mc)
Se observa en el semivariograma una autocorrelación espacial muy fuerte indicando que las mediciones cercanas, presentan temperaturas similares en contraste con las distantes.
Ahora vamos a realizar el ajuste del modelo teorico:
##Ajuste del Modelo de Semivarianza
ini.vals = expand.grid(seq(8,12,l=10), seq(6e-04,8e-04,l=10))
plot(variograma)
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: 174964.991544732model_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: 130506.318494236model_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: 107137.155200497lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="purple")
Podemos observar tanto graficamente como por el valor de SCE que el modelo gaussiano es el que mejor se ajusta a los puntos del semivariograma muestral.
Ahora vamos a realizar la interpolación (predicción espacial)
max(geo_datos[,3])## [1] -76.61254##Predicción Espacial Kriging
geodatos_grid=expand.grid( lon=seq(-76.61415,-76.61254,l=100),lat=seq(2.380799 ,2.382694 ,l=100))
plot(geodatos_grid)
points(geo_datos[,3:2],col="red")
geodatos_ko=krige.conv(geo_temp, loc=geodatos_grid,
krige= krige.control(nugget=0,trend.d="cte",
trend.l="cte",cov.pars=c(sigmasq=12.7080, phi=0.0007 )))## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhoodpar(mfrow=c(1,2))
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)
Ahora vamos a tranformar la imagen en raster:
require(rasterVis)## Loading required package: rasterVis## Loading required package: terra## terra version 1.1.4## Loading required package: lattice## Loading required package: latticeExtra##
## Attaching package: 'latticeExtra'## The following object is masked from 'package:ggplot2':
##
## layerpred=cbind(geodatos_grid,geodatos_ko$predict)
temp_predict=rasterFromXYZ(cbind(geodatos_grid,geodatos_ko$predict))
plot(temp_predict)
levelplot(temp_predict,par.settings =BuRdTheme)
temp_error=rasterFromXYZ(cbind(geodatos_grid,sqrt(geodatos_ko$krige.var)))
levelplot(temp_error,par.settings =BuRdTheme)
por ultimo vamos a cortar los raster con el shape de la finca:
finca=shapefile("~/Desktop/Fresno_area.shp")
plot(temp_predict)
plot(finca,add=T)
temp_predict_finca=mask(temp_predict,finca)
plot(temp_predict_finca)
plot(finca,add=T)