Reading ROMS files in R (R-project)

[wilkin]

Can anyone point me to help on code/utilities that will correctly interpret the ROMS ocean_s_coordinate_g1 and ocean_s_coordinate_g2 (defined in CF conventions) when reading ROMS output netcdf files in "R" (https://www.r-project.org)?

[mdsumner]

At base you can use the ncdf4 package, the system deps are in the source tar ball:

https://cran.rstudio.com/web/packages/ncdf4/index.html

A bit of simplistic pseudo code.

Code: Select all

## install.packages("ncdf4") ## assuming system deps are installed
library(ncdf4)
rfile <- "/path/to/roms/r0001.nc"
nc <- nc_open(rfile)
## standard var-get, optional arguments start/count
myvar <- nvar_get(nc, "ocean_s_coordinate_g1") 

But, I don't know exactly what "ocean_s_coordinate_g1" means here. My ROMS output examples don't have this as an explicit variable. Is this the depth of every cell in [XI,ETA, W]?

From what I've seen we'd need to take the "h" variable and multiply it out by the "Cs_r" stretching curve to get that (??)

I'm writing R tools for ROMS, still in early development but I can help you get what you need. I prefer to use the ncdf4 tools via the raster package, which takes a bit of care (since it assumes a regular affine-grid, like GIS) but is really powerful.

Warning: I'm not proficient in either ROMS or CF, so I probably will use the wrong terms here. I'm pretty au fait with NetCDF generally though, and definitely a native in R. We're working on tools for integrating ROMS output with ecosystem models.

[ocefpaf]

Can anyone point me to help on code/utilities that will correctly interpret the ROMS ocean_s_coordinate_g1 and ocean_s_coordinate_g2 (defined in CF conventions) when reading ROMS output netcdf files in "R" (https://www.r-project.org)?

I do not know of any R tool that does that, but if R is a necessity, you might want to investigate the possibility of using Jupyter notebooks to "glue" R and Python. There are many Python packages that read ROMS files and interpret the vertical coordinate, then you could call R on the results without saving/loading data to change languages.

Here is a quick example:

http://nbviewer.jupyter.org/gist/ocefpa ... a46987b6b4

If people are interested in trying that route I can help with the installation instructions to reproduce that notebook.

Edit: The same can be done with Matlab/Octave and R. The Jupyter notebook has kernels for many languages. See: https://github.com/ipython/ipython/wiki ... -languages

[shchepet]

Speaking CF-conventions....

Would it be hard to adapt and universally enforce a convention that time dimension and
the associated timing variable in any ROMS-related netCDF file have the same name?

This should apply to ROMS code itself as well as all pre- and post-processing software,
Matlab, and Python scripts.

Doing so abolishes the need to keep track of timing variable name and simplifies bookkeeping.

As of right now cat varinfo.dat | grep time > xxx results in a 584-line long output xxx.

If the convention is adopted, a sequence of calls

Code: Select all

! Determine netCDF variable ID for timing variable by assuming that 
! the file is "CF-compliant", so the name of timing variable should 
! match the name of the time dimension (as it should be a "coordinate 
! variable"), while the time dimension itself is expected to be the 
! last dimension of the primary variable.

               ierr=nf_inq_varndims(ncid, varid, ndims)
               if (ierr == nf_noerr) then
                 ierr=nf_inq_vardimid(ncid, varid, dimids)
                 if (ierr == nf_noerr) then
                   ierr=nf_inq_dimname(ncid, dimids(ndims), tname)
                   if (ierr == nf_noerr) then
                     ltvr=lenstr(tname)
                     ierr=nf_inq_varid(ncid, tname(1:ltvr), timevarid)

becomes universally applicable and there is no need to know the name of timing
variable a priori.

[JDTilley]

In addition to the netcdf4 library, rgdal works great with many spatial datatypes.