<div dir="ltr"><div><div><div><div><div>Olá pessoal da lista!<br><br></div>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. <br><br></div>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. <br><img src="cid:ii_jey65rey0_1623e1949d5c5b87" style="margin-right: 0px;" width="171" height="172"><br></div>​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.<br></div>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. <br></div>Segue o código reproduzível e em anexo os arquivos referente ao shapefile que representa as linhas.<br><div><div><div><div><div><div><div><div><br><br>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 "<br></div><div>## Colocar os arquivis anexados na mesma pasta da código R<br></div><div>Links2l <- spTransform(readOGR("Links2l.shp"), CRS(proj))<br>Rmean <- c(2.9424, 4.3982)<br><br>## Criando uma imagem raster vazia<br>ras <- raster(Links2l, res=0.1, ext=extent(Links2l)); ras[] <- rep(NA, ncell(ras))<br><br>## Celulas (pixels) interceptadas pelas links<br>cells_links <- cellFromLine(ras, Links2l)<br>cells <- unique(unlist(cells_links))<br><br>## Celulas compartilhadas entre os links<br>cells_links[[1]][which(cells_links[[1]] %in% cells_links[[2]])]<br><br>## Simulando valores a partir da media nas celulas interceptadas pelas links<br>set.seed(234)<br>for(i in 1:nrow(Links_ams)){<br>    ras[cells_links[[i]]] <- apply(grf(length(cells_links[[i]]), grid=xyFromCell(ras, cells_links[[i]]), nsim=10,<br>                                       mean = Rmean[i],cov.pars = c(3.7, 18.7), nugget = 0.37)$data, 1, mean)<br>}<br><br>plot(ras, col = gray.colors(10, start = 0.3, end = 0.9, gamma = 2.2, alpha = NULL))<br>lines(Links2l)<br><br>## Funcao de otimizacao<br>opt <- function(par, Rmean){<br>    res <- c()<br>    for(i in 1:length(Rmean)){<br>        res[i] <- sum(mean(par[cells %in% cells_links[[i]]])-Rmean[i])^2<br>    }<br>    sum(res)<br>}<br><br>## Definir espaco de iteracao<br>lower <- list()<br>upper <- list()<br><br>for(i in 1:nrow(Links_ams)){<br>    lower[[i]] <- rep(mean(ras[cells_links[[i]]])-3.7^0.5, ncell(cells_links[[i]]))<br>    upper[[i]] <- rep(mean(ras[cells_links[[i]]])+3.7^0.5, ncell(cells_links[[i]]))<br>}<br><br>lower <- ifelse(unlist(lower)<0, 0, unlist(lower))<br>upper <- ifelse(unlist(upper)<0, mean(unlist(upper))+sd(unlist(upper)), unlist(upper))<br><br>## Otimizar e atualizar as celulas para resolver a condicao<br>cell_upd <- optim(par=ras[cells], fn=opt,lower = lower, upper=upper, Rmean=Rmean,method="L-BFGS-B")$par<br><br>ras_old <- ras<br>ras[cells] <- cell_upd<br><br>## Comparando celulas antigas com as novas, apenas a célula compartilhada entre os links nao é linear<br>plot(ras_old[cells_links[[1]]], ras[cells_links[[1]]], ann=T)<br>abline(a=0, b=1)<br><br>plot(ras_old[cells_links[[2]]], ras[cells_links[[2]]], ann=T)<br>abline(a=0, b=1)<br><br>## Verificando se as condicoes foram satisfeitas<br>mean(ras[cells_links[[1]]])<br>Rmean[1]<br><br>mean(ras[cells_links[[2]]])<br>Rmean[2]<br><br>## Resultado, valores iguais "SUDOKU" resolvido!<br><br>######################################################<br>## Interpolacao das celulas<br><br>pts <- rasterToPoints(ras)<br>plot(pts)<br>boxcox(pts~1)<br><br>geo_link <- as.geodata(pts, coords.col = 1:2, data.col = 3)<br><br>plot(geo_link, low =T)<br>plot(geo_link, trend='1st', low =T)<br>plot(geo_link, trend='2nd', low =T)<br><br>vario <- variog(geo_link,max.dist=10,uvec=seq(0, 10, by=0.01))<br>plot(vario)<br><br>link_fit <- list()<br>link_fit$cte <- likfit(geo_link, cov.model="exp", trend ="cte",ini=c(0.5,5), nugget = 0.1)<br>link_fit$lon <- likfit(geo_link, cov.model="exp", trend =geo_link$coords[, 1],ini=c(0.5,5), nug=0.1)<br>link_fit$st <- likfit(geo_link, cov.model="exp", trend ="1st",ini=c(0.5,5), nug=0.1)<br>link_fit$nd <- likfit(geo_link, cov.model="exp", trend ="2nd",ini=c(0.5,20), nug=0.1)<br><br>sapply(link_fit, AIC)<br><br>summary(link_fit$st)<br><br>## Predição na área krigagem<br>grid <- xyFromCell(ras, 1:ncell(ras))<br>link_fit$st<br>kr.link <- krige.conv(geo_link, loc=grid, krige=krige.control(obj=link_fit$st))<br><br>ras[] <- kr.link$predict<br><br>image(ras)<br>contour(ras, add=T)<br>lines(Links2l, col=3, lwd=2)<br><br>## Krigagem ordinaria é BLUE, consequentemente condicoes sao satisfeitas<br>mean(ras[cells_links[[1]]])<br>Rmean[1]<br><br>mean(ras[cells_links[[2]]])<br>Rmean[2]<br><br clear="all"></div><div>Att. <br></div><div>Abraco<br></div><div>-- <br><div class="gmail_signature"><div dir="ltr"><span style="font-size:medium"><b><i>Wagner Wolff, </i></b></span><i><b>PhD</b></i><br>"<b>Luiz de Queiroz</b><b><span>"</span> College of Agriculture,</b><br>University of São Paulo<br>Pádua Dias avenue11 | 13418-900| Piracicaba-SP| Brazil<br>Phone:  <a href="tel:+55%2019%2098238-5582" value="+5519982385582" target="_blank">+55 19 982385582</a>  <br><span><span><a href="http://orcid.org/0000-0003-3426-308X" target="_blank">http://orcid.org/0000-0003-3426-308X</a><br><a href="https://github.com/wwolff7" target="_blank">https://github.com/wwolff7</a><br><a href="http://buscatextual.cnpq.br/buscatextual/visualizacv.do?id=K4463141A1" target="_blank">http://buscatextual.cnpq.br/buscatextual/visualizacv.do?id=K4463141A1</a></span></span></div></div>
</div></div></div></div></div></div></div></div></div>