[R] Writing a NetCDF file in R

Steve Murray smurray444 at hotmail.com
Tue Aug 4 17:28:51 CEST 2009


Dear all,

I am attempting to convert 10 NetCDF files into a single NetCDF file, due to the data input requirements of a model I hope to use. I am using the ncdf package, version 1.6. The data are global-scale water values, on a monthly basis for 10 years (ie. 120 months of data in total; at present the data are separated by year, with 12 months of data in each file - mrunoff_1986 through to mrunoff_1995). For each month there are 720 longitude x 360 latitude values, each with a corresponding runoff value (although some of these may be NAs).

Here is my code so far:


# READ IN NetCDF FILES FROM DISK

library(ncdf)

year <- 1986:1995
file_list <- paste("mrunoff-",year,".nc", sep="")
file_list

start <- 1986

for (i in file_list) {
        assign(paste("netcdf_",start,"_temp", sep=""),open.ncdf(i))
        start = start+1
          }

# Start converting 10 files into 1

latitude <- get.var.ncdf(netcdf_1986_temp, "Lat")
longitude <- get.var.ncdf(netcdf_1986_temp, "Lon")
month <- get.var.ncdf(netcdf_1986_temp, "Mon")

year <- 1986:1995

mrunoff_1986 <- get.var.ncdf(netcdf_1986_temp, "mrunoff")
mrunoff_1987 <- get.var.ncdf(netcdf_1987_temp, "mrunoff")
mrunoff_1988 <- get.var.ncdf(netcdf_1988_temp, "mrunoff")
mrunoff_1989 <- get.var.ncdf(netcdf_1989_temp, "mrunoff")
mrunoff_1990 <- get.var.ncdf(netcdf_1990_temp, "mrunoff")
mrunoff_1991 <- get.var.ncdf(netcdf_1991_temp, "mrunoff")
mrunoff_1992 <- get.var.ncdf(netcdf_1992_temp, "mrunoff")
mrunoff_1993 <- get.var.ncdf(netcdf_1993_temp, "mrunoff")
mrunoff_1994 <- get.var.ncdf(netcdf_1994_temp, "mrunoff")
mrunoff_1995 <- get.var.ncdf(netcdf_1995_temp, "mrunoff")


# Define variable dimensions

dimx <- dim.def.ncdf("Lon", "deg E", as.double(longitude))
dimy <- dim.def.ncdf("Lat", "deg N", as.double(latitude))
month <- dim.def.ncdf("Mon", "Months: Jan 86=1, Dec 95=120)", 1:120)
year <- dim.def.ncdf("Year", "year", year)

# Assign data: extract mrunoff from each of the 10 files and put into one place

mrunoff_data <- dim.def.ncdf("mrunoff", "mm/month", c(mrunoff_1986, mrunoff_1987, mrunoff_1988, mrunoff_1989, mrunoff_1990, mrunoff_1991, mrunoff_1992, mrunoff_1993, mrunoff_1994, mrunoff_1995))

# Define runoff variable

mrunoff_dims <- var.def.ncdf("mrunoff_out", "mm/month", list(dimx, dimy, month), -9999.0, "Global Monthly Runoff for 1986-1995", "double")

# Create file

mrunoff_file <- create.ncdf("mrunoff.nc", mrunoff_dims)

# Put mrunoff data into the file

put.var.ncdf(mrunoff_file, mrunoff_dims, mrunoff_data)

# Write to disk

# close.ncdf(mrunoff_file)




However, when I run the code, I get the following error message:

> put.var.ncdf(mrunoff_file, mrunoff_dims, mrunoff_data)
Error in put.var.ncdf(mrunoff_file, mrunoff_dims, mrunoff_data) : 
  put.var.ncdf: error: you asked to write 31104000 values, but the passed data array only has 8 entries!


