[R-br] Criar arquivo de saida a partir de dados NetCDF

[Éder Comunello]

Olá, novamente!

Encontrei um erro no código anterior e que está presente também na primeira
solução que sugeri. Ao usar {raster} a definição de extent() estava
incorreta. Os valores x e y correspondem ao centro do pixel e para o valor
de extent() ficar correto deve-se considerar a resolução da imagem. É pouco
perceptível na pequena escala, mas pode causar problemas ao aumentar a
escala. Segue o código corrigido...

### <code r>
### Dados
require(raster)
data(volcano); volcano <- 10*volcano
image(volcano)
x <- 100*(1:nrow(volcano))
y <- 100*(1:ncol(volcano))

########## {ncdf4}
###########################################################
# browseURL("https://www.getdatajoy.com/learn/Read_and_Write_NetCDF_from_R")
library(ncdf4)
nc1_dim1  <- ncdim_def('EW', 'm', as.double(x))
nc1_dim2  <- ncdim_def('SN', 'm', as.double(y))
nc1_varz  <- ncvar_def('Elevation', 'm', list(nc1_dim1, nc1_dim2), -1,
longname = 'Data about a Volcano')
nc1_out   <- nc_create('volcano.nc', nc1_varz, force_v4 = TRUE)
ncvar_put(nc1_out, nc1_varz, volcano)
nc_close(nc1_out)

print(nc_open("volcano.nc"))
plot(raster("volcano.nc"))

########## {RNetCDF}
#########################################################
# browseURL("https://journal.r-project.org/archive/2013-2/michna-woods.pdf")
require(RNetCDF)
nc2_out <- create.nc("volcano2.nc")
dim.def.nc(nc2_out, "EW", length(x))
dim.def.nc(nc2_out, "SN", length(y))
var.def.nc(nc2_out, "EW", "NC_DOUBLE", "EW")
var.def.nc(nc2_out, "SN", "NC_DOUBLE", "SN")
var.def.nc(nc2_out, "Elevation", "NC_DOUBLE", c("EW", "SN"))
att.put.nc(nc2_out, "Elevation", "_FillValue", "NC_DOUBLE", -999.9)
var.put.nc(nc2_out, "EW", as.double(x))
var.put.nc(nc2_out, "SN", as.double(y))
var.put.nc(nc2_out, "Elevation", volcano)
close.nc(nc2_out)

print(nc_open("volcano2.nc"))    ### {ncdf4}
print.nc(open.nc("volcano2.nc")) ### {RNetCDF}
var.get.nc(nc2, "EW")
var.get.nc(nc2, "Elevation")
plot(raster("volcano2.nc"))

########## {raster}
##########################################################
require(raster)
tmp <- raster(volcano)
r <- t(flip(tmp, 1))
plot(r)

resx <- round(mean(x[-1]-x[-length(x)]), 6) # 100
resy <- round(mean(y[-1]-y[-length(y)]), 6) # 100

extent(r) <- c(range(x), range(y)) + c(-resx, +resx, -resy, +resy)/2
plot(r)
writeFormats()
writeRaster(r, "volcano3.nc", "CDF", over=T)

print(nc_open("volcano3.nc")) ### {ncdf4}
plot(raster("volcano3.nc"))

r1 <- raster("volcano.nc")
r2 <- raster("volcano2.nc")
r3 <- raster("volcano3.nc")

# extents
cbind(sapply(paste0("r", 1:3), function(x) as.vector(extent(get(x)))))
#        r1   r2   r3
# [1,]   50   50   50
# [2,] 8750 8750 8750
# [3,]   50   50   50
# [4,] 6150 6150 6150

{
windows()
par(mfrow=c(2,2))
sapply(paste0("r", 1:3), function(x) plot(get(x), main=x))
par(mfrow=c(1,1))
}

### </code>

​
================================================
Éder Comunello
Agronomist (UEM), MSc in Environ. Sciences (UEM)
DSc in Agricultural Systems Engineering (USP/Esalq)
Brazilian Agricultural Research Corporation (Embrapa)
Dourados, MS, Brazil |<O>|
================================================
GEO, -22.2752, -54.8182, 408m
UTC-04:00 / DST: UTC-03:00
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20160219/9adbe170/attachment.html>