0

I am trying to extract all the levels from a particular NetCDF file in R. I can do this manually by extracting each level as one line of code then combining them as a data frame. But this is very long when I have many files. Is it possible to extract all 43 layers in one file?

I have used this How to extract all levels from a netcdf file using the raster package? and Plotting netcdf file with levels in R as a guidance

In essence the nitrate data from
https://www.nodc.noaa.gov/cgi-bin/OC5/woa18/woa18oxnu.pl has the concentration at 43 different depths. Is it possible to extract all the depths for a particular location?

I can do this for one level. But each level represents a depth. Is it possible to get all levels?

I also do not understand the 3rd warning message: In .getCRSfromGridMap4(atts) : cannot process these parts of the CRS: epsg_code=EPSG:4326

I get a different result (0.5 for Jan level 1) but my colleague gets 1.4 for Jan level 1). Is my error due to the above warning?

#this works
Nit_Jan <- brick("~woa18_all_n01_01.nc", stopIfNotEqualSpaced =    
FALSE, varname = "n_an", level = 1)

#this doesn't
Nit_Jan <- brick("~woa18_all_n01_01.nc", stopIfNotEqualSpaced = 
FALSE, varname = "n_an", level = 1:43)

Warning messages:
1: In if (level <= 0) { :
the condition has length > 1 and only the first element will be used 
2: In if (oldlevel != level) { :
the condition has length > 1 and only the first element will be used
3: In .getCRSfromGridMap4(atts) : cannot process these parts of the   
CRS:epsg_code=EPSG:4326

I would like to plot nitrate by depth

liv
  • 93
  • 10

1 Answers1

2

There is some confusion here because the file has "levels" (a fourth dimension) but the number of levels is one (so no fourth dimension). The code should probably detect that, but for now you have to add lvar=4 to get the desired object.

library(raster)
f <- "woa18_all_n01_01.nc"
b <- brick(f, var="n_oa", lvar=4)
b
#class      : RasterBrick 
#dimensions : 180, 360, 64800, 43  (nrow, ncol, ncell, nlayers)
#resolution : 1, 1  (x, y)
#extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#source     : woa18_all_n01_01.nc 
#names      : X0, X5, X10, X15, X20, X25, X30, X35, X40, X45, X50, X55, X60, X65, X70, ... 
#meters     : 0, 800 (min, max)
#varname    : n_oa 
#level      : 1

And now you can do

pt <- cbind(121.5, 27.5)
e <- extract(b, pt)
e[1:5]
#[1] 10.43725 10.37617 10.23662 13.76292 13.65862

Warning #3

3: In .getCRSfromGridMap4(atts) : cannot process these parts of the CRS:epsg_code=EPSG:4326

can be ignored; but I will fix that in the next version. I think the best thing here is 

crs(b) <- "+init=EPSG:4326"

PS: The development version of raster now behaves better:

f <- "woa18_all_n01_01.nc"
brick(f, var="n_oa")
#class      : RasterBrick 
#dimensions : 180, 360, 64800, 43  (nrow, ncol, ncell, nlayers)
#resolution : 1, 1  (x, y)
#extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#crs        : +init=EPSG:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
#source     : woa18_all_n01_01.nc 
#names      : X0, X5, X10, X15, X20, X25, X30, X35, X40, X45, X50, X55, X60, X65, X70, ... 
#depth (meters): 0, 800 (min, max)
#varname    : n_oa 
#level      : 1
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • I tried to open it as suggested above but I get just one value for 43 levels. The data spans -180, 180, -90, 90 and I want one location stn1 = matrix(c(-150, -30), ncol = 2) stn1 extract(e, stn1) Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘extract’ for signature ‘"list", "matrix"’ – liv May 22 '19 at 21:23
  • Could it be because of the way the data is structured? test <- brick("/Users/Desktop/woa18_all_n01_01.nc", stopIfNotEqualSpaced = FALSE, varname = "n_an") print(test) class : RasterBrick dimensions: 180, 360, 64800, 1 (nrow, ncol, ncell, nlayers) resolution: 1, 1 (x, y) extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 data source : /Users/Desktop/woa18_all_n01_01.nc names : X696.5 z-value: 696.5 varname: n_an level : 1 I thought names meant variable and z-value depth. I feel like im incorrect.. – liv May 22 '19 at 21:30
  • I tried cbind(-150, -30) but the data is still incorrect. Is there an error because of the z-value? When the value is extracted is says X696.5. Which is the same as the z-value. The max depth is 800 (shouldn't this be the z-value?). Sorry for all the confusion! – liv May 22 '19 at 21:39
  • That location is on land in Australia --- all `NA`. The z-value label "X696.5" is not the value (and may not be correct as it does not change --- I will check that – Robert Hijmans May 22 '19 at 22:46
  • Terribly sorry! I actually wanted it off the coast of Australia (Pacific Ocean) side so 160W, 30E (-160, -30). Im guessing that all the 'grids' don't have data. I downloaded the .csv file and there are big gaps. – liv May 22 '19 at 23:02
  • Also how did you get 57 levels when the data only has 43 levels? I.e 43 depth levels from 0 to 800? – liv May 22 '19 at 23:31
  • The file I used (oxygen?) has 57 levels. I do not know why. – Robert Hijmans May 23 '19 at 00:45
  • Thank you @RobertHijmans! Can I just clarify what the z-value and names mean when it says X696.5? I can't seem to find it anywhere online names : X696.5 z-value : 696.5. I also found out why my values were different then a colleagues - the data doesn't exist for that exact coordinate. On ferret, when the coordinate is selected, it provides the coordinates to the closest data point. When we averaged the grid points in ferret, we get the same value I get here using R. Is there a way to know the lat and long the data is extracted from in R? – liv May 23 '19 at 11:12
  • I looked at the file and updated my answer. The "level" in the file represents time, and the value is 696.5 --- I do not know what that represents. – Robert Hijmans May 23 '19 at 17:22
  • @RobertHijmans Is there a function to automatically stack all levels from a ncdf file? I am meaning when you have a time dimension (nlayers) and a z-dimension (levels). – Herman Toothrot Aug 16 '20 at 20:12