6

I have been trying plot the following gridded netcdf file: "air.1999.nc" found at the following website:

http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.html

I have tried the code below based on answers I have found here and elsewhere, but no luck.

library(ncdf);
temp.nc <- open.ncdf("air.1999.nc");
temp <- get.var.ncdf(temp.nc,"air");

temp.nc$dim$lon$vals -> lon
temp.nc$dim$lat$vals -> lat

lat <- rev(lat)
temp <- temp[nrow(temp):1,]

temp[temp==-32767] <- NA
temp <- t(temp)

image(lon,lat,temp)
library(maptools)
data(wrld_simpl)
plot(wrld_simpl, add = TRUE)

This code was from modified from the one found here: The variable from a netcdf file comes out flipped

Does anyone have any ideas or experience with using these type of netcdf files? Thanks

Community
  • 1
  • 1
fjd
  • 65
  • 1
  • 1
  • 4
  • Does this example help ? [rworldmap RJournal paper](http://journal.r-project.org/archive/2011-1/RJournal_2011-1_South.pdf) page 41 : Using `rworldmap` to map netCDF data in tandem with `ncdf` – Andy Dec 11 '13 at 15:11
  • It goes someway to giving an example, but I was rather hoping I would be able to see a working example of the code required to plot the data I mentioned. But thanks for the link, shall explore it further. – fjd Dec 11 '13 at 15:22
  • Having had a quick look what I would do is use the package `raster` which is now better than mine for doing raster and netcdf operations. Try this stackoverflow post and links from it : [http://stackoverflow.com/questions/19330710/why-do-i-get-different-results-in-plotting-a-netcdf-layer-with-imagex-y-z-an](http://stackoverflow.com/questions/19330710/why-do-i-get-different-results-in-plotting-a-netcdf-layer-with-imagex-y-z-an) – Andy Dec 11 '13 at 18:15

1 Answers1

10

In the question you linked the whole part from lat <- rev(lat) to temp <- t(temp) was very specific to that particular OP dataset and have absolutely no universal value.

temp.nc <- open.ncdf("~/Downloads/air.1999.nc")
temp.nc
[1] "file ~/Downloads/air.1999.nc has 4 dimensions:"
[1] "lon   Size: 144"
[1] "lat   Size: 73"
[1] "level   Size: 12"
[1] "time   Size: 365"
[1] "------------------------"
[1] "file ~/Downloads/air.1999.nc has 2 variables:"
[1] "short air[lon,lat,level,time]  Longname:Air temperature Missval:32767"
[1] "short head[level,time]  Longname:Missing Missval:NA"

As you can see from these informations, in your case, missing values are represented by the value 32767 so the following should be your first step:

temp <- get.var.ncdf(temp.nc,"air")
temp[temp=="32767"] <- NA

Additionnaly in your case you have 4 dimensions to your data, not just 2, they are longitude, latitude, level (which I'm assuming represent the height) and time.

temp.nc$dim$lon$vals -> lon
temp.nc$dim$lat$vals -> lat
temp.nc$dim$time$vals -> time
temp.nc$dim$level$vals -> lev

If you have a look at lat you see that the values are in reverse (which image will frown upon) so let's reverse them:

lat <- rev(lat)
temp <- temp[, ncol(temp):1, , ] #lat being our dimension number 2

Then the longitude is expressed from 0 to 360 which is not standard, it should be from -180 to 180 so let's change that:

lon <- lon -180

So now let's plot the data for a level of 1000 (i. e. the first one) and the first date:

temp11 <- temp[ , , 1, 1] #Level is the third dimension and time the fourth.
image(lon,lat,temp11) 

And then let's superimpose a world map:

library(maptools)
data(wrld_simpl)
plot(wrld_simpl,add=TRUE)

enter image description here

plannapus
  • 18,529
  • 4
  • 72
  • 94
  • Thanks for the reply and example of the code that works, much appreciated. – fjd Dec 17 '13 at 09:41
  • 4
    The lon adjustment is incorrect. You can see Africa in the Pacific ocean. It should be ``lon[lon > 180] <- lon[lon > 180] - 360``, I believe. – kennyB Jun 12 '16 at 02:04