Oi Yury,
Pelo que entendi, você quer ter certeza que o código está realmente extraindo
as informações de que você precisa.
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.
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.
Antes, algumas observações sobre o seu código:
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).
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.
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.
Veja a comparação entre os dados extraídos pelo pacote raster e pelo
visualizador ncview:
--------------------------------------
#carrega pacote
library(raster)
#carrega arquivo netcdf, especificando a variável desejada (util caso haja mais
de uma)
s <- stack("~/Downloads/Rs_daily_UT_Brazil_v2_20070101_20131231.nc",
varname="Rs")
#criar lista de pontos a serem extraidos os valores de Rs
longitude <- c(-47.5,-47.5,-44.5,-40)
latitude <- c(-18.5,-19.5,-20.5,-22)
xy <- cbind(longitude,latitude)
#extração dos valores
e <- extract(s, xy, method= "bilinear", df=T)
# Melhorar aparencia do data frame
df <- data.frame(t(e)[-1,])
head(df)
# Abre output de texto do ncview para comparacao
ncv <- read.csv('~/Downloads/ncview.dump', header=F, skip=4, sep='')
head(ncv)
# Cria um plot simples comparando os valores, apenas para inspecao visual
avalia <- data.frame(df$X1, ncv$V2)
plot(df$X1, ncv$V2)
--------------------------------------
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.
Vale a pena dar uma olhada no paper que descreve esse conjuntos de dados:
http://onlinelibrary.wiley.com/doi/10.1002/joc.4518/abstract.
O primeiro autor, Alexandre, é quem produziu esse conjunto de dados.
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).
Era essa a validação do código que você buscava?
Hope this helps,
-- Thiago V. dos Santos
PhD student
Land and Atmospheric Science
University of Minnesota
On Monday, February 15, 2016 7:09 AM, Yury Duarte <[email protected]>
wrote:
Bom dia colegas programadores.
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.
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.
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.
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.
Desde já, agradeço pela ajuda de todos!
Link para download do arquivo:
https://drive.google.com/file/d/0B6tnaf2fmoCcOXhnRDIxWnM5UDA/view
Código do R:
#seleção de diretorio
setwd("C:/Users/user/Desktop/Novo_BD")
#seleção de pacote
library(raster)
#seleção do arquivo de interesse
r=raster("Rs_daily_UT_Brazil_v2_19800101_19891231.nc")
#seleção da variavel de interesse dentro do arquivo
s=stack("Rs_daily_UT_Brazil_v2_19800101_19891231.nc", varname="Rs")
#criar lista de pontos a serem extraidos os valores de Rs
longitude=c(-47.5,-47.5,-44.5,-40)
latitude=c(-18.5,-19.5,-20.5,-22)
xy=cbind(longitude,latitude)
sp=SpatialPoints(xy)
#extração dos valores
e=extract(s, sp, method= "bilinear", df=T) #ou method="simple"
#View(e)
# Consertar data frame
df=data.frame(t(e)[,])
View(df)
Yury Duarte
Engenheiro Agrônomo - ESALQ/USP
_______________________________________________
R-br mailing list
[email protected]
https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forne� c�igo m�imo
reproduz�el.
ncview.dump
Description: Binary data
_______________________________________________ R-br mailing list [email protected] https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forne�a c�digo m�nimo reproduz�vel.
