1

My goal is to plot nitrate (no3) data on a world map, using the correct longitude and latitude for these data.

There are two netcdf files:
1. with the data
2. with the grid information

Summary info on the data: no3 is an array of length x*y*sigma no3_df is 'x*y obs. of 3 variables' x = integer [180] y = integer [193] sigma = array[53]

I want to look at sigma ('depth') 20. I therefore did the following:

# Load the needed libraries to handle netcdf files
library(ncdf)
library(akima)

# Open data and grid files
file1 <- open.ncdf(file.choose())
grid  <- open.ncdf(file.choose())

# Read relevant variables/parameters from data file1
x <- get.var.ncdf(file1,varid="x")
y <- get.var.ncdf(file1,varid="y")
sigma <- get.var.ncdf(file1,varid="sigma")
no3 <- get.var.ncdf(file1,varid="no3")
sigma_plot <- no3[,,sigma=20]

# Read relevant variables/parameters from grid file
plon <- get.var.ncdf(grid,varid="plon")
plat <- get.var.ncdf(grid,varid="plat")

# Each cell of sigma_plot corresponds to one cell of plon and plat.
A <- array(c(plon,plat,sigma_plot),dim=c(180,193,3))

# Now B is an array containing for each row: (longitude, latitude, value).
B <- apply(A, 3, cbind)

# But it is not a regular grid, so interpolate to a regular grid. akima library
C <- interp(B[,1],B[,2],B[,3], 
            xo=seq(-180,180,1),yo=seq(-90,90,by=1), # tweak here the resolution
            duplicate='mean') # extra y values are duplicates

#########
# PLOTTING
#########

# This one works, but doesn't have a correct longitude and latitude:
filled.contour(x,y,sigma_plot, col=rich.colors(18))

# Try to plot with lon and lat
filled.contour(C, col=rich.colors(30))

Since the filled.contour plot doesn't have correct longitude and latitude, I would like to use ggplot. However, I don't know how to do this...

# And the plotting with ggplot
ggplot(aes(x=plon_datafrm,y=plat_datafrm),data=no3_df) +
  geom_raster() +
  coord_equal() +
  scale_fill_gradient()

This doesn't seem to work. I am net to ggplot so that might be the reason, I would truly appreciate any help.

