TP 7

Ejercicio 6

Cargar los paquetes necesarios

require(sp)

## Loading required package: sp

require(spatial)

## Loading required package: spatial

require(gstat)

## Loading required package: gstat

## 
## Attaching package: 'gstat'

## The following object is masked from 'package:spatial':
## 
##     variogram

require(lattice) # este es para graficar

## Loading required package: lattice

Importar objeto

perforaciones <- read.delim("~/Proyectos_R/geoest_TP/TP_7/perforaciones.TXT")

Convertirlo en un objeto de la clase SpatialPointsDataFrame

perf<-SpatialPointsDataFrame(perforaciones[,1:2],data.frame(cota = perforaciones[,3]))

Genero la grilla para graficar la predicción

bbox(perf)

##           min  max
## Easting  15.1 42.0
## Northing  0.1 26.2

expand.grid(Easting=seq(15,42.1, by=0.1), Northing=seq(0,26.3,by=0.1), KEEP.OUT.ATTRS = FALSE)->grid.c
coordinates(grid.c)<-c("Easting","Northing")

Grafico el h-scatterplot

hscat(cota~1,perf, c(0,3,6,9,12) )

Genero el variograma empírico

vario_emp<-variogram( cota ~1, perf)
plot(vario_emp)

Ajusto el variograma teórico (se muestran dos modelos teóricos ajustados)

vario_teo<-fit.variogram(vario_emp,vgm(psill=600000, "Lin",range = 12, nugget = 90000 ))
vario_teo

##   model    psill  range
## 1   Nug      0.0  0.000
## 2   Lin 537447.2 11.282

plot(vario_emp,vario_teo)

plot(vario_emp,vgm(470000, "Sph",range = 11, nugget = 7000))

vario_teo<-fit.variogram(vario_emp,vgm(470000, "Sph",range = 11, nugget = 7000))

## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges
## = fit.ranges, : No convergence after 200 iterations: try different initial
## values?

vario_teo

##   model   psill    range
## 1   Nug       0   0.0000
## 2   Sph 9434659 294.6927

plot(vario_emp,vario_teo) 

ts1 <- krige(cota~Easting+Northing, loc=perf, newdata=grid.c, model=NULL)

## [ordinary or weighted least squares prediction]

levelplot(var1.pred~Easting+Northing, as.data.frame(ts1))

ts2 <- krige(cota~Easting+Northing+I(Easting^2)+I(Northing^2), loc=perf, newdata=grid.c, model=NULL)

## [ordinary or weighted least squares prediction]

levelplot(var1.pred~Easting+Northing, as.data.frame(ts2))

m <- krige(cota~Easting+Northing, loc=perf, newdata=grid.c, model=vario_teo) 

## [using universal kriging]

levelplot(var1.pred~Easting+Northing, as.data.frame(m))

El último modelo con los detalles gráficos

levelplot(var1.pred~Easting+Northing, as.data.frame(m),
          col.regions=bpy.colors(64),
          asp="iso", main="cota",
          sub="sitios muestreados",
          xlab="Este", ylab="Norte",
          panel=function(x,y,z, ...) {
            panel.levelplot(x, y, z, ...)
            panel.grid(h=-1,v=-1, col="gray", lty=1)
            panel.points(coordinates(perf),
                         cex=1*perf$cota/max(perf$cota),
                         pch=20, col="white")
          })