I can understand where the 31104000 comes from ((720*360)*12)*10, but am confused as to why only 8 values are being passed to put.var.ncdf. I therefore tried doing a couple of tests to shed some light on this:

> length(dimx); length(dimy); length(mrunoff_dims); length(mrunoff_data)
[1] 8
[1] 8
[1] 9
[1] 8


> str(dimx); str(dimy)
List of 8
 $ name         : chr "Lon"
 $ units        : chr "deg E"
 $ vals         : num [1:720] 0.25 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75 ...
 $ len          : int 720
 $ id           : num -1
 $ unlim        : logi FALSE
 $ dimvarid     : num -1
 $ create_dimvar: logi TRUE
 - attr(*, "class")= chr "dim.ncdf"

List of 8
 $ name         : chr "Lat"
 $ units        : chr "deg N"
 $ vals         : num [1:360] -89.8 -89.2 -88.8 -88.2 -87.8 ...
 $ len          : int 360
 $ id           : num -1
 $ unlim        : logi FALSE
 $ dimvarid     : num -1
 $ create_dimvar: logi TRUE
 - attr(*, "class")= chr "dim.ncdf"


> str(mrunoff_dims); str(mrunoff_data)
List of 9
 $ name    : chr "mrunoff_out"
 $ units   : chr "mm/month"
 $ missval : num -9999
 $ longname: chr "Global Monthly Runoff for 1986-1995"
 $ id      : num -1
 $ prec    : chr "double"
 $ dim     :List of 3
  ..$ :List of 8
  .. ..$ name         : chr "Lon"
  .. ..$ units        : chr "deg E"
  .. ..$ vals         : num [1:720] 0.25 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75 ...
  .. ..$ len          : int 720
  .. ..$ id           : num -1
  .. ..$ unlim        : logi FALSE
  .. ..$ dimvarid     : num -1
  .. ..$ create_dimvar: logi TRUE
  .. ..- attr(*, "class")= chr "dim.ncdf"
  ..$ :List of 8
  .. ..$ name         : chr "Lat"
  .. ..$ units        : chr "deg N"
  .. ..$ vals         : num [1:360] -89.8 -89.2 -88.8 -88.2 -87.8 ...
  .. ..$ len          : int 360
  .. ..$ id           : num -1
  .. ..$ unlim        : logi FALSE
  .. ..$ dimvarid     : num -1
  .. ..$ create_dimvar: logi TRUE
  .. ..- attr(*, "class")= chr "dim.ncdf"
  ..$ :List of 8
  .. ..$ name         : chr "Mon"
  .. ..$ units        : chr "Months since Jan 1986 (Jan 86 =1)"
  .. ..$ vals         : int [1:120] 1 2 3 4 5 6 7 8 9 10 ...
  .. ..$ len          : int 120
  .. ..$ id           : num -1
  .. ..$ unlim        : logi FALSE
  .. ..$ dimvarid     : num -1
  .. ..$ create_dimvar: logi TRUE
  .. ..- attr(*, "class")= chr "dim.ncdf"
 $ ndims   : int 3
 $ unlim   : logi FALSE
 - attr(*, "class")= chr "var.ncdf"

List of 8
 $ name         : chr "mrunoff"
 $ units        : chr "mm/month"
 $ vals         : num [1:31104000] NA NA NA NA NA NA NA NA NA NA ...
 $ len          : int 31104000
 $ id           : num -1
 $ unlim        : logi FALSE
 $ dimvarid     : num -1
 $ create_dimvar: logi TRUE
 - attr(*, "class")= chr "dim.ncdf"



As you can see from this, the length of the objects and the $len values of the objects differ - the latter being the intended values - plus all the data (e.g. 1:31104000) appear to be present in the str of the object. I'm struggling to see therefore why only 8 'bits' of data are being passed on.

Any help would be much appreciated.

Thanks for your suggestions,

Steve


_________________________________________________________________

[[elided Hotmail spam]]




More information about the R-help mailing list