Writing 3D data frames to netcdf with write_mdim()

[jens-daniel-mueller]

Hi stars team,

To my understanding, write_mdim() is the only function available in R that would allow me to directly write multidimensional grids to netcdf files in a single line of code. This is great, thank you!

I tried to implement a reprex example (below) to write ocean interior data (i.e. including a z-coordinate depth, in addition to lon and lat). The approach works in general, and I'm able to read the same object back into R with read_mdim().

However, when writing with write_mdim() or reading the data with another stars function (read_stars() or read_ncdf()), I receive error messages that suggest that the dimensions lack an indexing variable, and that the z-coordinate depth is not a vertical dimension. I assume that this might be due to my conversion from a tibble to stars object, but couldn't figure out the issue.

Can someone help?
Note: A related question was posted on gis.statexchange and I can take care to reproduce recommendations given here over there.

Thanks much in advance for your help! Cheers, Jens
PS: I hope the reprex works OK, I did not quite know how to set it up for an I/O operation.

library(tidyverse)
library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.9.1, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUE

# create a simplified ocean 3D grid, with depth as z-coordinate
ocean_tibble <- expand_grid(
  lon = seq(-175, 175, 10),
  lat = seq(-85, 85, 10),
  depth = seq(0, 1000, 100)
)

# create an arbitrary variable (named temperature) from coordinates
ocean_tibble <- ocean_tibble %>%
  mutate(temperature = abs((abs(lat)-85))*0.25 - depth/100 + 12)

# convert from tibble to stars object
ocean_stars <- st_as_stars(ocean_tibble,
                           dims = c("lon", "lat", "depth"))

# assign coordinate system to horizontal (lon, lat) coordinates
st_crs(ocean_stars) <- 4326

# plot variable per depth level
ggplot() +
  geom_stars(data = ocean_stars) +
  facet_wrap(~depth) +
  scale_fill_viridis_c()

# write stars object to .nc file
ocean_stars %>%
  write_mdim("data/ocean_raster_write_mdim.nc")
  
#> Warning messages:
#> 1: In CPL_write_mdim(file, driver, dimx, cdl, wkt, xy, root_group_options,  :
#>   GDAL Message 1: Dimension lat lacks a indexing variable
#> 2: In CPL_write_mdim(file, driver, dimx, cdl, wkt, xy, root_group_options,  :
#>   GDAL Message 1: Dimension lon lacks a indexing variable
#> 3: In CPL_write_mdim(file, driver, dimx, cdl, wkt, xy, root_group_options,  :
#>   GDAL Error 6: SetIndexingVariable() not implemented
#> 4: In CPL_write_mdim(file, driver, dimx, cdl, wkt, xy, root_group_options,  :
#>   GDAL Error 6: SetIndexingVariable() not implemented
#> 5: In CPL_write_mdim(file, driver, dimx, cdl, wkt, xy, root_group_options,  :
#>   GDAL Error 6: SetIndexingVariable() not implemented

# reopen .nc file with all available stars functions
ocean_read_stars <- read_stars("data/ocean_raster_write_mdim.nc")

#> Warning messages:
#> 1: In CPL_get_metadata(file, NA_character_, options) :
#>   GDAL Message 1: dimension #0 (depth) is not a Time or Vertical dimension.
#> 2: In CPL_get_metadata(file, NA_character_, options) :
#>   GDAL Message 1: dimension #0 (depth) is not a Time or Vertical dimension.
#> 3: In CPL_read_gdal(as.character(x), as.character(options), as.character(driver),  :
#>   GDAL Message 1: dimension #0 (depth) is not a Time or Vertical dimension.
#> 4: In CPL_read_gdal(as.character(x), as.character(options), as.character(driver),  :
#>   GDAL Message 1: dimension #0 (depth) is not a Time or Vertical dimension.
  
ocean_read_ncdf <- read_ncdf("data/ocean_raster_write_mdim.nc")

#> no 'var' specified, using temperature
#> other available variables:
#>  crs, lon, lat, depth
#> Will return stars object with 7128 cells.
#> Warning messages:
#> 1: In CPL_crs_from_input(x) :
#>   GDAL Error 1: PROJ: proj_create: Error -7 (unknown unit conversion id)
#> 2: In value[[3L]](cond) : failed to create crs based on grid mapping
#> and coordinate variable units. Will return NULL crs.
#> Original error: 
#> Error in st_crs.character(base_gm): invalid crs: +proj=longlat +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs +units=

ocean_read_mdim <- read_mdim("data/ocean_raster_write_mdim.nc")

# plot variable in object reopened with read_mdim
ggplot() +
  geom_stars(data = ocean_read_mdim) +
  facet_wrap(~depth) +
  scale_fill_viridis_c()
  

[edzer]

One thing is that your stars object has no information that depth is actually about the vertical dimension; it has for instance not a unit associated with elevation / height (such as m, but also pressure level). I tried to make write_dim write that the dimension, when called "depth", is VERTICAL (see r-spatial/sf@fbca696) but that doesn't change anything about the written netcdf, which makes me suspect that this is an omission either of the (rather new) GDAL mdim interface, or of the NetCDF GDAL driver code. That would take time to sort out.

Despite the warnings,

ocean_read_stars <- read_stars("data/ocean_raster_write_mdim.nc")

seems to return the object you're looking for? The error in read_ncdf() remains, but we may have to deprecate that given lack of active maintenance.

[jens-daniel-mueller]

Dear Edzer,

thanks much for your immediate reply and trying to modify the write_dim() function such that it automatically write depth as a VERTICAL dimension.

You mention that in my reprex data frame, depth is not defined as a vertical dimension. Would that be possible in stars, and if yes, how? I searched but couldn't find instructions how to achieve this.

Many thanks, Jens
PS: I also tried reading the 3D file with tidync::tidync() now and that works as well. So in terms of content, the written netcdf appears OK.

[edzer]

You mention that in my reprex data frame, depth is not defined as a vertical dimension. Would that be possible in stars, and if yes, how? I searched but couldn't find instructions how to achieve this.

Neither do I; a GDAL-centric solution would be a vertical CRS in WKT, but I'm not sure if such a thing exists; GDAL SpatialReference do have an IsVertical() member, but it is also possible to write WKT CRS with horzontal + vertical CRS, not sure what makes sense (never tried). stars objects register which dimensions are x and y, and whether they are simple feature geometries, implying x+y combined.

I also looked into the CF conventions but that didn't bring much enlightenment. and to which things the GDAL NetCDF driver is sensitive, e.g. here.

But, as a heuristic, if a dimension is not [x] or [y], and has length dimension (e.g. m), a good guess would be to take it as the depth or height, right? Currently stars is not sensitive to this, and I could make it sensitive to it but, as written above, have the feeling that the upstream GDAL mdim interface is not sensitive to it (at least for the NetCDF driver).

[jens-daniel-mueller]

OK, thanks for clarifying. I'll proceed with what we have (which is superb already).