NetCDF – Plotting NetCDF File Using Lat and Lon Contained in Variables

netcdfr

I am trying to plot this NetCDF file (and similar) as a 2d contour map.
The file contains a single layer of one variable, "pr".
The file has 413×229 X-Y pixel coordinates, and, stored inside the "xlon" and "xlat" variables, the exact coordinates of each grid point.
I can easily plot using X-Y coordinates:

filename <- "file.nc"
varname <- "pr"
library(raster)
variable <- raster(filename, varname=varname)
plot(variable)

Hpwever, the plot comes out "squeezed", and moreover I want to plot using lat-lon coordinates. How can I do that?

This is somehow similar, I think, to this question:
https://gis.stackexchange.com/questions/48518/how-to-re-project-ease-equal-area-scalable-earth-grid-with-a-25-km-cylindrica?newreg=0eac93c7cbf04b6386d4e62e64e747ef

However, the solution given there is not suited to my case because:

  • I would like to avoid using the ncdf4 package, I'd prefer to use raster only as this would allow me to expand the procedure to bigger files
  • I get several errors in the convertGrid() function, such as:
    coordinates(xp)=~x+y
    Error in .subset(x, j) : only 0's may be mixed with negative subscripts

I should reverse the projection, I think. All I know of the projection is contained in the file ncdump:

projection = "LAMCON"
grid_size_in_meters = 12000.
latitude_of_projection_origin = 39.
longitude_of_projection_origin = 14.
standard_parallel = 35., 51. 
grid_factor = 0.684241343018562

The file does not have a CRS specification.
I know I should be able to do so using the sp, raster and rasterVis packages but so far I've not succeeded.

Trying to plot as a data.frame does not work either because the grid is not constant.

Best Answer

I guess that my answer from your other StackOverflow question did not lead you in the right direction?

Here is a more detailed answer that may nudge you in the right direction.

First, we need to know the projection of your data and the extent to be able to project to the Long/Lat grid correctly. Unfortunately, we do not have the PROJ4 CRS or the extent, so I try my best to estimate those.

# First, estimate the extent.
# We start with the lat and long, then project it to the Lambert Conformal projection
library(raster)
inputfile <- "pr_averaged_Med_CORDEX_ATM.1980-2008.nc"

# Grab the lat and lon from the data
lat <- raster(inputfile, varname="xlat")
lon <- raster(inputfile, varname="xlon")

# Convert to points and match the lat and lons
plat <- rasterToPoints(lat)
plon <- rasterToPoints(lon)
lonlat <- cbind(plon[,3], plat[,3])

# Specify the lonlat as spatial points with projection as long/lat
lonlat <- SpatialPoints(lonlat, proj4string = CRS("+proj=longlat +datum=WGS84"))

# Need the rgdal package to project it to the original coordinate system
library("rgdal")

# My best guess at the proj4 string from the information given
mycrs <- CRS("+proj=lcc +lat_1=35 +lat_2=51 +lat_0=39 +lon_0=14 +k=0.684241 +units=m +datum=WGS84 +no_defs")
plonlat <- spTransform(lonlat, CRSobj = mycrs)
# Take a look
plonlat
extent(plonlat)

# Yay! Now we can properly set the coordinate information for the raster
pr <- raster(inputfile, varname="pr")
# Fix the projection and extent
projection(pr) <- mycrs
extent(pr) <- extent(plonlat)
# Take a look
pr
plot(pr)

# Project to long lat grid
r <- projectRaster(pr, crs=CRS("+proj=longlat +datum=WGS84"))
# Take a look
r
plot(r)
# Add contours
contour(r, add=TRUE)

# Add country lines
library("maps")
map(add=TRUE, col="blue")

I think that looks correct. The country borders near the ocean seem to match up with the colours of the plot. I think you can take it from there. You can make prettier plots with rasterVis or ggplot.

Related Question