Krigagem Ordinária no R

O código abaixo exemplifica a utilização da Krigagem Ordinária utilizando o R dentro de uma sessão do GRASS. O intuito desta postagem é mostrar a codificação básica necessária para se chegar a um primeiro resultado e não será abordado aqui como realizar o ajuste do variograma.

	#No terminal em que o GRASS está em execução, abre-se o R:
	R
	#Bibliotecas necessárias:
	library(gstat)
	library(spgrass6)
	#Leitura das variáveis de ambiente do GRASS:
	G <- gmeta6()
	#Cópia da camada vetorial na variável abaixo (pontos a serem interpolados):
	pcs <- readVECT6('pontos_cotados')
	#Objeto grid baseado na region ativa:
	grid <- gmeta2grd()
	#Objeto new_data, necessário para realizar a krigagem:
	new_data <- SpatialGridDataFrame(grid, data=data.frame(k=rep(1,G$cols*G$rows)), proj4string=CRS(pcs@proj4string@projargs))
	#Criação do objeto gstat:
	g <- gstat(id="elev", formula=elev ~ 1, data=pcs)
	#Ajuste do variograma, utilizando os parâmetros default:
	v.fit <- fit.variogram(variogram(g) ,vgm(model="Lin") )
	#Atualização do objeto gstat, com o variograma ajustado:
	g <- gstat(g, id="elev", model=v.fit )
	#Krigagem ordinária:
	pcs.krig <- predict(g, id='elev', newdata=new_data)
	#Criação da camada raster no GRASS:
	writeRAST6(pcs.krig, 'pcs.krig', zcol='elev.pred')
	#saída do R:
	q()

view raw exemplo_krigagem.R hosted with ❤ by GitHub

Foram utilizados 300 pontos cotados para gerar a superfície mostrada na Figura 1.

Figura 1 - Resultado da Krigagem