[R-sig-Geo] reading into R IPCC data files
Robert Hijmans
r.hijmans at gmail.com
Tue Jan 18 00:04:04 CET 2011
I agree with Andy, although dealing with ncdf files has its own problems
((although you can open most with raster: brick(filename))
So, for the record, reading these ASCII files is not that hard if you have
R. See below for an example that works on the one file I looked at. You
could easily extend this to automate downloading and processing all the
files with functions download.file() and uzip()
Robert
library(raster)
filename <- HADCM3_A1F_DSWF_1980.mea
x <- readLines(filename)
x <- apply(matrix(x),1, function(x) gsub(" ", " ", x))
x <- apply(matrix(x),1, function(x) trim (gsub(" ", " ", x)))
start = which(substr(x,1,5)=="7008 ")+1
end = c(which(substr(x,1,4)=="IPCC")[-1] - 1, length(x))
nc = as.integer(substr(x[2], 8, 11))
nr = as.integer(substr(x[2], 13, 16))
# "Grid is 96 * 73 Month is Jan"
if (length(start) != 12) {
stop('something wrong')
} else {
d = matrix(nrow=nc *nr, ncol=12)
month = NULL
for (m in 1:12) {
d[,m] <- as.numeric(unlist(strsplit(x[start[m]:end[m]], " ")))
mm <- x[start[m]-4]
mm <- trim(substr(mm, nchar(mm)-3, nchar(mm)))
month <- c(month, mm)
}
}
b <- brick(nrow=nr, ncol=nc)
b <- setValues(b, d)
b <- rotate(b)
ln <- paste(month, ";", x[1], ";", x[3], ";", x[4])
layerNames(b) <- ln
plot(b, 1)
plot(b)
b <- writeRaster(b, 'filename.tif')
--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/reading-into-R-IPCC-data-files-tp5925687p5933842.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
More information about the R-sig-Geo
mailing list