Krigagem

Atividade da Especialização em Geoprocessamento e Georreferênciamento Neste projeto, foram usados dados de estações meterologicas do estado de Minas Gerais/BR. Onde então estudou-se a temperatura superficial da mesma.

Autor: João Ataíde; Portifolio: https://www.joaoataide.com; Github: https://github.com/jvataidee/atividade_Geoestatistica

Instalando Packeges

install.packages("geoR")
install.packages("splancs")
install.packages("e1071")

Importando Packeges

library(geoR)
library(e1071)
library(splancs)

Importando dados (Escolha dados Anuais)

dados = read.geodata(file =  "anual.txt", sep = ' ', header = T)

dados

$coords
      Coordenada_X Coordenada_Y
 [1,]    1779363.1      1866667
 [2,]    1917183.0      1743463
 [3,]    1644171.3      1691338
 [4,]    1824790.9      1653488
 [5,]    1334213.4      1663299
 [6,]    1732933.9      1549995
 [7,]    1866684.0      1534184
 [8,]    1048687.8      1515874
 [9,]    1355469.9      1500911
[10,]    1652439.3      1512304
[11,]    1811659.5      1424039
[12,]    1500722.7      1414913
[13,]    1204360.5      1379653
[14,]    1304633.6      1392301
[15,]    1754715.3      1878732
[16,]    1394579.9      1334039
[17,]    1572480.5      1378406
[18,]    1598646.2      1321875
[19,]    1659615.0      1420473
[20,]    1888319.0      1340905
[21,]    1585652.0      1312729
[22,]    1794456.3      1230586
[23,]    1694857.4      1215287
[24,]    1382474.5      1139880
[25,]    1477593.9      1182225
[26,]    1598063.9      1167414
[27,]    1414037.7      1565444
[28,]    1632604.6      1102693
[29,]    1465614.1      1081517
[30,]    1776723.5      1317055
[31,]    1239981.1      1838058
[32,]    1303890.6      1861736
[33,]     821383.9      1285369
[34,]    1245387.0      1282935
[35,]    1865499.2      1183149
[36,]    1224628.0      1175973
[37,]    1784017.6      1152091
[38,]    1404004.5      1013691
[39,]    1187204.0      1115865
[40,]    1664009.6      2004230
[41,]    2032323.2      1612207
[42,]    1594320.6      1041196

$data
 [1] 1367.734 1173.866 1278.491 1167.419 1181.332 1001.410 1027.926 1228.853 1172.484 1062.085 1083.234 1101.362 1190.618 1098.009 1510.332
[16] 1042.998 1156.135 1144.668 1018.784 1257.934 1067.023 1010.569  993.191 1011.884 1049.461  937.951 1249.710  962.760  964.234 1037.695
[31] 1199.366 1163.055 1174.166 1151.219 1136.087 1205.206 1110.348  798.245 1198.475 1424.981 1054.153 1013.840

$call
read.geodata(file = "anual.txt", header = T, sep = " ")

attr(,"class")
[1] "geodata"

summary(dados)

Number of data points: 42 

Coordinates summary
    Coordenada_X Coordenada_Y
min     821383.9      1013691
max    2032323.2      2004230

Distance summary
       min        max 
  15890.42 1254271.48 

Data summary
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
 798.245 1030.368 1123.217 1123.316 1188.296 1510.332 

Other elements in the geodata object
[1] "call"

Estátisticas Básicas

names(dados)

$coords
[1] "Coordenada_X" "Coordenada_Y"

$data
[1] "data"

$other
[1] "call"

sd(dados$data)

[1] 132.8226

var(dados$data)

[1] 17641.83

Moda

getmode = function(v) {
  uniqv = unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}
moda = getmode(dados$data)
moda

[1] 1367.734

skewness(dados$data)

[1] 0.5129477

kurtosis(dados$data)

[1] 0.8814792

Histograma

hist(dados$data, main = NULL, xlab = 'Data', ylab = 'Frequência')

Histograma com curva normal

dif = diffpairs(dados$coords, dados$data)
hist(dif$dif, xlab = 'Classe de ET0', ylab = 'Frequência',, main = NULL, prob = T)
dados.d = density(dif$diff)
lines(dados.d, col = "red")

Boxplot

boxplot(dados$data, horizontal=TRUE)

Shapiro test (TESTE DE NORMALIDADE)

plot(dados)

shapiro.test(dados$data)

    Shapiro-Wilk normality test

data:  dados$data
W = 0.96282, p-value = 0.1862

plot(dados$data, dados$coords[, 1], xlab = 'ET0', ylab = 'Latitude')

plot(dados$data,dados$coords[,2], xlab = 'ET0', ylab = 'Latitude')

Variogáfica Semivariograma empírico ou experimental (criar outro objeto)

binET0 = variog(dados,uvec = 12,max.dist = 1000000, pairs.min = 10)

variog: computing omnidirectional variogram

Semivariograma modelado máxima verossimilhança (ajuste do modelo)

