plotting netcdf data for selected region

Question

I am finding it difficult to work with the netcdf data I have. I went through the examples here (Plotting netcdf file with levels in R) but I think im missing something.

I am trying to plot mixed layer depth for region south of 35 degrees. The data can be found here (the last file at the bottom of the page): http://www.ifremer.fr/cerweb/deboyer/mld/Surface_Mixed_Layer_Depth.php

the file has 7 variables each containing the lat, lon, time (12 months), and the value which is the mixed layer depth.

So far I have:

    MLD <- "mld_DReqDTm02_c1m_reg2.0.nc"
    MLD <- nc_open(MLD)
    print(MLD)

there are 7 variables and I only want 'mld'

lon <- ncvar_get(MLD, varid = "lon")
lat <- ncvar_get(MLD, varid = "lat")

summary(lon)
summary(lat)



MLD$dim$time$units

MLD_1.array <- ncvar_get(MLD, "mld")
dim(MLD_1.array)

length(lon)
length(lat)

ndvi.slice <- MLD_1.array[, , 12] 
dim(ndvi.slice)

mld.vec.long <- as.vector(MLD_1.array)
length(mld.vec.long)

nlon <- dim(lon)
nlat <- dim(lat)
lonlat <- expand.grid(lon, lat)

t <- ncvar_get(MLD, "time")
tunits <- ncatt_get(MLD, "time", "units")
nt <- dim(t)

dname <- "mld" 

tmp.mat <- matrix(mld.vec.long, nrow = nlon * nlat,    ncol = nt)
dim(tmp.mat)

head(na.omit(tmp.mat))

lonlat <- expand.grid(lon, lat)
tmp.df02 <- data.frame(cbind(lonlat, tmp.mat))
 names(tmp.df02) <- c("lon", "lat", "Jan", "Feb", "Mar", "Apr", "May", 
"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

library(reshape)
tmp.df03 <- melt(tmp.df02, id=c("lat","lon"))

tmp.df04 <- subset(tmp.df03, lat >= "-35" & lat <="-80")

 tmp.df04[tmp.df04 ==1.000e+09] <- NA
 tmp.df04
 summary(tmp.df04)


 target = c(-180, 180, -90, -20)
 w <- crop(rgeos::gBuffer(spTransform(countriesLow, CRS(pprj)), width = 0))

I can plot this using ggplot after extracting the data into a .csv, removing the mask values, reimporting it then replotting it (it takes a while...).

Is there a way I can extract one variable (mld) then find the maximum value for each month (time) and plot that? I know my code is really messy and choppy...

any help would be much appreciated! thank you!

@dww: these are not (regular) raster data, so `brick(MLD, layer = 'mld')` gives this error: `cells are not equally spaced; you should extract values as points` — Robert Hijmans, Oct 28 '18 at 22:59
Thanks @RobertHijmans. Is there a better way I could do this? — , Oct 29 '18 at 20:46
You can also try using the new package `stars`. It has many great features and integrates well with `sf`. It is, however, still in constant development, so things can change. — AF7, Nov 05 '18 at 08:02

rar · Answer 1 · 2018-10-30T07:19:27.090

Hope this helps!

library(ncdf4)
library(raster)

#Reading netcdf file and extracting lat, lon and variable
MLD <- nc_open("D:/Personal/test/mld_DReqDTm02_c1m_reg2.0.nc")
lon <- ncvar_get(MLD, varid = "lon")
lat <- ncvar_get(MLD, varid = "lat")
mld <- ncvar_get(MLD, "mld")

# using 'raster' package to read the temporal values of variable into a raster
#r1<-flip(raster(t(matrix(mld[,,1], nrow = 180,    ncol = 90))),direction="y")
e<-extent(min(lon),max(lon),min(lat),max(lat))
#extent(r1)<-e
R<-stack()
# creating raster stack of the time series data
for(i in 1:12){
  r1<-flip(raster(t(matrix(mld[,,i], nrow = 180,    ncol = 90))),direction="y")
  extent(r1)<-e
  R<-stack(R,r1)
} # IMPORTANT: This raster stack could be additionally cropped to extract the user's area of interest.
plot(R)

# Extracting max/min values for each time (raster layer) into a dataframe.
df<-data.frame(Months=month.abb)
df$Months <- factor(df$Months, levels = df$Months)
df$months<-c(1:12)
df$MLD_max<-maxValue(R)
df$MLD_min<-minValue(R)
# using 'ggplot2' package
ggplot(df)+geom_point(aes(Months,MLD_min))+geom_line(aes(months,MLD_min))

score 1 · Answer 2 · answered Oct 30 '18 at 16:59

The reason why brick(f) failed is that there is one latitude that is not regular (the last one, which should have been 90, or absent)

library(ncdf4)
library(raster)
MLD <- nc_open("c:/temp/mld_DReqDTm02_c1m_reg2.0.nc")
lat <- ncvar_get(MLD, varid = "lat")
tail(lat)
#[1] 80.0 82.0 84.0 86.0 88.0 89.5

As it is only one row (the most northern one), it is perhaps OK to ignore that, and do as @rar shows. Here with a slightly different approach

mld <- ncvar_get(MLD, "mld")
a <- aperm(mld, c(2,1,3))
b <- flip(brick(a), 'y')
extent(b) <- c(0,360,-90,90)
b <- rotate(b)
b <- reclassify(b, cbind( 1e+09, NA))
names(b) = month.abb
plot(b,1)

You can now use crop

bb <- crop(b, extent(-180,180,-90,35))
plot(bb,1)

Or extract by coordinates

extract(b, cbind(0,0))
#          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep      Oct      Nov      Dec
#[1,] 27.80343 23.17406 25.27729 19.04638 16.96231 18.96719 18.79073 21.04775 30.04144 38.48611 27.82375 20.18804
>

thank you!! that was incredible helpful! I tried extracting the entire dataset to make a box plot but extract(b) gave me an error: Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘extract’ for signature ‘"RasterBrick", "missing"’ Is there a way I can extract the entire dataset (lat, lon, time and value?) — , Nov 05 '18 at 15:08
You can do `values(b)` for the values. See `rasterToPoints` if you want lon/lat included (or do it yourself with `xyFromCell`. In this case you could consider `boxplot(b)` — Robert Hijmans, Nov 05 '18 at 18:57

plotting netcdf data for selected region

2 Answers2