如何将多个波段从 netcdf 写入二进制文件？

Question

我有 netcdf 文件 我打开它并读取了一个变量：

 K=open.ncdf("C:\\hiba_history.nc")
Smonthly= get.var.ncdf(nc=K,varid="evap",verbose=TRUE)

[1] "vobjtodimname: is a character type varid.  This file has 9 dims"
[1] "vobjtodimname: no cases found, returning FALSE"
[1] "get.var.ncdf: isdimvar: FALSE"
[1] "vobjtovarid: entering with varid=evap"
[1] "Variable named evap found in file with varid= 10"
[1] "vobjtovarid: returning with varid deduced from name; varid= 10"
[1] "get.var.ncdf: ending up using varid= 10"
[1] "ndims: 3"
[1] "get.var.ncdf: varsize:"
[1] 34 30 12
[1] "get.var.ncdf: start:"
[1] 1 1 1
[1] "get.var.ncdf: count:"
[1] 34 30 12
[1] "get.var.ncdf: totvarsize: 12240"
[1] "Getting var of type 3  (1=short, 2=int, 3=float, 4=double, 5=char, 6=byte)"
[1] "get.var.ncdf: C call returned 0"
[1] "count.nodegen: 34    Length of data: 12240" "count.nodegen: 30    Length of data: 12240"
[3] "count.nodegen: 12    Length of data: 12240"
[1] "get.var.ncdf: final dims of returned array:"
[1] 34 30 12
[1] "varid: 10"

如您所见，这个变量有 30 像素和 34 行和 12 个波段（月） 我只想写 12 的总和，所以我终于得到了一个文件来计算所有 12 个月的总和（每年）

      apply(Smonthly, c(1,2), sum) -> Sannual 
  to.write = file(paste("C:\\annual.bin",sep=""),"wb")

   writeBin(as.double(Sannual),to.write,size=4)

当我通过另一个程序打开文件时，我发现地图（文件）是颠倒的

Answer 1

...一旦您有了每月的总和，您就可以致电writeBin 将它们保存到文件中：

a = array(runif(10*10*10), dim = c(10,10,10))
a_sum = apply(a, c(1,2), sum)
# Write stuff
writeBin(as.numeric(a_sum), "/tmp/test.bin", size = 4)
# read stuff back in
test = readBin("/tmp/test.bin", numeric(), 10*10, size = 4)
# ...succes???
> head(data.frame(test, as.numeric(a_sum)))
      test as.numeric.a_sum.
1 5.581374          5.581374
2 5.974429          5.974429
3 4.854637          4.854637
4 5.040194          5.040193
5 3.709209          3.709210
6 6.119048          6.119048
>     all.equal(test, as.numeric(a_sum))
[1] "Mean relative difference: 2.313248e-08"
>     all.equal(test, as.numeric(a_sum), tolerance = 1e-7)
[1] TRUE

注意：为什么你需要设置tolerance 来让all.equal 返回TRUE 我留给读者作为练习。

Answer 2

就 12 个月的总和而言，方法如下：

library(ncdf)
K <- open.ncdf("math.nc")
Smonthly <- get.var.ncdf(nc=K,varid="evap")
apply(Smonthly, c(1,2), sum) -> Sannual 
# Since the months are represented by dimension 3, you apply sum on dimensions 1 and 2

根据请求，这里尝试使用raster：

library(raster)
library(ncdf)
Smonthly <- raster("math.nc", varname="evap", band=12)
Sannual <- calc(Smonthly, sum)