I am trying to plot a Hovmoller diagram using NetCDFs from Copernicus database. I am trying to plot it for one specific point (lon,lat). When I create a database for this with all depth values and all time values (I have monthly data for multiple years), it is consistently showing an error when I convert it to a raster. I have tried different packages and even different NetCDFs, yet I get the same problem. The code and error is below. ******* THE ERROR *******
> DataRast <- raster(SelData) Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘raster’ for signature ‘"numeric"’
******* THE CODE ********
FileNC <- "cmems_mod_glo_phy_my_0.083deg_P1M-m_1744109649217.nc"
#Requested vertical level for analysis
ReqLev <- 1 # 1 = surface
#Location (lat and lon) of Hovmoller diagram
ReqLat <- 10 # Degrees North
ReqLon <- 40 # Degrees East
#2. Data Access
#Opening NETCDF, using RNetCDF package
DataFile <- open.nc(FileNC)
print.nc(DataFile)
#Selecting working variable
VarName <- 'thetao' # SST
#Accessing Data
Lat <- var.get.nc(DataFile, 'latitude')
print(Lat)
Lon <- var.get.nc(DataFile, 'longitude')
print(Lon)
Depth <- var.get.nc(DataFile, 'depth')
print(Depth)
Time <- var.get.nc(DataFile, 'time')
print(Time)
TimeUnit <- att.get.nc(DataFile, "time", "units") #you get the units
Date <- utcal.nc(TimeUnit, Time, "c") # Transform times to POSIXct
DateCh <- format(Date, format="%Y-%m")
print(DateCh)
#Reading variable attributes
VarLongName <- att.get.nc(DataFile, VarName, "long_name")
VarUnits <- att.get.nc(DataFile, VarName, "units")
Credit <- att.get.nc(DataFile, 'NC_GLOBAL', 'credit') #to include in the plots
#checking the order of the code
var.inq.nc(DataFile, 'thetao')
dim.inq.nc(DataFile, 0) # Dimension 0 is time
dim.inq.nc(DataFile, 1) # Dimension 1 is depth
dim.inq.nc(DataFile, 2) # Dimension 2 is lat
dim.inq.nc(DataFile, 3) # Dimension 3 is lon
#Plotting
#Closest lat and lon indices
ReqLatI <- which.min(abs(Lat-ReqLat))
ReqLonI <- which.min(abs(Lon-ReqLon))
#requesting data from file
#changed to fit the arrangement pattern of lon,lat,depth,time
SelData <- var.get.nc(DataFile, VarName,
start=c(ReqLonI, ReqLatI, 1, 1),
count=c( 1, 1, NA, NA))
#Making dataframe for plot
Grid <- meshgrid(Time, Depth)
TimeRast <- raster(Grid$X)
names(TimeRast) <- 'Time'
DepthRast <- raster(Grid$Y)
names(DepthRast) <- 'Depth'
DataRast <- raster(SelData)
names(DataRast) <- 'Data'
DepthDF <- as.data.frame(DepthRast)
DataDF <- as.data.frame(DataRast)
TimeDF <- as.data.frame(TimeRast)
TimeDF <- mutate(TimeDF, DateCT=utcal.nc(TimeUnit, Time, "c"))
TimeDF <- mutate(TimeDF, Date=as.Date(DateCT))
DataDF <- bind_cols(TimeDF, DepthDF, DataDF)
head(DataDF)
#plot
ggplot(data=DataDF, aes(x=Date, y=Depth, z=Data)) +
geom_contour_filled() +
geom_point(shape=3, color='azure3', size=0.01) +
coord_cartesian(expand=FALSE) +
scale_y_reverse() +
labs(fill = VarUnits) +
theme_light() +
scale_x_date(date_breaks = "1 year",
date_labels = "%Y") +
labs(title = VarLongName,
subtitle = paste('Hovemöller Diagram (Lat=', Lat[ReqLatI],
'º, Lon=', Lon[ReqLonI],'º)', sep=''),
x = "Date",
y = "Depth",
caption = paste('Source:', Credit))
I also tried converting the 'SelData' into a matrix and then proceed. But I recieve an error and the matrix that is generated has only NA values.
The code I used to convert it into a matrix -->
DataRast<- raster(as.matrix(SelData))
The netCDF I am using has only one depth value (i.e, the surface) and monthly data from year 1993 - 2021/06
It is unclear what you mean by the term "raster". It most likely refers to the older raster package, which is for geospatial analysis. That won't really help you here because you have 4D data and you want to plot non-spatial data.
You are using the RNetCDF package to read the netCDF file. That is almost heroic because it is so painful. While RNetCDF is very good at what it does, it lacks an easy interface. You are much better off using the ncdfCF package, which is built on top of RNetCDF. You should also install the data.table package.
library(ncdfCF)
library(data.table)
# This file is slightly different from yours. It is for the Mediterranean, but
# otherwise it has the same structure.
fn <- "~/path/to/cmems_mod_med_phy-tem_anfc_4.2km_P1M-m_1743406956144.nc"
# Open the data set
ds <- open_ncdf(fn)
#> Warning: Unmatched `coordinates` value 'lat' found in variable 'bottomT'.
#> Warning: Unmatched `coordinates` value 'lon' found in variable 'bottomT'.
#> Warning: Unmatched `coordinates` value 'lat' found in variable 'thetao'.
#> Warning: Unmatched `coordinates` value 'lon' found in variable 'thetao'.
# These warnings are related to a problem in the production of the file:
# The coordinates are called "latitude" and "longitude", as you can see below.
# Get the data variable from the data set
(thetao <- ds[["thetao"]])
#> <Variable> thetao
#> Long name: Sea temperature
#>
#> Axes:
#> id axis name long_name length values unit
#> 2 X longitude Longitude 338 [7.000002 ... 21.04167] degrees_east
#> 1 Y latitude Latitude 307 [33.229168 ... 45.979168] degrees_north
#> 3 Z depth Depth 140 [1.018237 ... 5646.199219] m
#> 0 T time Time 13 [2024-01-01 ... 2025-01-01] seconds since 1970-01-01 00:00:00
#>
#> Attributes:
#> id name type length value
#> 0 units NC_STRING 1 degrees_C
#> 1 _FillValue NC_FLOAT 1 NaN
#> 2 standard_name NC_STRING 1 sea_water_potential_temperature
#> 3 long_name NC_STRING 1 Sea temperature
#> 4 coordinates NC_STRING 1 time depth lat lon
#> 5 missing_value NC_DOUBLE 1 1.00000002004088e+20
#> 6 valid_max NC_FLOAT 1 40
#> 7 valid_min NC_FLOAT 1 1
(dt <- thetao$profile(longitude = 10, latitude = 40, .names = "Cala_Gonone", .as_table = TRUE))
#> Warning in `[.data.table`(dt, , `:=`(.value, private$values)): 13 column matrix
#> RHS of := will be treated as one vector
#> depth time longitude latitude .variable .value
#> <num> <char> <num> <num> <char> <num>
#> 1: 1.018237 2024-01-01 10 40 Cala_Gonone 16.06493
#> 2: 3.165747 2024-01-01 10 40 Cala_Gonone 16.06098
#> 3: 5.464963 2024-01-01 10 40 Cala_Gonone 16.05880
#> 4: 7.920377 2024-01-01 10 40 Cala_Gonone 16.05823
#> 5: 10.536604 2024-01-01 10 40 Cala_Gonone 16.05782
#> ---
#> 1816: 5224.744629 2025-01-01 10 40 Cala_Gonone NaN
#> 1817: 5328.593750 2025-01-01 10 40 Cala_Gonone NaN
#> 1818: 5433.460938 2025-01-01 10 40 Cala_Gonone NaN
#> 1819: 5539.333496 2025-01-01 10 40 Cala_Gonone NaN
#> 1820: 5646.199219 2025-01-01 10 40 Cala_Gonone NaN
With this data.table you can do the plotting as you please. To get the coordinates along the axes do depths <- ds[["depth"]]$coordinates and months <- ds[["time"]]$time()$format("%Y-%m").