user2003479
  • 121
  • 4
  • Did you mean `filled.contour(plon,plat,sigma_plot, col=rich.colors(18))`? – plannapus Mar 10 '14 at 13:21
  • No, the filled.contour just uses the normal x[180] and y[193] since it gives an error when I feed it lon and lat. – user2003479 Mar 10 '14 at 13:24
  • Is the result of `dim(sigma_plot)` c(180,193) or c(193,180)? – plannapus Mar 10 '14 at 13:25
  • Might be that you need to transpose `sigma_plot` first: `filled.contour(plon,plat,t(sigma_plot), col=rich.colors(18))`. What error does it raise? – plannapus Mar 10 '14 at 13:25
  • So plon and plat are x*y matrices as described in the script. So no3 is a dataset of x*y*sigma dimension. The 'parameter' grid consists of (a.o.) two variables that provide a relationship between x,y and lon (plon) and x,y and lat (plat). So for example plon contains a x*y matrix with for every x,y combination the relevant longitude. – user2003479 Mar 10 '14 at 13:30
  • dim(sigma_plot) = c(180,193) – user2003479 Mar 10 '14 at 13:31
  • filled.contour(plon,plat,t(sigma_plot), col=rich.colors(18)) Error in filled.contour(plon, plat, t(sigma_plot), col = rich.colors(18)) : increasing 'x' and 'y' values expected – user2003479 Mar 10 '14 at 13:32
  • Could you please clean up your code in your question then, because it is not really posible to understand it (you don't declare x and y, you use no3 before declaring it,...). Also linking to the netcdf files would help considerably. – plannapus Mar 10 '14 at 13:33
  • But I think I would like to use ggplot instead of filled.contour because of its advanced possibilities. – user2003479 Mar 10 '14 at 13:33
  • I understand you prefer ggplot to base plots but it will eventually all come down to having the proper x y z (here longitude, latitude and nitrate) for plotting. Feeding it to filled.contour, ggplot or image+contour afterwards will be the trivial part I think. – plannapus Mar 10 '14 at 13:39
  • Yes it is way clearer, thanks. You cannot really upload it here but you can link to somewhere else. – plannapus Mar 10 '14 at 14:09
  • You didn't need to put my answer in your question, i m working on it. Your file was really not as i thought so it s just plain wrong. – plannapus Mar 10 '14 at 15:08
  • Ok sorry! Thank you ver much for your time. I thought I will summarize the latest full script in my question – user2003479 Mar 10 '14 at 15:12
  • Can you provide a basic link to the file that's not via such a seedy looking site? I can get the data file but not the grid.nc. Is it public data from somewhere? – mdsumner Mar 10 '14 at 20:27
  • I addded Dropbox links, the data are not from somewhere else on the internet. Thank you – user2003479 Mar 11 '14 at 09:52
  • @mdsumner: is this useful to you? Does anyone have any suggestions on my approach with ggplot? – user2003479 Mar 12 '14 at 16:11

1 Answers1

3
library(ncdf)
data <- open.ncdf(file1)
no3 <- get.var.ncdf(data,varid="no3")
sigma_plot <- no3[,,20]
grid <- open.ncdf(file2)
plon <- get.var.ncdf(grid,varid="plon")
plat <- get.var.ncdf(grid,varid="plat")

Contrary to what I previously understood, each cell of sigma_plot corresponds to one cell of plon and plat.

A <- array(c(plon,plat,a),dim=c(180,193,3))
B <- apply(A, 3, cbind)

Now B is an array containing for each row: (longitude, latitude, value). But it is not a regular grid, so you need to interpolate a regular grid. Easiest way would be using interp from package akima:

library(akima)
C <- interp(B[,1],B[,2],B[,3], 
            xo=seq(-180,180,1),yo=seq(-90,90,by=1), #you can tweak here the resolution
            duplicate='mean') #for some reasons some entries are duplicates, i don t know how you want to handle it.

image(C) #for instance, or filled.contour if you prefer
library(maptools)
data(wrld_simpl)
plot(wrld_simpl, add=TRUE, col="white") #To add a simple world map on top
plannapus
  • 18,529
  • 4
  • 72
  • 94
  • Thank you for your answer! It doesn't work however.. I am trying to find a way to share my files. apply(cbind(x,y),1,function(i)plon[i[1],i[2]]) Warning message: In cbind(x, y) : number of rows of result is not a multiple of vector length (arg 1) – user2003479 Mar 10 '14 at 14:26
  • then x and y are not indexes. What are they? – plannapus Mar 10 '14 at 14:26
  • Ah ok, wait a minute i ll update with maybe something more relevant. (In the meantime you can have a look at [this relevant q/a](http://stackoverflow.com/questions/13700254/the-variable-from-a-netcdf-file-comes-out-flipped/13701358#13701358)) – plannapus Mar 10 '14 at 14:27
  • Edited. I'm assuming here that the dimensions in your file2 are actually named x and y. – plannapus Mar 10 '14 at 14:31
  • x and y are the 'coordinates' at which the no3 concentration is available. plon and plat are thus an equivalent of x and y only expressed in longitude-latitude. So I would think I can just provide either plon,plat and the data or x,y and the data. The second worked in filled.contour. So maybe it goes wrong in ggplot (what format is needed to be to make ggplot happy..).. – user2003479 Mar 10 '14 at 14:31
  • Mmm the dimensions of the grid file have x and y, but it returns the following: Warning message: In cbind(x, y) : number of rows of result is not a multiple of vector length (arg 1) – user2003479 Mar 10 '14 at 14:40
  • In my last code, x and y are still the same x and y you had since the beginning so there is no reason it wouldn t work now but did before. – plannapus Mar 10 '14 at 14:43
  • data file http://speedy.sh/vF3Cy/NOINYOC-lr08.micom.hbgcm.0001-08.nc and then two lines would need to be added: grid <- open.ncdf(file.choose()) and file1<-open.ncdf(...) – user2003479 Mar 10 '14 at 14:45
  • What I see is just latitudinal distorsion (plus the fact that filled.contour(x,y,sigma_plot) is not centered on greenwich meridian but around india). Other than those two things the color pattern is exactly similar – plannapus Mar 10 '14 at 15:44
  • Yes this looks good! The map looks like it has correct lon and lat and the same values as with my initial plotting using x,y,sigma_plot. I am not sure though what is happening after the C<- command. Why is it 'not regular' (there is one data point for every lon (-180 to +180) and lat (-90 to +90) isn't that OK?). What do you mean with extra y values? – user2003479 Mar 10 '14 at 15:51
  • Honestly, you have to ask the person who did the ncdf i cannot possibly know that. Look at your plon and plat matrices and you will see that x and y were not strictly a vector of increasing latitude/longitude. They go back and forth something (shipboard analysis?). See in particular the last few rows of plon and plat. – plannapus Mar 10 '14 at 15:52
  • Ok I will ask around here. Thank you for all your help. I am still curious about a ggplot solution for this though; hope that someone else can help with that. – user2003479 Mar 10 '14 at 15:59
  • I cannot find the double data you are referring to @plannapus. The plon at least is constant in every row (as it should be) and goes from -179 to +179 (though not starting at -179 in row 1. It starts at -109, goes to +179 from there and finishes with the last part of the sequence from -179 to -111). I am looking into plat. Can you maybe explain me in words wat is happening behind the C<- command? – user2003479 Mar 12 '14 at 16:27
  • it interpolates lineraly what the value at the nodes of a regular grid (of coordinates x0=seq(-180,180,1) and y0=seq(-90,90,1)) would be based on the values at your actual datapoints. – plannapus Mar 13 '14 at 06:51