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.