Interpolacion
Andres Lopez
19/4/2021
El presente informe muestra la metodología empleada para obtener un modelo digital de elevación “MDE” mediante técnicas de geoestadística haciendo uso del software RStudio.
Primero, cargamos la librerías necesarias.
-(readxl): Para la carga de tablas de Excel. -(geoR): Para tratamiento de datos espaciales en Geoestadística. -(raster): Para la manipulación de información en formato RASTER. -(ggplot): Para las gráficas. -(Leaflet): Para creación de aplicaciones que pueden ser representadas en páginas HTML. -(rasterVis): Para mejorar la visualización e interacción con archivos raster.
Haciendo uso de la librería (readxl) Cargamos la tabla Excel “geo_datos” y visualizamos su contenido:
getwd()## [1] "C:/Users/Andres/Desktop/ESPECIALIZACION UNIVALLE/1. TRATAMIENTO DE DATOS ESPACIALES - LUNES 7-10 AM/CLASE 10 - 19-04-2021"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)
require(rasterVis)## Loading required package: rasterVis## Warning: package 'rasterVis' was built under R version 4.0.5## Loading required package: terra## Warning: package 'terra' was built under R version 4.0.5## terra version 1.1.4## Loading required package: lattice## Loading required package: latticeExtra## Warning: package 'latticeExtra' was built under R version 4.0.5##
## Attaching package: 'latticeExtra'## The following object is masked from 'package:ggplot2':
##
## layer##Cargamos la tabla excel "geo_datos" y visualizamos su contenido:
geo_datos = read_excel("geo_datos.xlsx")
geo_datos## # A tibble: 394 x 7
## id_arbol Latitude Longitude `Relative Humid~ 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 de localización (coordenadas geográficas "Latitud - Longitud) que corresponde a 394 árboles, a los cuales se les asocia variables climáticas como, (4) humedad relativa, (5) temperatura, (6) velocidad del viento y (7) altitud.
Haciendo uso de la librería (leaflet) graficamos los puntos en sistema de referencia geográfico:
leaflet() %>% addTiles() %>% addCircleMarkers(lng = geo_datos$Longitude,lat = geo_datos$Latitude,radius = 0.2,color = "black")# para prender polígonos de Thiessen utilizamos "clusterOptions = markerClusterOptions()":
# 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 zona rural al sur del municipio de Popayán, se aprecia que la distribución espacial de los árboles es bastante homogénea.
##COMENZAMOS CON EL ANÁLISIS GEOESTADÍSTICO##
Como primer paso realizamos la creación de la variable regionalizada (columna 7. Altitud), esperando que los resultados sean satisfactorios debido a que la topografía o elevación del terreno tiene una correlacion espacial muy fuerte:
as.geodata se utiliza para poner datos como variable regionalizada, con el fin de que pueda ser utilizada para hacer predicciones, donde, geo_datos es la tabla de Excel importada, las columnas de las coordenadas son la 3 (long) y la 2 (lat) y el dato de la variable regionalizada está en la columna 7 (Altitud).
geo_alt=as.geodata(geo_datos,coords.col = 3:2,data.col = 7)
plot(geo_alt)
La gráfica anterior muestra un histograma bi-modal con alturas predominantes de 1.87-190 y 1.93-1.94, además, muestra tendencias de niveles de terreno homogéneas por sectores.
Sacamos la matriz de distancias para verificar en que rango se va a mover el variograma y dividimos el valor máximo (3er cuartil) sobre 10 para crear 10 puntos en la representación.
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_alt,option = "bin",uvec=seq(0,0.0009703,9.703e-05))## variog: computing omnidirectional variogramUtilizamos la simulación de Montecarlo, para permutar los datos de forma aleatoria y con esto romper cualquier correlación espacial.
La finalidad de ejecutar este paso es darnos cuenta si la variable tiene correlación espacial (con esto obtenemos las bandas de no correlación), en el caso que el espacio comprendido entre el variograma y las bandas sea el mismo o sea muy parecido, descartaría este método para interpolar y hacer predicciones por la no correlación espacial en la variable analizada.
variograma_mc=variog.mc.env(geo_alt,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(variograma_mc)
La grafica anterior demuestra que existe una correlación espacial muy fuerte en la variable analizada (Altitud), el patrón es muy diferente a las bandas de no correlación, indicando que las mediciones cercanas, presentan alturas similares en contraste con las más lejanas.
Realizamos el ajuste del modelo teórico de Semivarianza
Utilizamos el rango de valores del alcance (de 7e-04 a 11e-04) y la meseta (de 6e-04 a 8e-04) analizando la gráfica anterior, creamos una grilla de 10 y posteriormente corremos los modelos Exponencial, Gaussiano y Esférico.
ini.vals = expand.grid(seq(7e-04,11e-04,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 "0" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 0.00231708894611235model_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 "0" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 0.000625213373359565model_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 "0" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 0.00409035977432688lines(model_mco_exp,col="green")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="purple")
Los valores Sigma y Phi para cada uno de los modelos se presentan a continuación:
Exponential: SCE: 0.0023 - Sigma:0.0009 - Phi: 0.0008 Gauss: SCE: 0.0006 - Sigma:0.0011 - Phi: 0.0008 Spherical: SCE: 0.0041 - Sigma:0.0007 - Phi: 0.0008
Podemos observar tanto gráficamente como por el valor de SCE (sumas de cuadrados del error minima ponderada) 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)
Empezamos sacando los mínimos y máximos de las coordenadas para crear el Data Frame con una grilla de 100 para tener una predicción algo detallada y luego corremos la interpolación con Kriging:
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_alt, loc=geodatos_grid,
krige= krige.control(nugget=0,trend.d="cte",
trend.l="cte",cov.pars=c(sigmasq=0.0007, phi=0.0008 )))## 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)
En la gráfica anterior (Varianza de predicción) se puede ver que la precisión del modelo en cuanto a alturas es más fuerte entre más cerca se esté de las posiciones de los árboles y menor entre más se aleje de ellos, por fuera del polígono que encierra todos los árboles, entre más lejos, los valores se van tornando menos confiables, por esto, es importante tomar registros en grilla sobre toda el área de estudio para poder interpolar de una mejor manera.
Ahora vamos a transformar la imagen en raster
geodatos_grid = Coordenadas de los 10000 ptos. geodatos_ko$predict = Predicción para cada punto cbind = Crea una matriz con las coordenadas y los valores predichos de altitud para cada pixel. rasterFromXYZ = Crea una imagen raster con la matriz creada con cbind.
pred=cbind(geodatos_grid,geodatos_ko$predict)
alt_predict=rasterFromXYZ(cbind(geodatos_grid,geodatos_ko$predict))
plot(alt_predict)
levelplot(alt_predict,par.settings =BuRdTheme)
alt_error=rasterFromXYZ(cbind(geodatos_grid,sqrt(geodatos_ko$krige.var)))
levelplot(alt_error,par.settings =BuRdTheme)
Cortamos el ráster generado con la interpolación de Kriging para las alturas con el Shape (polígono de la finca) asociándole el sistema de referencia geográfico.
finca=shapefile("Fresno_area.shp")
plot(alt_predict)
plot(finca,add=T)
alt_predict_finca=mask(alt_predict,finca)
plot(alt_predict_finca)
plot(finca,add=T)
crs(alt_predict_finca)=crs("+proj=longlat +datum=WGS84")
leaflet() %>% addTiles() %>% addRasterImage(alt_predict_finca,opacity = 0.6)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 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 +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 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 +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definitionFinalmente realizamos la validación cruzada (quita un valor, lo reemplaza con el predicho y lo compara con el valor real)
valida_exp=xvalid(geodata = geo_alt,model = model_mco_exp) MAE=mean(abs(valida_exp$error)) MAE
Ocurrió un error al procesar con el modelo Gaussiano (insertando los valores de Sigma y Phi en la interpolación con Kriging), por consiguiente, usamos el modelo exponencial que presentó la segunda SCE más baja.
Y creamos la imagen raster “alt.tif”.
writeRaster(alt_predict_finca, "alt.tif", overwrite=TRUE)Que corresponde a un modelo digital de elevación MDE, interpolado por el método de Kriging, cortado con el perímetro del área de estudio, teniendo como base la localización y altitud de 394 árboles en una finca al sur de la ciudad de Popayán - Departamento del Valle del Cauca.