[R-br] Modelo geoestatístico sobre condição, jogo SUDOKU

Wagner Wolff wwolff em usp.br
Segunda Março 19 09:07:53 -03 2018 

Olá pessoal da lista!

Estou lidando com um problema de certa forma parecido com o jogo SUDOKU.
Resumindo nesse jogo tem uma área subdividida em quadriculas ou células que
devem ser preenchidas seguindo condições de que os números de 1 a 9 não se
repitam nas colunas e linhas.

No meu caso, eu preciso interpolar espacialmente uma variável chamada de
intensidade da chuva, entretanto, eu tenho a média dessa variável entre
dois pontos que ligando eles formam uma linha, figura abaixo.

​A questão é a seguinte, como a partir da média eu obtenho os dados
pontuais sobre as linhas que satisfaça a condição de que a média permaneça
até mesmo nos pontos onde as linhas se cruzam. Ou seja, onde elas se cruzam
precisa ser um valor que satisfaça a média da linha 1 e da linha 2.
Depois disso preciso interpolar espacialmente esses pontos na área toda. Eu
obtive uma solução, mas não sei se é a mais adequada e talvez para uma
banco de dados maior ela pode ser impraticável. Resumindo, eu simulei
valores nas linhas com a função grf() do geoR e depois com a função optim()
variei as células (pixels) até que as condições sejam satisfeitas.
Segue o código reproduzível e em anexo os arquivos referente ao shapefile
que representa as linhas.


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
-- 
*Wagner Wolff, **PhD*
"*Luiz de Queiroz**" College of Agriculture,*
University of São Paulo
Pádua Dias avenue11 | 13418-900| Piracicaba-SP| Brazil
Phone:  +55 19 982385582 <+55%2019%2098238-5582>
http://orcid.org/0000-0003-3426-308X
https://github.com/wwolff7
http://buscatextual.cnpq.br/buscatextual/visualizacv.do?id=K4463141A1
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20180319/ef85fa62/attachment.html>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: R Graphics: Device 2 (ACTIVE)_004.png
Tipo: image/png
Tamanho: 20558 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20180319/ef85fa62/attachment.png>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: Links2l.dbf
Tipo: application/x-dbf
Tamanho: 85 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20180319/ef85fa62/attachment.bin>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: Links2l.prj
Tipo: application/octet-stream
Tamanho: 433 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20180319/ef85fa62/attachment.obj>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: Links2l.qpj
Tipo: application/octet-stream
Tamanho: 706 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20180319/ef85fa62/attachment-0001.obj>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: Links2l.shp
Tipo: application/x-esri-shape
Tamanho: 340 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20180319/ef85fa62/attachment-0001.bin>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: Links2l.shx
Tipo: application/x-esri-shape
Tamanho: 116 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20180319/ef85fa62/attachment-0002.bin>

Mais detalhes sobre a lista de discussão R-br