proj="+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=km +no_defs "

## Colocar os arquivis anexados na mesma pasta da código R

Links2l <- spTransform(readOGR("Links2l.shp"), CRS(proj))
Rmean <- c(2.9424, 4.3982)

## Criando uma imagem raster vazia
ras <- raster(Links2l, res=0.1, ext=extent(Links2l)); ras[] <- rep(NA, ncell(ras))

## Celulas (pixels) interceptadas pelas links
cells_links <- cellFromLine(ras, Links2l)
cells <- unique(unlist(cells_links))

## Celulas compartilhadas entre os links
cells_links[[1]][which(cells_links[[1]] %in% cells_links[[2]])]

## Simulando valores a partir da media nas celulas interceptadas pelas links
set.seed(234)
for(i in 1:nrow(Links_ams)){
    ras[cells_links[[i]]] <- apply(grf(length(cells_links[[i]]), grid=xyFromCell(ras, cells_links[[i]]), nsim=10,
                                       mean = Rmean[i],cov.pars = c(3.7, 18.7), nugget = 0.37)$data, 1, mean)
}

plot(ras, col = gray.colors(10, start = 0.3, end = 0.9, gamma = 2.2, alpha = NULL))
lines(Links2l)

## Funcao de otimizacao
opt <- function(par, Rmean){
    res <- c()
    for(i in 1:length(Rmean)){
        res[i] <- sum(mean(par[cells %in% cells_links[[i]]])-Rmean[i])^2
    }
    sum(res)
}

## Definir espaco de iteracao
lower <- list()
upper <- list()

for(i in 1:nrow(Links_ams)){
    lower[[i]] <- rep(mean(ras[cells_links[[i]]])-3.7^0.5, ncell(cells_links[[i]]))
    upper[[i]] <- rep(mean(ras[cells_links[[i]]])+3.7^0.5, ncell(cells_links[[i]]))
}

lower <- ifelse(unlist(lower)<0, 0, unlist(lower))
upper <- ifelse(unlist(upper)<0, mean(unlist(upper))+sd(unlist(upper)), unlist(upper))

## Otimizar e atualizar as celulas para resolver a condicao
cell_upd <- optim(par=ras[cells], fn=opt,lower = lower, upper=upper, Rmean=Rmean,method="L-BFGS-B")$par

ras_old <- ras
ras[cells] <- cell_upd

## Comparando celulas antigas com as novas, apenas a célula compartilhada entre os links nao é linear
plot(ras_old[cells_links[[1]]], ras[cells_links[[1]]], ann=T)
abline(a=0, b=1)

plot(ras_old[cells_links[[2]]], ras[cells_links[[2]]], ann=T)
abline(a=0, b=1)

## Verificando se as condicoes foram satisfeitas
mean(ras[cells_links[[1]]])
Rmean[1]

mean(ras[cells_links[[2]]])
Rmean[2]

## Resultado, valores iguais "SUDOKU" resolvido!

######################################################
## Interpolacao das celulas

pts <- rasterToPoints(ras)
plot(pts)
boxcox(pts~1)

geo_link <- as.geodata(pts, coords.col = 1:2, data.col = 3)

plot(geo_link, low =T)
plot(geo_link, trend='1st', low =T)
plot(geo_link, trend='2nd', low =T)

vario <- variog(geo_link,max.dist=10,uvec=seq(0, 10, by=0.01))
plot(vario)

link_fit <- list()
link_fit$cte <- likfit(geo_link, cov.model="exp", trend ="cte",ini=c(0.5,5), nugget = 0.1)
link_fit$lon <- likfit(geo_link, cov.model="exp", trend =geo_link$coords[, 1],ini=c(0.5,5), nug=0.1)
link_fit$st <- likfit(geo_link, cov.model="exp", trend ="1st",ini=c(0.5,5), nug=0.1)
link_fit$nd <- likfit(geo_link, cov.model="exp", trend ="2nd",ini=c(0.5,20), nug=0.1)

sapply(link_fit, AIC)

summary(link_fit$st)

## Predição na área krigagem
grid <- xyFromCell(ras, 1:ncell(ras))
link_fit$st
kr.link <- krige.conv(geo_link, loc=grid, krige=krige.control(obj=link_fit$st))

ras[] <- kr.link$predict

image(ras)
contour(ras, add=T)
lines(Links2l, col=3, lwd=2)

## Krigagem ordinaria é BLUE, consequentemente condicoes sao satisfeitas
mean(ras[cells_links[[1]]])
Rmean[1]

mean(ras[cells_links[[2]]])
Rmean[2]

Att.

Abraco

--