Netcdf
NetCDF files
NetCDF is a sophisticated self describing file format that is much used in climatology. A formal description and tools and more can be found here. Sophistication has a price and NetCDF files are difficult to use at first. You can use the Intel array viewer to visually inspect these files.
Below are two examples of R scripts that you can use. You will need to make many change to make it work for your cdf file though! The first file reads "station type" data. The second script reads "grid type" data (and computes averages). You can also use Perl, Python and other scripting languages that can link the netCDF C++ or Fortan libraries. Or you can incorporate it into your Java or Delphi applications.
##Start of script##
#removes all object from memory
rm(list=ls(all=TRUE))
#load the ncdf library
library(ncdf)
# now you should read up on this library a bit:
?ncdf
#Set the working directory
setwd("E:/PROJECTS/Climate/")
#opens the cdf data file
datafile<-"data.cdf"
nc<-open.ncdf(datafile)
# now you can inspect the contents of nc (by typing "nc")
nc
#This is just an example for the specific files that was used here.
#load the data for IWMO, latitude,longitude and elevation to user define variables
dLat <-get.var.ncdf(nc, "lat")
dLon <-get.var.ncdf(nc, "lon")
dElev <-get.var.ncdf(nc, "elev")
dIWMO <-get.var.ncdf(nc, "IWMO")
#close the CDF file
close.ncdf(nc)
#binds the data together column-wise and rediret the output to the variable 'stations'
stations <- cbind(dIWMO, dLat, dLon, dElev)
#create an text output file
fname <- "stations.txt"
#write the content of the variable 'stations' to the external output file
write.table(stations, fname, append = FALSE, quote = FALSE, sep = " ",
eol = "\n", na = "-9999", dec = ".", row.names = FALSE,
col.names = FALSE)
##end of script
#######################################################
# R script to export NetCDF climate data from e.g. PMIP
# to Arc-ASCII grids
#
# Robert Hijmans, March 2005
#
# This scipt assumes that the input data are of the form:
# climatevariable(lon, lat, time)
# and that time is in units of months
# there should thus be x * 12 time elements with x = (1..n)
# data can be stored in multiple files (nfiles)
# but the filename conventions are somewhat hard-coded,
# assuming 10 years per file, starting at startyear
#
#######################################################
# first clean up (remove all objects)
#######################################################
rm(list=ls(all=TRUE))
library(ncdf)
#######################################################
# user must change the following
#######################################################
modelname <- "CCSM"
timeperiod <- "21k"
startyear <- 200
nfiles <- 10
workpath <- "e:/pmip/work/"
varlist <- list("tasmin", "tasmax", "pr")
#######################################################
# end user changes section #
#######################################################
setwd(workpath)
for (vv in varlist) {
varname <- paste(vv, "_A_MO_pmip2_", timeperiod, "_oa_", modelname, sep="")
outname <- paste(vv, "_", timeperiod, "_", modelname, sep="")
datapath <- "../pmip2_dbext/"
datafile <- paste(varname,"_",sep="")
iters <- 0
for (f in 1:nfiles) {
if (nfiles > 1) {
#######################################################
# somewhat hard-coded filename conventions start here #
#######################################################
x <- f - 1
start <- startyear + ((f-1) * 10)
fin <- start + 9
if (start < 10) {start <- paste("00", start, sep="")}
else if (start < 100) {start <- paste("0", start, sep="")}
if (fin < 10) {fin <- paste("00", fin, sep="")}
else if (fin < 100) {fin <- paste("0", fin, sep="")}
dataf <- paste(datafile, "0", start, "-", "0", fin, ".nc", sep="")
vname <- paste(varname, "0", start, "-", "0", fin, ".nc", sep="")
#######################################################
# end hard-coded filenames #
#######################################################
} else {
dataf <- datafile
vname <- varname
}
input <- paste(datapath, dataf, sep="")
nc <- open.ncdf(input)
datain <- get.var.ncdf(nc)
#######################################################
# if you want to inspect lat and lon values
#######################################################
# lon <- nc$var$lon
# ndims <- lon$ndims
# dlon <- get.var.ncdf(nc, lon$dim[[ndims]]$name)
# filename <- paste(varname,"_lon.txt", sep="")
# write.table(dlon, filename)
# lat <- nc$var$lat
# ndims <- lat$ndims
# dlat <- get.var.ncdf(nc, lat$dim[[ndims]]$name)
# filename <- paste(varname,"_lat.txt", sep="")
# write.table(dlat, filename)
#######################################################
close.ncdf(nc)
#######################################################
# create dataout (and other arrays) only the first time!
#######################################################
if (f == 1) {
nlon <- length(datain[,1,1])
nlat <- length(datain[1,,1])
ntime <- length(datain[1,1,])
mnt <- array(0,12)
dataout <- array(0, c(nlon, nlat, 12, nfiles))
}
#######################################################
# BEGIN for each cell (lat, long)
#######################################################
for (i in 1:nlon) {
for (j in 1:nlat) {
t <- 0
for (m in 1:12) {mnt[m] <- 0}
#######################################################
# sum values for each month across "years" (ntime)
#######################################################
while (t < ntime) {
for (m in 1:12) {
mnt[m] <- mnt[m] + datain[i,j,t+m]
}
t <- t + 12
}
#######################################################
# now average the data over over "years" (ntime) and store
# note the "* 12" because ntime = years * months
#######################################################
for (m in 1:12) {
dataout[i,j,m,f] <- 12 * mnt[m] / ntime
}
}
}
#######################################################
# END for each cell (lat, long)
#######################################################
# test output
# for (m in 1:12) {
# filename <- paste(vname, "_", m, ".txt", sep="")
# write.table(dataout[,,m,f], filename)
# }
iters <- iters + 1
}
#######################################################
# we have read the data for all files
# now average the data over the files
#######################################################
# check to avoid that this part continues if a .cdf file was missing.
if (iters < nfiles) {mywarn <- "not good"} else {
mywarn <- "ok"
alldataout <- array(0, c(nlon, nlat, 12))
for (i in 1:nlon) {
for (j in 1:nlat) {
for (m in 1:12) {
for (f in 1:nfiles) {
alldataout[i, j, m] <- alldataout[i, j, m] + dataout[i, j, m, f]
}
alldataout[i, j, m] <- alldataout[i, j, m] / nfiles
}
}
}
#######################################################
# change lat to go from N to S instead of S to N
# change lon so that it goes from -180 to 180 instead of from 0 to 360
#######################################################
finaldata <- array(0, c(nlat, nlon, 12))
halflon <- nlon / 2
for (i in 1:nlat) {
ii <- (nlat - i + 1)
for (j in 1:nlon) {
if (j > halflon) {jj <- j - halflon} else {jj <- j + halflon}
for (m in 1:12) {
finaldata[i, j, m] <- alldataout[jj, ii, m]
}
}
}
#######################################################
# Write an arc-ascii file for each month
#######################################################
for (m in 1:12) {
filename <- paste(outname, "_", m, ".asc", sep="")
thefile <- file(filename, "w") # open an output file connection
cat("ncols ",nlon, "\n", file = thefile, sep = " ")
cat("nrows ",nlat, "\n", file = thefile, sep = " ")
cat("xllcorner -180", "\n", file = thefile, sep = " ")
cat("yllcorner -90", "\n", file = thefile, sep = " ")
cat("cellsize ", 360/nlon, "\n", file = thefile, sep = " ")
cat("NODATA_value -9999", "\n", file = thefile, sep = " ")
close(thefile)
write.table(finaldata[,,m], filename, append = TRUE, quote = FALSE, sep = " ",
eol = "\n", na = "-9999", dec = ".", row.names = FALSE,
col.names = FALSE)
}
}
}
mywarn
iters
nfiles
nc$file
Last modified June 5, 2008 12:12 am