I have a netCDF file (in this LINK) which contains monthly average values for 4 parameters at 11 different depth levels (0-200m). I would like to extract the monthly value for temperature "votemper" at the maximum valid depth (i.e. when the water is 20 m deep there are only NA vlaues for greater depths) however, the last valid depth varies spatially based on the water bathymetry. I have got as far as turning the netCDF into an array in R but am having trouble with how to extract just the values from the maximum valid depth value.
1 Answer
You can use the raster
package and its functionalities regarding netCDF files to address your problem. Make sure the ncdf
package is available.
# Required packages
library(raster)
library(ncdf)
library(RColorBrewer)
You can manually retrieve some information about the netCDF file you are dealing with using open.ncdf()
. This step is important for subsequent loops along available depth levels and time steps included in the file.
# Information about ncdf file, e.g. variables, dimensions, no. of time steps
nc <- open.ncdf("oceandata.nc")
nc
[1] "file oceandata.nc has 4 dimensions:"
[1] "time Size: 3"
[1] "depth Size: 12"
[1] "lat Size: 23"
[1] "lon Size: 15"
[1] "------------------------"
[1] "file oceandata.nc has 4 variables:"
[1] "short vosaline[lon,lat,depth,time] Longname:Sea Water Salinity Missval:-32768"
[1] "short vomecrty[lon,lat,depth,time] Longname:Northward Current Velocity in the Water Column Missval:-32768"
[1] "short votemper[lon,lat,depth,time] Longname:Sea Water Potential Temperature Missval:-32768"
[1] "short vozocrtx[lon,lat,depth,time] Longname:Eastward Current Velocity in the Water Column Missval:-32768"
Now that you know how much time steps (3) and depth levels (12) are included, you can start importing the data into R. However, I am still wondering whether you can directly import daily RasterBricks consisting of 12 layers, i.e. one for each depth level (resulting in 3 RasterBricks in this particular case), rather than RasterBricks consisting of all available RasterLayers per depth level (status quo, resulting in 12 RasterBricks comprising three layers (3 time steps) each).
# Import ncdf file as 'Raster*' object by looping through all available
# depth levels (Size: 12), returning a list containing one RasterBrick with
# all available time steps per depth level
lst_by_depth <- lapply(1:12, function(i) {
brick("oceandata.nc", varname = "votemper", lvar = "depth", level = i)
})
rst_by_depth <- brick(lst_by_depth)
Now that you have the RasterBricks for each depth level, you can start reordering your data by manually concatenating all available depth levels that correspond to the same day, resulting in a list of daily RasterBrick objects.
# Loop through no. of time steps (Size: 3), returning a list containing one
# RasterBrick with all available depth levels per day
lst_by_time <- lapply(1:3, function(i) {
indices <- seq(i, nlayers(rst_by_depth), 3)
rst_by_depth[[indices]]
})
Next, you need to define a function that you can later on pass to calc
for identifying the last valid (i.e. non-NA) temperature value on a pixel basis. If no NA values exist, it will return the deepest possible value. If no non-NA values exist (I suppose this should correspond to land, then), it returns NA.
# Define a function that, for each cell,
# - returns NA in case all depth level values are NA
# - returns the value of the last available depth level in case no NA occur
# - else returns the value of the last non-NA depth level
deepestValid <- function(x) {
na <- is.na(x)
if (all(na)) {
return(NA)
} else if (all(!na)) {
return(x[length(x)])
} else {
first_na <- which(na)[1]
last_valid <- first_na - 1
return(x[last_valid])
}
}
Finally, apply this function to each element, i.e. to each daily RasterBrick, in the previously generated list. This results in a list of three daily RasterLayers holding the temperature value of the deepest possible point per pixel.
# Apply custom function to each daily RasterBrick, returning daily RasterLayers
# with deepest available value
deepest_valid <- lapply(lst_by_time, function(i) {
calc(i, fun = deepestValid)
})
col <- colorRampPalette(rev(brewer.pal(9, "RdBu")))
spplot(deepest_valid[[1]], col.regions = col(100), cuts = 50)
However, I have to confess that the approach seems rather long-winded to me, and I am still wondering whether there is a more convenient solution to your problem.
-
Hi flowla, i think that is a very well documented answer - thanks! I tried a different approach by converting to a data frame, removing NA values and then selecting the maximum depth level but the process of converting from netCDF to large matrix didn't feel very efficient.– xyzblueOct 2, 2014 at 10:27