The file Description of VERA file format.pdf indicates that the VERA test case files have 46 lines of header,
and although the column names are not given (in a form that R will understand), they are given in the
information file, where it can also be seen that some of the columns are not used. The code below demonstrates
how to read a VERA file into R, remove the unused columns, add a header, assign one column (precipitation) to
a separate R object, read in the longitude and latitude coordinate information, and make a nice plot of the
file. For this tutorial, the VERA files VERA_8km_2007062007_01_bog.dat, VERA_8km_2007062008_01_bog.dat and
VERA_8km_2007062009_01_bog.dat will be used. [path] below indicates the path to these data sets on your own
computer, and should be changed accordingly. It should also end with a "/", and be wrapped in quotes
(e.g., [path] might be "/home/user/").
In the code below, the map function is from the package maps (Original S code by Richard A. Becker and
Allan R. Wilks. R version by Ray Brownrigg. Enhancements by Thomas P. Minka R package
version 2.3-6), and poly.image and image.plot are from the package fields (Nychka et al. 2014). The first
allows maps to be added to plots, and the latter allows for plotting projections and adding the color scale,
respectively. Both of these packages are automatically installed and loaded when you install and load
SpatialVx ("R>" indicates that the command is typed from the R prompt, and should not be typed/copied,
in cases where a command goes to more than one line, the "R>" is not given, and the following lines are
indented).
R> path <- [path]
R> Y <- read.table( paste( path, "VERA_8km_2007062007_01_bog.dat", sep = "" ), skip = 46 )
R> Y <- Y[, -c(3, 4, 10, 11, 13:17) ]
R> colnames( Y ) <- c("x-coordinate (km, distance from origin)",
"y-coordinate (km, distance from origin)",
"precipitation (mm/h)",
"10m u-wind (m/s)",
"10m v-wind (m/s)",
"2m potential temperature (deg C)",
"2m equivalent potential temperature (deg C)",
"msl - pressure (hPa)",
"moisture flux divergence (kg / (kg * 10 # (-3)), post processing)",
"moisture flux divergence (kg / (kg * 10 # (-3)), post processing)",
"mixing ratio (kg / kg * 10#(-3), post processing)" )
R> summary( Y )
R> loc <- read.table( paste(path, VERA_8km_coordinates_lam_phi.txt", sep = ""))
R> mloc <- list( x = matrix( loc[,1], 209, 193, byrow = TRUE ),
y = matrix( loc[,2], 209, 193, byrow = TRUE ) )
R> lr <- apply( loc, 2, range, finite = TRUE )
R> Precip <- matrix( Y[, 3], 209, 193, byrow = TRUE )
R> Precip[ Precip < 0 ] <- 0
R> map( xlim = lr[,1], ylim = lr[, 2], type = "n" )
R> poly.image( mloc$x, mloc$y, Precip, col = c("gray", tim.colors(64) ), add = TRUE )
R> map( add = TRUE )
R> image.plot( Precip, legend.only = TRUE, horizontal = TRUE, col = c("gray", tim.colors(64) ) )
The COSMO2 (00 UTC) files valid at the times of the above VERA analysis are co2_8km_2007062007+07.dat,
co2_8km_2007062008+08.dat and co2_8km_2007062009+09.dat, and those for CMC-GEMH (06 UTC) are
cmh_8km_2007062007+01.dat, cmh_8km_2007062008+02.dat and cmh_8km_2007062009+03.dat. Using information from the
model file description (Description of model file format.pdf), the following code will read them into R and
grab only the precipitation fields from them. Additionally, values of 9999 have been added to the model fields
in order to be blown up to the VERA model domain. Here, they are set to NA for simplicity.
R> M1a <- read.table(paste([path], "co2_8km_2007062007+07.dat", sep = ""), skip = 46)
R> M1b <- read.table(paste([path], "co2_8km_2007062008+08.dat", sep = ""), skip = 46)
R> M1c <- read.table(paste([path], "co2_8km_2007062009+09.dat", sep = ""), skip = 46)
R> M2a <- read.table(paste([path], "cmh 8km 2007062007+01.dat", sep = ""), skip = 46)
R> M2b <- read.table(paste([path], "cmh 8km 2007062008+02.dat", sep = ""), skip = 46)
R> M2c <- read.table(paste([path], "cmh 8km 2007062009+03.dat", sep = ""), skip = 46)
R> M1a <- M1a[, -c(3, 4, 10, 11, 14:17) ]
R> M1b <- M1b[, -c(3, 4, 10, 11, 14:17) ]
R> M1c <- M1c[, -c(3, 4, 10, 11, 14:17) ]
R> M2a <- M2a[, -c(3, 4, 10, 11, 14:17) ]
R> M2b <- M2b[, -c(3, 4, 10, 11, 14:17) ]
R> M2c <- M2c[, -c(3, 4, 10, 11, 14:17) ]
R> colnames( M1a) <- colnames( M1b ) <- colnames( M1c ) <-
R> colnames( M2a ) <- colnames( M2b ) <- colnames( M2c ) <-
c( "x-coordinate (km, distance from origin)",
"y-coordinate (km, distance from origin)",
"precipitation (mm/h)",
"10m u-wind (m/s)",
"10m v-wind (m/s)",
"2m potential temperature (deg C)",
"2m equivalent potential temperature (deg C)",
"msl - pressure (hPa) reduced with model formulae",
"msl - pressure (hPa) reduced with std pressure reduction formulae",
"mixing ratio ( kg / kg * 10#(-3), post processing)",
"moisture flux divergence (kg / (kg * s#(-1)* 10#(-4)), post processing)" )
R> COSMOprecip1 <- matrix( M1a[,3], 209, 193, byrow = TRUE )
R> COSMOprecip2 <- matrix( M1b[,3], 209, 193, byrow = TRUE )
R> COSMOprecip3 <- matrix( M1c[,3], 209, 193, byrow = TRUE )
R> GEMHprecip1 <- matrix( M2a[,3], 209, 193, byrow = TRUE )
R> GEMHprecip2 <- matrix( M2b[,3], 209, 193, byrow = TRUE )
R> GEMHprecip3 <- matrix( M2c[,3], 209, 193, byrow = TRUE )
R> COSMOprecip1[ COSMOprecip1 == 9999 ] <- NA
R> COSMOprecip2[ COSMOprecip2 == 9999 ] <- NA
R> COSMOprecip3[ COSMOprecip3 == 9999 ] <- NA
R> GEMHprecip1[ GEMHprecip1 == 9999 ] <- NA
R> GEMHprecip2[ GEMHprecip2 == 9999 ] <- NA
R> GEMHprecip3[ GEMHprecip3 == 9999 ] <- NA
The above can now be plotted in similar fashion as the VERA data. Note that wherever an NA exists, the plot will just be white in color, so they may at first glance not look very similar to the VERA.