ml = likfit(dados,  ini = c(17000, 267000), nugget = 0, cov.model = 'exponential')

kappa not used for the exponential correlation function
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
         arguments for the maximisation function.
        For further details see documentation for optim.
likfit: It is highly advisable to run this function several
        times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.

ml1 = likfit(dados, ini = c(17000, 267000), nugget = 0, cov.model = 'spherical')

kappa not used for the spherical correlation function
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
         arguments for the maximisation function.
        For further details see documentation for optim.
likfit: It is highly advisable to run this function several
        times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.

ml2 = likfit(dados, ini = c(17000, 267000), nugget = 0, cov.model = 'gauss')

kappa not used for the gaussian correlation function
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
         arguments for the maximisation function.
        For further details see documentation for optim.
likfit: It is highly advisable to run this function several
        times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.

ml

likfit: estimated model parameters:
      beta      tausq    sigmasq        phi 
"  1164.7" "   690.1" " 15566.6" "267000.0" 
Practical Range with cor=0.05 for asymptotic range: 799860.5

likfit: maximised log-likelihood = -249.5

plot(binET0)
lines(ml, type = 'l', lty = 2, lwd = 1,  col = "red")

lines(ml1, type = 'l', lty = 2, lwd = 1,  col = "blue",)
lines(ml2, type = 'l', lty = 2, lwd = 1,  col = "green")

legend("topleft", legend = c("exponencial", "esferica","gaussiana" ),
       col = c("red", "blue", "green"), lty=2, cex=0.8)

Semivariograma modelado máxima verossimilhança (ajuste do modelo)

rl = variofit (binET0, ini = c(17000, 267000), nugget = 0, cov.model = 'exponential')

variofit: covariance model used is exponential 
variofit: weights used: npairs 
variofit: minimisation function used: optim 

rl1 = variofit (binET0, ini = c(17000, 267000), nugget = 0, cov.model = 'spherical')

variofit: covariance model used is spherical 
variofit: weights used: npairs 
variofit: minimisation function used: optim 

rl2 = variofit (binET0, ini = c(17000, 267000), nugget = 0, cov.model = 'gauss')

variofit: covariance model used is gaussian 
variofit: weights used: npairs 
variofit: minimisation function used: optim 

plot(binET0)
lines(rl, type = 'l', lty = 2, lwd = 1,  col = "red",)

lines(rl1, type = 'l', lty = 2, lwd = 1,  col = "blue",)
lines(rl2, type = 'l', lty = 2, lwd = 1,  col = "green")

legend("topleft", legend = c("exponencial", "esferica","gaussiana" ),
       col = c("red", "blue", "green"), lty=2, cex=0.8)

Validação cruzada (estima os valores e compara com real)

valid_ml = xvalid(dados, model = ml)

xvalid: number of data locations       = 42
xvalid: number of validation locations = 42
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, 
xvalid: end of cross-validation

valid_ml1 = xvalid(dados, model = ml1)

xvalid: number of data locations       = 42
xvalid: number of validation locations = 42
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, 
xvalid: end of cross-validation

valid_ml2 = xvalid(dados, model = ml2)

xvalid: number of data locations       = 42
xvalid: number of validation locations = 42
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, 
xvalid: end of cross-validation

valid_rl = xvalid(dados, model = rl)

xvalid: number of data locations       = 42
xvalid: number of validation locations = 42
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, 
xvalid: end of cross-validation

valid_rl1 = xvalid(dados, model = rl1)

xvalid: number of data locations       = 42
xvalid: number of validation locations = 42
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, 
xvalid: end of cross-validation

valid_rl2 = xvalid(dados, model = rl2)

xvalid: number of data locations       = 42
xvalid: number of validation locations = 42
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, 
xvalid: end of cross-validation

Exportando os dados da Validação

write.csv(valid_ml$data, file = "validação1")
write.csv(valid_ml1$data, file = "validação2")
write.csv(valid_ml2$data, file = "validação3")
write.csv(valid_rl$data, file = "validação11")
write.csv(valid_rl1$data, file = "validação21")
write.csv(valid_rl2$data, file = "validação31")

Krigagem novo pacote

Importação do Poligono

limite = read.table(file = "areaminas.txt", header = T)
limite

Limites do estado

plot(limite)

loci = expand.grid( seq( 800000, 2500000, l = 100), seq(800000, 2000000, l = 100))

Modelo escolhido foio gaussiano

kc = krige.conv(dados, loc = loci, border = limite, krige = krige.control(obj = ml))

krige.conv: results will be returned only for prediction locations inside the borders
krige.conv: model with constant mean
krige.conv: Kriging performed using global neighbourhood 

Mapa de krigagem (curva de nível)

contour( kc, xlab = 'Longitude', ylab = 'Latitude')
title(main = "Temperatura Superficial", font.main = 4)

Mapa de krigagem

image(kc, loc = loci, border = limite, val = kc$predict, xlab= 'Longitude', ylab='Latitude')
title(main = "Temperatura Superficial", font.main = 4)