[R-sig-Geo] Problems with ncdf package

[Luke Knowles]

Hi,

I'm trying to make a netcdf file with three variables in. Two  
variables have the same dimensions, [96,73,5,12], and the other has  
different dimensions [96,73,1,12].

I have set verbose to true;

[1] "create.ncdf: input was a list of vars"
[1] "Calling R_nc_create for file  /Users/luke/Ancil_Files/Nc_Files/ 
land_fraction_SDGVM.nc"
[1] "back from R_nc_create for file  /Users/luke/Ancil_Files/Nc_Files/ 
land_fraction_SDGVM.nc"
[1] "create.ncdf: about to create the dims"
[1] "var.add.ncdf: entering with indefine= TRUE"
[1] "var.add.ncdf: ncid of file to add to= 13    filename= /Users/ 
luke/Ancil_Files/Nc_Files/land_fraction_SDGVM.nc     writable= TRUE"
[1] "var.add.ncdf: varname to add= field1393"
[1] "var.add.ncdf: creating 4 dims for var field1393"
[1] "var.add.ncdf: working on dim > longitude < (number 1 ) for var  
field1393"
[1] "create.ncdf: creating dim longitude"
[1] "dim.create.ncdf: entering for ncid= 13"
[1] "dim.create.ncdf: entering for dim longitude"
[1] "dim.create.ncdf: about to call R_nc_def_dim for dim longitude"
[1] "dim.create.ncdf: making dimvar for dim longitude"
[1] "dim.create.ncdf: about to call R_nc_def_var_double for dimvar  
longitude"
[1] "dim.create.ncdf: about to call R_nc_put_vara_double dimvals for  
dimvar longitude"
[1] "dim.create.ncd: exiting for ncid= 13   dim= longitude"
[1] "var.add.ncdf: working on dim > latitude < (number 2 ) for var  
field1393"
[1] "create.ncdf: creating dim latitude"
[1] "dim.create.ncdf: entering for ncid= 13"
[1] "dim.create.ncdf: entering for dim latitude"
[1] "dim.create.ncdf: about to call R_nc_def_dim for dim latitude"
[1] "dim.create.ncdf: making dimvar for dim latitude"
[1] "dim.create.ncdf: about to call R_nc_def_var_double for dimvar  
latitude"
[1] "dim.create.ncdf: about to call R_nc_put_vara_double dimvals for  
dimvar latitude"
[1] "dim.create.ncd: exiting for ncid= 13   dim= latitude"
[1] "var.add.ncdf: working on dim > height < (number 3 ) for var  
field1393"
[1] "create.ncdf: creating dim height"
[1] "dim.create.ncdf: entering for ncid= 13"
[1] "dim.create.ncdf: entering for dim height"
[1] "dim.create.ncdf: about to call R_nc_def_dim for dim height"
[1] "dim.create.ncdf: making dimvar for dim height"
[1] "dim.create.ncdf: about to call R_nc_def_var_double for dimvar  
height"
[1] "dim.create.ncdf: about to call R_nc_put_vara_double dimvals for  
dimvar height"
[1] "dim.create.ncd: exiting for ncid= 13   dim= height"
[1] "var.add.ncdf: working on dim > time < (number 4 ) for var  
field1393"
[1] "create.ncdf: creating dim time"
[1] "dim.create.ncdf: entering for ncid= 13"
[1] "dim.create.ncdf: entering for dim time"
[1] "dim.create.ncdf: about to call R_nc_def_dim for dim time"
[1] "dim.create.ncdf: making dimvar for dim time"
[1] "dim.create.ncdf: about to call R_nc_def_var_double for dimvar time"
[1] "dim.create.ncdf: about to call R_nc_put_vara_double dimvals for  
dimvar time"
[1] "dim.create.ncd: exiting for ncid= 13   dim= time"
[1] "create.ncdf: creating single precision var field1393"
[1] "create.ncdf: C call returned value 0"
[1] "var.add.ncdf: entering with indefine= TRUE"
[1] "var.add.ncdf: ncid of file to add to= 13    filename= /Users/ 
luke/Ancil_Files/Nc_Files/land_fraction_SDGVM.nc     writable= TRUE"
[1] "var.add.ncdf: varname to add= field1392"
[1] "var.add.ncdf: creating 4 dims for var field1392"
[1] "var.add.ncdf: working on dim > longitude < (number 1 ) for var  
field1392"
[1] "var.add.ncdf: working on dim > latitude < (number 2 ) for var  
field1392"
[1] "var.add.ncdf: working on dim > height < (number 3 ) for var  
field1392"
[1] "var.add.ncdf: working on dim > time < (number 4 ) for var  
field1392"
[1] "create.ncdf: creating single precision var field1392"
[1] "create.ncdf: C call returned value 0"
[1] "var.add.ncdf: entering with indefine= TRUE"
[1] "var.add.ncdf: ncid of file to add to= 13    filename= /Users/ 
luke/Ancil_Files/Nc_Files/land_fraction_SDGVM.nc     writable= TRUE"
[1] "var.add.ncdf: varname to add= field1384"
[1] "var.add.ncdf: creating 4 dims for var field1384"
[1] "var.add.ncdf: working on dim > longitude < (number 1 ) for var  
field1384"
[1] "var.add.ncdf: working on dim > latitude < (number 2 ) for var  
field1384"
[1] "var.add.ncdf: working on dim > height < (number 3 ) for var  
field1384"
Error in var.add.ncdf(nc, vars[[ivar]], verbose = verbose, indefine =  
TRUE) :
   This is not allowed.


