[R-br] Obtenção de Valores em Pontos de Arquivo Raster
Éder Comunello
comunello.eder em gmail.com
Terça Fevereiro 16 10:24:49 BRST 2016
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 em 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 em 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 em 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 em 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 em 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.
>
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20160216/6bacdf5c/attachment.html>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: plot01.png
Tipo: image/png
Tamanho: 12928 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20160216/6bacdf5c/attachment.png>
Mais detalhes sobre a lista de discussão R-br