Plot NetCDF variable-grid data file using ggplot2:

问题:

I need to plot some data from this NetCDF file (1.1MB).

The file contains a 413x229 grid (94577 points). Every point has a precipitation value which I need to print in the correct LAT-LON. The grid is not necessarily constant, so that we have two additional 413x229 variables (xlat and xlon) containing the lat-lon values for each grid point.

Thanks to the hints from these questions...

• How to plot contours on a map with ggplot2 when data is on an irregular grid?
• R - Plotting netcdf climate data
• Plotting netcdf in R with correct grid
• How to plot latitude/longitude points over netcdf map
• Plotting netcdf in R with correct grid
• how to display a projected map on an R::lattice::layerplot?

...plotting the data on an arbitrary grid using ggplot is fairly easy:

library(ncdf4)
library(ggplot2)
library(reshape)
ncfile <- nc_open(inputfile.nc)
pr <- ncvar_get(ncfile, "pr")
pr <- pr*86400
mpr <- melt(pr)
ggplot(aes(x=X1, y=X2, fill=value), data=mpr) + geom_raster() + coord_equal()

This will produce a plot with a grid which is not lat-lon of course. However, when trying to plot the same data with the correct grid:

xlon <- ncvar_get(ncfile, "xlon")
xlat <- ncvar_get(ncfile, "xlat")
df <- data.frame(as.vector(lat), as.vector(lon), as.vector(pr))
ggplot(aes(x=lat, y=lon, fill=pr), data=df) + geom_raster() + coord_equal()

ggplot will return a "vector is too large" error and plot nothing.
pr, lat and lon are 413x229 arrays.

So, questions:

1. What am I doing wrong?
2. How can I easily customize fixed non-continuous contour levels in ggplot?
3. How can I overlay a political map of the plotted area? (Alps) I've not looked yet into this as it is of secondary importance. It may be answered in one of the questions listed above, so feel free to ignore it.

EDIT: Answers suggest I should use the raster and rasterVis packages. However, I'm not sure on how to do that. The file does not provide any more info. I know it's a Lambert Conformal projection with :

latitude_of_projection_origin = 39.
longitude_of_projection_origin = 14.
standard_parallel = 35., 51. 
grid_factor = 0.684241343018562

You can see this by "ncdump -h".

However I don't understand how to plot this using levelplot() on the correct lat-lon grid. It really looks troublesome to me. I know how to do this in GrADS quite simply, but GrADS is very limiting and I'd prefer to avoid using it honestly.

回答1:

I remember having a similar problem as you.

Here is what I learned:

As far as I know, ggplot's geom_raster will only plot on a regularly spaced grid, and if you are pulling the unprojected latitudes and longitudes from the raster, you may not have a regularly spaced grid to plot. If you are set on using ggplot to plot your data and you only have the latitudes and longitudes, you will have to find the projection information and project it to its native, regularly spaced grid prior to plotting. ggplot uses projections provided by the mapproject package and does not use PROJ4 projection libraries, which I find to be more powerful and universally supported across other spatial packages in R.

If you are working with a regular grid, which I assume you are, I suggest using the raster package, which can directly read netcdf4 files. For a simple plot, try:

library(raster)
r <- raster("your_nc_file_path", varname = "pr")
plot(r)
contour(r, add = TRUE)

If you have a properly made netcdf file, you can also extract the projection information as well with proj4string(r) or crs(r) or projection(r). You may be able to translate this to work with ggplot's projection system.

Also, if you really, really want to plot this in ggplot, try using the rasterVis package, which conveniently prepares rasters to be used with ggplot.

library(rasterVis)
# Notice `g`plot (from rasterVis) and not `gg`plot
gplot(r) + 
  geom_tile(aes(fill = value))

Hope that helps.