Senhores, bom dia!

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.

### <code r>

setwd("D:/Temp")

sapply(c("raster", "sp", "ncdf4"), require, char=T)

url0 <- "https://drive.google.com/file/d/0B6tnaf2fmoCcOXhnRDIxWnM5UDA/view"

browseURL(url0) ### 266 Mb

dir(patt="\\.nc")

fn <- "Rs_daily_UT_Brazil_v2_20070101_20131231.nc"

download.file(url1, fn, mode="wb")

# nc <- nc_open(fn); names(nc[['var']]) ### inspeção de conteúdo

# r <- raster(fn); dim(r); object.size(r) ### importará somente primeiro layer

s <- stack(fn, varname="Rs")

dim(s); print(object.size(s), units="Mb") ### 2557 layers (datas), 30.4 Mb

POI <- SpatialPoints(cbind(lon=c(-47.5,-47.5,-44.5,-40), lat=c(-18.5,-19.5,-20.5,-22)))

proj4string(POI) <- proj4string(s)

e <- extract(s, POI, method= "bilinear", df=T) #ou method="simple"

dim(e); head(e[,1:5])

e2 <- t(e[-1]) ### inverte, retirando a primeira coluna (ID)

dim(e2) ### [1] 2557 4

head(e2, 3)

# [,1] [,2] [,3] [,4]

# X2007.01.01 11.67222 12.36980 8.660182 NA

# X2007.01.02 17.28591 14.08362 9.963414 NA

# X2007.01.03 15.78652 12.72562 10.278209 NA

tail(e2, 3)

# [,1] [,2] [,3] [,4]

# X2013.12.29 23.22835 23.38377 19.95641 NA

# X2013.12.30 18.80989 21.45013 16.14789 NA

# X2013.12.31 19.01340 17.74833 20.40014 NA

### Comparação de métodos de extração (utilizando primeiro layer)

sp1 <- as(s[[1]], "SpatialPixelsDataFrame")

over(POI, sp1)[,1]

data.frame(POI@coords,

over = over(POI, sp1)[,1],

bil = extract(s[[1]], POI, method= "bilinear", df=T)[,2],

simp = extract(s[[1]], POI, method= "simple", df=T)[,2])

# lon lat over bil simp

# 1 -47.5 -18.5 11.586952 11.672215 11.639958

# 2 -47.5 -19.5 12.564141 12.369798 11.964512

# 3 -44.5 -20.5 8.624653 8.660182 7.975467

# 4 -40.0 -22.0 NA NA NaN

### Visualizando os dois primeiros pontos

crop <- crop(s[[1]], extent(-48, -47, -20, -18))

image(crop, asp=T)

text(crop, lab=round(getValues(crop), 6), cex=.5); points(POI, pch=20, cex=2)

### over pega à esquerda, simple à direita e bilinear interpola.

### </code>

================================================
Éder Comunello

PhD Student in Agricultural Systems Engineering (USP/Esalq)

Brazilian Agricultural Research Corporation (Embrapa)

Dourados, MS, Brazil [22 16.5'S, 54 49.0'W]

Em 15 de fevereiro de 2016 15:19, Thiago V. dos Santos <thi_veloso@yahoo.com.br> escreveu:

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 <yurynepomuceno@gmail.com> 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
R-br@listas.c3sl.ufpr.br
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.
_______________________________________________
R-br mailing list
R-br@listas.c3sl.ufpr.br
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.