<div dir="ltr"><div class="gmail_default" style="font-family:verdana,sans-serif">Senhores, bom dia!</div><div class="gmail_default" style="font-family:verdana,sans-serif"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif">Também tenho interesse no uso de dados dessa natureza, sendo assim, aproveitei o código do Yuri e as dicas do Thiago pra fazer um exercício. Deixo o script que utilizei, caso possa ser útil.</div><div class="gmail_default" style="font-family:verdana,sans-serif"><br></div><div class="gmail_default"><font face="monospace, monospace">### <code r></font></div><div class="gmail_default" style=""><div class="gmail_default"><font face="monospace, monospace">setwd("D:/Temp")</font></div><div class="gmail_default"><font face="monospace, monospace">sapply(c("raster", "sp", "ncdf4"), require, char=T)</font></div><div class="gmail_default"><font face="monospace, monospace">url0 <- "<a href="https://drive.google.com/file/d/0B6tnaf2fmoCcOXhnRDIxWnM5UDA/view">https://drive.google.com/file/d/0B6tnaf2fmoCcOXhnRDIxWnM5UDA/view</a>"</font></div><div class="gmail_default"><font face="monospace, monospace">browseURL(url0) ### 266 Mb</font></div><div class="gmail_default"><font face="monospace, monospace"><br></font></div><div class="gmail_default"><font face="monospace, monospace">dir(patt="\\.nc")</font></div><div class="gmail_default"><font face="monospace, monospace">fn <- "Rs_daily_UT_Brazil_v2_20070101_20131231.nc"</font></div><div class="gmail_default"><font face="monospace, monospace">download.file(url1, fn, mode="wb")</font></div><div class="gmail_default"><font face="monospace, monospace"><br></font></div><div class="gmail_default"><font face="monospace, monospace"># nc <- nc_open(fn); names(nc[['var']]) ### inspeção de conteúdo</font></div><div class="gmail_default"><font face="monospace, monospace"># r <- raster(fn); dim(r); object.size(r) ### importará somente primeiro layer</font></div><div class="gmail_default"><font face="monospace, monospace"><br></font></div><div class="gmail_default"><font face="monospace, monospace">s <- stack(fn, varname="Rs")</font></div><div class="gmail_default"><font face="monospace, monospace">dim(s); print(object.size(s), units="Mb") ### 2557 layers (datas), 30.4 Mb</font></div><div class="gmail_default"><font face="monospace, monospace"><br></font></div><div class="gmail_default"><font face="monospace, monospace">POI <- SpatialPoints(cbind(lon=c(-47.5,-47.5,-44.5,-40), lat=c(-18.5,-19.5,-20.5,-22)))</font></div><div class="gmail_default"><font face="monospace, monospace">proj4string(POI) <- proj4string(s)</font></div><div class="gmail_default"><font face="monospace, monospace"><br></font></div><div class="gmail_default"><font face="monospace, monospace">e <- extract(s, POI, method= "bilinear", df=T) #ou method="simple"</font></div><div class="gmail_default"><font face="monospace, monospace"><br></font></div><div class="gmail_default"><font face="monospace, monospace">dim(e); head(e[,1:5])</font></div><div class="gmail_default"><font face="monospace, monospace"><br></font></div><div class="gmail_default"><font face="monospace, monospace">e2 <- t(e[-1]) ### inverte, retirando a primeira coluna (ID)</font></div><div class="gmail_default"><font face="monospace, monospace">dim(e2) ### [1] 2557 4</font></div><div class="gmail_default"><font face="monospace, monospace"><br></font></div><div class="gmail_default"><font face="monospace, monospace">head(e2, 3)</font></div><div class="gmail_default"><font face="monospace, monospace"># [,1] [,2] [,3] [,4]</font></div><div class="gmail_default"><font face="monospace, monospace"># X2007.01.01 11.67222 12.36980 8.660182 NA</font></div><div class="gmail_default"><font face="monospace, monospace"># X2007.01.02 17.28591 14.08362 9.963414 NA</font></div><div class="gmail_default"><font face="monospace, monospace"># X2007.01.03 15.78652 12.72562 10.278209 NA</font></div><div class="gmail_default"><font face="monospace, monospace"><br></font></div><div class="gmail_default"><font face="monospace, monospace">tail(e2, 3)</font></div><div class="gmail_default"><font face="monospace, monospace"># [,1] [,2] [,3] [,4]</font></div><div class="gmail_default"><font face="monospace, monospace"># X2013.12.29 23.22835 23.38377 19.95641 NA</font></div><div class="gmail_default"><font face="monospace, monospace"># X2013.12.30 18.80989 21.45013 16.14789 NA</font></div><div class="gmail_default"><font face="monospace, monospace"># X2013.12.31 19.01340 17.74833 20.40014 NA</font></div><div class="gmail_default"><font face="monospace, monospace"><br></font></div><div class="gmail_default"><font face="monospace, monospace">### Comparação de métodos de extração (utilizando primeiro layer)</font></div><div class="gmail_default"><font face="monospace, monospace">sp1 <- as(s[[1]], "SpatialPixelsDataFrame")</font></div><div class="gmail_default"><font face="monospace, monospace">over(POI, sp1)[,1]</font></div><div class="gmail_default"><font face="monospace, monospace"><br></font></div><div class="gmail_default"><font face="monospace, monospace">data.frame(POI@coords, </font></div><div class="gmail_default"><font face="monospace, monospace"> over = over(POI, sp1)[,1], </font></div><div class="gmail_default"><font face="monospace, monospace"> bil = extract(s[[1]], POI, method= "bilinear", df=T)[,2],</font></div><div class="gmail_default"><font face="monospace, monospace"> simp = extract(s[[1]], POI, method= "simple", df=T)[,2])</font></div><div class="gmail_default"><font face="monospace, monospace"># lon lat over bil simp</font></div><div class="gmail_default"><font face="monospace, monospace"># 1 -47.5 -18.5 11.586952 11.672215 11.639958</font></div><div class="gmail_default"><font face="monospace, monospace"># 2 -47.5 -19.5 12.564141 12.369798 11.964512</font></div><div class="gmail_default"><font face="monospace, monospace"># 3 -44.5 -20.5 8.624653 8.660182 7.975467</font></div><div class="gmail_default"><font face="monospace, monospace"># 4 -40.0 -22.0 NA NA NaN</font></div><div class="gmail_default"><font face="monospace, monospace"><br></font></div><div class="gmail_default"><font face="monospace, monospace">### Visualizando os dois primeiros pontos</font></div><div class="gmail_default"><font face="monospace, monospace">crop <- crop(s[[1]], extent(-48, -47, -20, -18))</font></div><div class="gmail_default"><font face="monospace, monospace">image(crop, asp=T)</font></div><div class="gmail_default"><font face="monospace, monospace">text(crop, lab=round(getValues(crop), 6), cex=.5); points(POI, pch=20, cex=2)</font></div><div class="gmail_default"><font face="monospace, monospace"><br></font></div><div class="gmail_default"><font face="monospace, monospace">### over pega à esquerda, simple à direita e bilinear interpola.</font></div><div class="gmail_default"><span style="font-family:monospace,monospace">### </code></span><br></div><div class="gmail_default"><span style="font-family:monospace,monospace"><br></span></div><div class="gmail_default"><span style="font-family:monospace,monospace"><img src="cid:ii_ikpdlayd1_152ea04c66469442" width="339" height="508"><br><br></span></div></div><div class="gmail_extra"><br clear="all"><div><div class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><font face="arial, helvetica, sans-serif" style="font-size:small"><div style="font-family:'trebuchet ms',sans-serif;display:inline"></div>================================================<br>Éder Comunello</font><div style="font-size:small"><span style="font-family:arial,helvetica,sans-serif">PhD Student in Agricultural Systems Engineering (USP/Esalq)</span><br></div><div><span style="font-size:small">Brazilian Agricultural Research Corporation (</span><font face="arial, helvetica, sans-serif" style="font-size:small">Embrapa)</font><div style="font-size:small"><font face="arial, helvetica, sans-serif">Dourados, MS, Brazil [</font>22 16.5'S, 54 49.0'W<span style="font-family:arial,helvetica,sans-serif">]</span></div><div><div><br></div><div><br></div></div><div style="font-size:small"><br></div></div></div></div></div></div></div></div></div>
<br><div class="gmail_quote">Em 15 de fevereiro de 2016 15:19, Thiago V. dos Santos <span dir="ltr"><<a href="mailto:thi_veloso@yahoo.com.br" target="_blank">thi_veloso@yahoo.com.br</a>></span> escreveu:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Oi Yury,<br>
<br>
<br>
Pelo que entendi, você quer ter certeza que o código está realmente extraindo as informações de que você precisa.<br>
<br>
Nesse caso a melhor opção seria usar o ncview para extrair a serie temporal da coordenada que você precisa, e então comparar com os valores que o seu código recuperou do arquivo.<br>
<br>
Veja a seguir um código que adaptei do seu para comparar a série temporal extraída pelo seu código com a que eu extrai do ponto adjacente (-47.375 e -18.375) usando o ncview. Para rodar a parte final desse código você vai precisar do arquivo que eu anexei a essa mensagem.<br>
<br>
Antes, algumas observações sobre o seu código:<br>
<br>
1) Você não precisa usar o comando raster, precisa usar apenas o comando stack. A diferença é que o comando raster carrega apenas uma "fatia" do arquivo netcdf, ao passo que o comando stack carrega toda a "pilha". Para ilustrar: você está trabalhando com um arquivo de radiação, cujas medidas vão de 1 de janeiro de 2007 a 31 de dezembro de 2013. Usando o comando raster, você carrega apenas a primeira fatia dos dados, que seria o dia 1 de janeiro/07. Usando o comando stack vc "empilha" todas as fatias de dados - de 1 jan/07 a 31 dec/13 - em uma só variável (que tem o nome de 's' no seu código).<br>
<br>
2) Para definir os pontos de seu interesse, você precisa apenas do cbind... Isso já fornece para o comando "extract" as informações que ele precisa. O comando SpatialPoints serve apenas para converter aquelas coordenadas em pontos georeferenciados. (experimente dar um plot(sp) pra ver o que acontece). Qualquer um dos dois funciona para o comando "extract", mas vc não precisa necessariamente dessa linha no seu código.<br>
<br>
3) As coordenadas desse arquivo netcdf não coincidem exatamente com as coordenadas que você busca. Por exemplo, você busca o ponto (-47.5 e -18.5), mas as coordenadas mais próximas no arquivo são -47.375 e -47.625 na longitude e -18.375 e -18.625 na latitude. Como a coordenada que você busca não "cai" exatamente dentro de uma célula do arquivo, o comando extract vai fazer uma média das 4 células adjacentes e entregar esse valor. Note que última coordenada que você passou cai no oceano, cujos valores no arquivo são inválidos. Por isso a última coluna do dataframe mostra apenas NA's.<br>
<br>
Veja a comparação entre os dados extraídos pelo pacote raster e pelo visualizador ncview:<br>
<br>
--------------------------------------<br>
#carrega pacote<br>
library(raster)<br>
<br>
#carrega arquivo netcdf, especificando a variável desejada (util caso haja mais de uma)<br>
s <- stack("~/Downloads/Rs_daily_UT_Brazil_v2_20070101_20131231.nc", varname="Rs")<br>
<span class=""><br>
#criar lista de pontos a serem extraidos os valores de Rs<br>
</span>longitude <- c(-47.5,-47.5,-44.5,-40)<br>
latitude <- c(-18.5,-19.5,-20.5,-22)<br>
xy <- cbind(longitude,latitude)<br>
<br>
#extração dos valores<br>
e <- extract(s, xy, method= "bilinear", df=T)<br>
<br>
# Melhorar aparencia do data frame<br>
df <- data.frame(t(e)[-1,])<br>
head(df)<br>
<br>
# Abre output de texto do ncview para comparacao<br>
ncv <- read.csv('~/Downloads/ncview.dump', header=F, skip=4, sep='')<br>
head(ncv)<br>
<br>
# Cria um plot simples comparando os valores, apenas para inspecao visual<br>
avalia <- data.frame(df$X1, ncv$V2)<br>
plot(df$X1, ncv$V2)<br>
<br>
<br>
--------------------------------------<br>
<br>
Uma nota que vale a pena mencionar, sobre esse conjunto de dados que você está usando, é que ele já foi criado a partir de dados observados. Por sinal, junto de cada arquivo de dados existe um arquivo complementar com as medidas de erro associadas à espacialização (interpolação) dos dados. Se me lembro bem, essas medidas de erro são: distância da estação mais próxima e erro médio quadrático da medida do ponto em relação a essa estação.<br>
Vale a pena dar uma olhada no paper que descreve esse conjuntos de dados:<br>
<a href="http://onlinelibrary.wiley.com/doi/10.1002/joc.4518/abstract" rel="noreferrer" target="_blank">http://onlinelibrary.wiley.com/doi/10.1002/joc.4518/abstract</a>.<br>
<br>
O primeiro autor, Alexandre, é quem produziu esse conjunto de dados.<br>
Caso você esteja planejando "validar" esses dados, pode valer mais a pena investir tempo na leitura do paper e utilização das medidas de erro. (e eu digo "pode" porque não conheço o objetivo final).<br>
<br>
<br>
Era essa a validação do código que você buscava?<br>
<br>
Hope this helps,<br>
-- Thiago V. dos Santos<br>
<br>
PhD student<br>
Land and Atmospheric Science<br>
University of Minnesota<br>
<div><div class="h5"><br>
<br>
On Monday, February 15, 2016 7:09 AM, Yury Duarte <<a href="mailto:yurynepomuceno@gmail.com">yurynepomuceno@gmail.com</a>> wrote:<br>
<br>
<br>
<br>
Bom dia colegas programadores.<br>
<br>
Estive trabalhando na obtenção de de dados climáticos a partir de um arquivo de extensão .nc, com um gride de 0,25 grau para todo o território nacional.<br>
Minha necessidade é de extrair essas informações climáticas, dia a dia, em formato numérico para alguns pontos específicos do mapa, que irão coincidir com os pontos das estações meteorológicas de onde coletei os dados observados.<br>
Após algumas tentativas e com algumas ajudas dos colegas, cheguei no código que segue ao final do texto junto com o link para baixar o arquivo analisado.<br>
Gostaria de validar esse código, pois não estou familiarizado com esse tipo de arquivo e com o pacote do R que usei para abri-lo. Como não conheço o conteúdo do arquivo, também não sei se houve uma limitação na visualização da saída do comando View(e) pelo R.<br>
Desde já, agradeço pela ajuda de todos!<br>
<br>
Link para download do arquivo:<br>
<br>
<a href="https://drive.google.com/file/d/0B6tnaf2fmoCcOXhnRDIxWnM5UDA/view" rel="noreferrer" target="_blank">https://drive.google.com/file/d/0B6tnaf2fmoCcOXhnRDIxWnM5UDA/view</a><br>
<br>
<br>
Código do R:<br>
<br>
#seleção de diretorio<br>
setwd("C:/Users/user/Desktop/Novo_BD")<br>
#seleção de pacote<br>
library(raster)<br>
#seleção do arquivo de interesse<br>
r=raster("Rs_daily_UT_Brazil_v2_19800101_19891231.nc")<br>
#seleção da variavel de interesse dentro do arquivo<br>
s=stack("Rs_daily_UT_Brazil_v2_19800101_19891231.nc", varname="Rs")<br>
#criar lista de pontos a serem extraidos os valores de Rs<br>
longitude=c(-47.5,-47.5,-44.5,-40)<br>
latitude=c(-18.5,-19.5,-20.5,-22)<br>
xy=cbind(longitude,latitude)<br>
sp=SpatialPoints(xy)<br>
#extração dos valores<br>
e=extract(s, sp, method= "bilinear", df=T) #ou method="simple"<br>
#View(e)<br>
# Consertar data frame<br>
df=data.frame(t(e)[,])<br>
View(df)<br>
<br>
<br>
<br>
<br>
Yury Duarte<br>
Engenheiro Agrônomo - ESALQ/USP<br>
<br>
</div></div>_______________________________________________<br>
R-br mailing list<br>
<a href="mailto:R-br@listas.c3sl.ufpr.br">R-br@listas.c3sl.ufpr.br</a><br>
<a href="https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br" rel="noreferrer" target="_blank">https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br</a><br>
Leia o guia de postagem (<a href="http://www.leg.ufpr.br/r-br-guia" rel="noreferrer" target="_blank">http://www.leg.ufpr.br/r-br-guia</a>) e forne� c�igo m�imo reproduz�el.<br>_______________________________________________<br>
R-br mailing list<br>
<a href="mailto:R-br@listas.c3sl.ufpr.br">R-br@listas.c3sl.ufpr.br</a><br>
<a href="https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br" rel="noreferrer" target="_blank">https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br</a><br>
Leia o guia de postagem (<a href="http://www.leg.ufpr.br/r-br-guia" rel="noreferrer" target="_blank">http://www.leg.ufpr.br/r-br-guia</a>) e forneça código mínimo reproduzível.<br></blockquote></div><br></div></div>