Here is the code that i have written - if there is anything obvious  
that i have done wrong please let me know - maybe i can't do it like  
this? Do i need to close the ncdf file, reopen it and append the  
different sized array???

library(ncdf)
library(fields)
library(abind)
library(maps)
# HEADER
SDGVM_TYPE <- c 
("Bare_Soil","C3_grass","C4_grass","Evergreen_Broadleaf","Evergreen_Need 
leaf","Deciduous_BroadLeaf","Deciduous_NeedleLeaf","Crops")
SDGVM_TYPES <- c("Bare Soil","C3 grass","C4 grass","Evergreen  
Broadleaf","Evergreen Needleaf","Deciduous Broad Leaf","Deciduous  
Needle Leaf","Crops")
UM_TYPES <- c("Broadleaf","Needleaf","C3 Grasses","C4  
grasses","Shrub","Urban","Water","Bare Soil","Ice")
month <- c 
("January","February","March","April","May","June","July","August","Sept 
ember","October","November","December")
# BODY
land_fraction <- open.ncdf("/Users/luke/Ancil_Files/ 
qrparm_land_frac.nc")
# field1392(LAI) (96,73,5,12)
# field1393(Canopy Height) (96,73,5,12)
SDGVM <- open.ncdf("/Users/luke/Ancil_Files/SDGVM_Nc/VegPresentR8.nc")
# field1391(Fractions Of Surface Types) (96,73,9,1)
IGBP_Func <- open.ncdf("/Users/luke/Ancil_Files/IGBP_Veg_Func.nc")
# field1384(Canopy Conductance) (96,73,1,12)
IGBP_CanCon <- open.ncdf("/Users/luke/Ancil_Files/Nc_Files/ 
Veg_Func_Monthly.nc")
# lai (96,73,8,12)
# maxvegetfrac (96,73,8,12)

pseudo <- get.var.ncdf(land_fraction,"pseudo")
height <- get.var.ncdf(SDGVM,"veg")
lat <- get.var.ncdf(land_fraction,"latitude")
lat_sdgvm <- get.var.ncdf(SDGVM,"lat")
lon <- get.var.ncdf(land_fraction,"longitude")
lon_sdgvm <- get.var.ncdf(SDGVM,"lon")
tme <- get.var.ncdf(land_fraction,"t")
tme_sdgvm <- get.var.ncdf(SDGVM,"time")
#
numlevs<-length(pseudo)
numlevs_sdgvm<-length(height)
numlats<-length(lat)
numlats_sdgvm<-length(lat_sdgvm)
numlong<-length(lon)
numlong_sdgvm<-length(lon_sdgvm)
numtime<-length(tme)
numtime_sdgvm<-length(tme_sdgvm)
#
axis_x_sdgvm<-seq(from=-180, to=180, by=45)
axis_x <- seq(from =0, to=360, by=45)
axis_y_sdgvm<-seq(from=-90, to=90, by=30)
axis_y<-seq(from=-90, to=90, by=30)
#
land_frac_1 <- get.var.ncdf(land_fraction,"field1391",start=c 
(1,1,1,1),count=c(48,73,9,1))
land_frac_2 <- get.var.ncdf(land_fraction,"field1391",start=c 
(48,1,1,1),count=c(48,73,9,1))
land_frac <- abind(land_frac_2,land_frac_1,along=1)
print("Line 48")
assign("land_mask",land_frac)

igbp_urban <- do.call(abind,c(rep(list(land_mask[,,6]),12),along=3))
igbp_water <- land_mask[,,7]
igbp_soil <- land_mask[,,8]
igbp_ice <- do.call(abind,c(rep(list(land_mask[,,9]),12),along=3))
print("Line 55")
igbp_water <- do.call(abind,c(rep(list(igbp_water),12),along=3))
#
sdgvm_lai <- get.var.ncdf(SDGVM,"lai",start=c(1,1,1,1),count=c 
(96,73,8,12))
#
sdgvm_frac <- get.var.ncdf(SDGVM,"maxvegetfrac",start=c 
(1,1,1,1),count=c(96,73,8,12))
#
print("Line 62")
sdgvm_lai_5npft <- abind((sdgvm_lai[,,4,] + sdgvm_lai[,,6,]), 
(sdgvm_lai[,,5,] + sdgvm_lai[,,7,]), (sdgvm_lai[,,2,] + sdgvm_lai[,, 
8,]) ,(sdgvm_lai[,,3,]) , (sdgvm_lai[,,8,]* 0.0),along=2.5)
#
assign("LAI_LF",sdgvm_lai_5npft)
print("Line 66")
sdgvm_broadleaf <- sdgvm_frac[,,4,] + sdgvm_frac[,,6,]
sdgvm_needleaved <- sdgvm_frac[,,5,] + sdgvm_frac[,,7,]
sdgvm_c3 <- sdgvm_frac[,,2,] + sdgvm_frac[,,8,]
sdgvm_c4 <- sdgvm_frac[,,3,]
sdgvm_soil <- sdgvm_frac[,,1,]
sdgvm_crops <- sdgvm_frac[,,8,]
sdgvm_null <- sdgvm_frac[,,8,] * 0.0
#
print("Line 75")
soil_only <- sdgvm_soil[,,] - igbp_ice
#
sdgvm_soil_1 <- ifelse(sdgvm_soil[,1:12,] >= 0.0, 0.0)
sdgvm_soil_2 <- ifelse(sdgvm_soil[,61:73,] >= 0.0, 0.0)
#
print("Line 81")
sdgvm_frac_lf <- abind 
(sdgvm_broadleaf,sdgvm_needleaved,sdgvm_c3,sdgvm_c4,sdgvm_null,igbp_urba 
n,igbp_water,sdgvm_soil,igbp_ice,along=2.5)
#
# Fractional types of average canopy heights c(Broadleaf,Needleleaf,  
C3, C4, Shrubs)
#
canopy_height_average <- c(28,22,1.5,1.3,2.1)
#
CHT_LF <- abind((sdgvm_frac[,,1,]*0.0) , (sdgvm_frac[,,2,] 
*canopy_height_average[3]) , (sdgvm_frac[,,3,]*canopy_height_average 
[4]) , (sdgvm_frac[,,4,]*canopy_height_average[1]) , (sdgvm_frac[,,5,] 
*canopy_height_average[2]) , (sdgvm_frac[,,6,]*canopy_height_average 
[1]) , (sdgvm_frac[,,7,]*canopy_height_average[2]) , (sdgvm_frac[,,8,] 
*canopy_height_average[3]),along=2.5)
#
CHT_5NPFT <- abind((CHT_LF[,,4,] + CHT_LF[,,6,]),(CHT_LF[,,5,] +  
CHT_LF[,,7,]),(CHT_LF[,,2,] + CHT_LF[,,8,]),(CHT_LF[,,3,]),(CHT_LF[,, 
1,]*0.0),along=2.5)
#
igbp_cancon_1 <- get.var.ncdf(IGBP_CanCon,"field1384",start=c 
(1,1,1,1),count=c(50,73,1,12))
igbp_cancon_2 <- get.var.ncdf(IGBP_CanCon,"field1384",start=c 
(50,1,1,1),count=c(46,73,1,12))
#
igbp_cancon <- get.var.ncdf(IGBP_CanCon,"field1384",start=c 
(1,1,1,1),count=c(96,73,1,12))
#
CANCON_IGBP <- abind(igbp_cancon_2,igbp_cancon_1,along=1)
#
lon <- get.var.ncdf(IGBP_CanCon,"longitude")
tme <- get.var.ncdf(IGBP_CanCon,"t")
#
sewd <- "/Users/luke/Desktop/Dropbox/R_Output/CAS_NET/SDGVM/Can_Hgt/"
setwd(sewd)
#
filename <- "/Users/luke/Ancil_Files/Nc_Files/land_fraction_SDGVM.nc"
#
x <- dim.def.ncdf("longitude","Degrees East",seq(-176.25,180.00,3.75))
x_1 <- dim.def.ncdf("longitude","Degrees East",seq(-180.00,176.25,3.75))
y <- dim.def.ncdf("latitude","Degrees North",seq(-90,90,2.5))
z <- dim.def.ncdf("height","Model Level",seq(1,5,1))
z_cancon <- dim.def.ncdf("height","Model Level",1)
w <- dim.def.ncdf("time","days since 0000-01-01 00:00:00",14.00)
w_12 <- dim.def.ncdf("time","days since 0000-01-01 00:00:00",seq 
(14,344,30))
# missing value
mv <- 2.00000004008175E+20
#
field1392 <- var.def.ncdf("field1392","m2_m2",list 
(x,y,z,w_12),mv,longname="LEAF AREA INDEX OF PLANT FUNC TYPES")
field1393 <- var.def.ncdf("field1393","m",list 
(x,y,z,w_12),mv,longname="CANOPY HEIGHT OF PLANT FUNC TYPES")
field1384 <- var.def.ncdf("field1384","",list 
(x,y,z_cancon,w_12),mv,longname="CANOPY CONDUCTANCE AFTER TIMESTEP")
#
nc_file <- create.ncdf(filename,list 
(field1393,field1392,field1384),verbose=TRUE)
#
put.var.ncdf(nc_file,field1384,CANCON_IGBP,verbose="TRUE")
put.var.ncdf(nc_file,field1393,CHT_5NPFT,verbose="TRUE")
put.var.ncdf(nc_file,field1392,LAI_LF,verbose="TRUE")
#
close.ncdf(nc_file)

Thanks

Luke

PhD Student
Centre for Atmospheric Science
Department of Chemistry
University of Cambridge