4

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 413x229 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.

4

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.

  • This was exactly what I was looking for, thanks a lot. I really could not get my head around the proj problem. Just a couple short questions descending from the above. You'll see that using levelplot(r) or gplot(r) + geom_tile(aes(fill=value)) + coord_equal() followed by map(add=TRUE, col="blue") leads to some problems, namely, the map is somehow shifted. What's the reason for this and/or a possible fix? Simple plot(r) does not show this behavior, as you noted. – AF7 Nov 5 '14 at 10:06
  • Answer: Of course when using map() the overlay is relative to the plot(). For ggplot, one can do: world <-data.frame(map(plot=FALSE)[c("x","y")]) and then gplot(r*86400) + geom_tile(aes(fill=value)) + geom_path(aes(x,y), data=world) + coord_equal(xlim=c(xmin(r), xmax(r)), ylim=c(ymin(r), ymax(r))) – AF7 Nov 5 '14 at 11:51
  • Hello, I have similar problem, but sadly could not solve it given tips provided stackoverflow.com/questions/43785505/… – MIH May 5 '17 at 10:50
0

Sometimes, like a variable extract from a HDF (hdf5) file, you don`t have the regular cells spacing. The lat and lon extract from the file, for example, are irregular spacing. So you are going to have problems to make a correct raster or image. Therefore, you can interpolate the data. I use, for this, the function interp from package akima. This is the way how I solve the problem when there is a matriz or data set with irregular values of lat and lon. I think this is util for use to ggplot2 too.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.