Stars: Errors reading NetCDF files coming from climate model

[adrfantini]

Thanks for stars, it looks extremely promising. I'm using the latest github version.

I've tested it on a few files coming from our climate model, RegCM, which produces well-described NetCDF data as output.
The data is on a Lambert Conformal Conic grid and the CRS is described, as per CORDEX guidelines, by a crs variable.

Two kind of errors are recurring for different files. The first (input file link):

> sts <- st_stars("CORDEXRUN_4_STS.1990010100.nc")
xlon, xlat, mask, topo, ps, tsmax, tsmin, prmax, pr, sund, psmin, psavg, mrros, mrro, tasmax, tasmin, tas, sfcWindmax, time_bnds, 
Error in check_equal_dimensions(dots) : 
  object 1 and 5 have different dimensions
In addition: Warning messages:
1: In CPL_read_gdal(x, options, driver, TRUE) :
  GDAL Message 1: dimension #1 (time_bounds) is not a Longitude/X dimension.
2: In CPL_read_gdal(x, options, driver, TRUE) :
  GDAL Message 1: dimension #0 (time) is not a Latitude/Y dimension.

The second (input file link):

> rad <- st_stars("CORDEXRUN_4_RAD.1990010100.nc")
xlon, xlat, mask, topo, ps, rsns, rsnl, rtnscl, rsnscl, rtnlcl, rsnlcl, rts, rsnt, prw, clwvi, clivi, rtl, clh, clm, cll, cl, Error: inherits(value, "units") || inherits(value, "symbolic_units") is not TRUE

You can check the structure of those files using ncdf4::nc_open() from R or ncdump -h from the command line.

[edzer]

Thanks for expressing interest!

The first problem is caused by arrays not having the same dimensions; use the sub argument e.g.

x = st_stars("CORDEXRUN_4_STS.1990010100.nc", sub = 9:12)

to select a subset of arrays that has the same dimensions.

[adrfantini]

Thanks Edzer, I assume that opening files with variables with different dimension might be a problem? Is this planned, or that's a feature that will never come?
This is basically what NetCDF is about: a very flexible way of defining several datasets with any number/kind of dimensions in a single file.

If this is not readable in a single stars call (I guess the infinite flexibility of NetCDF would make that quite hard), I would then suggest that with NetCDF data st_stars() should accept variable names in sub=, just like raster() accepts an optional arg variable=. Or something along those lines. Having to check which number corresponds to which variable is quite annoying.

Also, once we have x, stars might want to have a method for plot(), like raster does, so that one can take a quick glance at the data, but that is another matter.

[edzer]

Yes I agree that sub should accept names; will take care of that.

stars objects are thought of as single, potentially multi-attribute, array (even if implemented as a list of arrays) with shared dimension properties; netcdf files are collections of different-dimensioned arrays. You can read in different arrays by two calls to st_stars:

> dim(st_stars("CORDEXRUN_4_RAD.1990010100.nc", sub = 1))
 x  y 
65 65 
> dim(st_stars("CORDEXRUN_4_RAD.1990010100.nc", sub = 9))
   x    y time 
  65   65  120 

similar to how raster reads either layers, or bricks. And yes, a quick overview plot is also high on my list; there is now image which only shows a single (first) slice.

[tim-salabim]

Re: plotting, see also mapview issue https://github.com/r-spatial/mapview/issues/115 currently working on supporting stars. Single layer image is already possible via addStarsImage (like addRasterImage), however, we need a longlat image and specifically transform the bounds to epsg 3857 - as I have not figured out how to transform a stars object so far.

[edzer]

What do you mean by a longlat image? Many netcdf files (like the ones above) come with a long and a lat array, giving {long,lat} for each (i,j); I consider these somewhat obsolete as we have the proj.4 string and can compute long/lat for each (x,y) from that. Would you like to have them imported, or computed?

As for transform, I can see two ways:

1. convert, as is, to sf (as polygons), then transform that (as shown in the blog with the affine grid); should be possible now
2. warp the grid (with some strategy) to a grid in 3857 (not yet implemented; maybe use gdalwarp, available through stars::gdal_utils?)

[tim-salabim]

I mean that for now, it only works for stars objects with crs 4326. An equivalent of projectRaster would solve that.

[adrfantini]

When I need to plot projected (non-longlat, crs!=4326) netcdf files "properly" with ggplot, I usually extract the polygons (with rasterToPolygons, or now with spex::polygonize since @mdsumner pointed it out - thanks!), transform those to sf, and plot the single polygons. This is extremely slow, but is the only way currently to properly plot projected gridded data without any alteration AFAIK.

Another option is to reproject, which I usually do with CDO outside of R since it can use several good weighting methods. gdalwarp is another option. However, the reprojecting ALWAYS introduces small mistakes and differences, which are unnoticeable most of the time, but can be important in some corner cases (which I experience all the time :().

In the climate world, common tools for plotting NetCDF gridded data are GrADS (which afaik automatically reprojects data on it's internal lat-lon grid) and NCL (which can deal with projected coordinates and is also quite fast). However I have no idea how NCL manages it.

[edzer]

@adrfantini your second dataset ("RAD") poses a problem for GDAL for those variables that vary over elevation:

$gdalinfo NETCDF:"CORDEXRUN_4_RAD.1990010100.nc":cl | head -170 | tail -7
  NETCDF_DIM_EXTRA={time,kz}
  NETCDF_DIM_kz_DEF={18,6}
  NETCDF_DIM_kz_VALUES={351390,351396,351402,351408,351414,351420,351426,351432,351438,351444,351450,351456,351462,351468,351474,351480,351486,351492,351498,351504,351510,351516,351522,351528,351534,351540,351546,351552,351558,351564,351570,351576,351582,351588,351594,351600,351606,351612,351618,351624,351630,351636,351642,351648,351654,351660,351666,351672,351678,351684,351690,351696,351702,351708,351714,351720,351726,351732,351738,351744,351750,351756,351762,351768,351774,351780,351786,351792,351798,351804,351810,351816,351822,351828,351834,351840,351846,351852,351858,351864,351870,351876,351882,351888,351894,351900,351906,351912,351918,351924,351930,351936,351942,351948,351954,351960,351966,351972,351978,351984,351990,351996,352002,352008,352014,352020,352026,352032,352038,352044,352050,352056,352062,352068,352074,352080,352086,352092,352098,352104}
  NETCDF_DIM_time_DEF={120,6}
  NETCDF_DIM_time_VALUES={351390,351396,351402,351408,351414,351420,351426,351432,351438,351444,351450,351456,351462,351468,351474,351480,351486,351492,351498,351504,351510,351516,351522,351528,351534,351540,351546,351552,351558,351564,351570,351576,351582,351588,351594,351600,351606,351612,351618,351624,351630,351636,351642,351648,351654,351660,351666,351672,351678,351684,351690,351696,351702,351708,351714,351720,351726,351732,351738,351744,351750,351756,351762,351768,351774,351780,351786,351792,351798,351804,351810,351816,351822,351828,351834,351840,351846,351852,351858,351864,351870,351876,351882,351888,351894,351900,351906,351912,351918,351924,351930,351936,351942,351948,351954,351960,351966,351972,351978,351984,351990,351996,352002,352008,352014,352020,352026,352032,352038,352044,352050,352056,352062,352068,352074,352080,352086,352092,352098,352104}
  time#bounds=time_bnds
  time#calendar=gregorian

this means that the values associated with kz are identical in value and in number to those of time; this cannot be the case. I can see in the ncdump that there is a time(time) array (vector), but not an kz(kz). I am not sure whether this is a problem of the gdal netcdf driver, or of the netcdf file (which doesn't seem to be CF-compliant!). I think there is little hope we can get to the elevation values using the current approach of stars. Could you pls verify that this file is CF-1 compliant?

[edzer]

@adrfantini for raster polygonizing, would you mind comparing performance between your options and stars::st_as_sf(x, as_points = FALSE)? It should be pretty fast, written in c++.

[edzer]

(or st_xy2sfc, which keeps the info in a stars object)

[tim-salabim]

@edzer it's a bit of weird approach that leaflet takes for displaying raster images. The grid itself needs to be in 4326, but the bounds (the 4 corner points) need to be passed in 3587. We can solve that internally if we have a st_bbox - which I think stars already provides - and a st_transform. Then we could easily provide a method for one band and hence multiple bands.

[edzer]

Shall I add an st_transform method for bbox objects?

[adrfantini]

@edzer
The level info is stored not into kz(kz) but into sigma(kz) since this is a sigma pressure coordinates model. This is CORDEX compliant which should be 100% CF-1 compliant. Are you sure that CF-1 conventions assume a variable of the same name of the dimension describing it? I'm just a poor PhD student (doing something else entirely), this is a bit over my head :)

I'll test the polygonize functions you provided, thanks.

[edzer]

I can see that, and I don't now. Sorry, I overlooked the

NC_GLOBAL#Conventions=CF-1.4

will have to ask on gdal-dev, then.

[tim-salabim]

@edzer maybe a better approach would be to include a as_sfc argument in st_bbox? That way the user can decide whether a valid sf geometry is returned and transform as needed. In any case I personally am fine with st_transform(st_as_sfc(st_bbox())))

[adrfantini]

@edzer I will ask our modeller tomorrow too.

[edzer]

@tim-salabim : that would break type safety ; yes, sf_as_sfc(bbox) is simple enough.

[tim-salabim]

OK, then the only thing is st_transform

[edzer]

as in: "the only thing still to do"? If so, st_transform of what?

[tim-salabim]

A stars object. Equivalent to projectRaster.

[edzer]

Ah, OK.

[edzer]

@adrfantini problem reported here: https://lists.osgeo.org/pipermail/gdal-dev/2017-November/047800.html

[adrfantini]

@edzer despite installing the latest stars from github I seem to be unable to have stars::st_as_sf() ? :/

I used devtools::install_github("r-spatial/stars"), maybe the current version is somewhere else? Sorry, I'm a noob.

Error: 'st_as_sf' is not an exported object from 'namespace:stars'

[edzer]

Sorry, you just use st_as_sf on a stars object, after loading pkg stars, that's all. My fault.

[adrfantini]

@edzer see the performance of different polygonizing methods in #13

[adrfantini]

@edzer I looked at CF-1.4 and it looks like it does not say much about coordinate variables. However the CF extend COARDS conventions, and they say:

Coordinate Variables:
1-dimensional netCDF variables whose dimension names are identical to their variable names are regarded as "coordinate variables" (axes of the underlying grid structure of other variables defined on this dimension).

So I'd say that our model is currently non-compliant.

[edzer]

The other interpretation is that kz is not a coordinate variable. In any case, the bug in GDAL has been reported (and, it seems, fixed): https://trac.osgeo.org/gdal/ticket/7165

[edzer]

@tim-salabim we now have an st_transform method for stars objects, doing warp/regrid for grids. It needs testing, and does only method='ngb', in raster-speak.

[adrfantini]

If in the future additional methods could be provided it'd be nice. Gdal supports different warping methods (seehttp://www.gdal.org/gdalwarp.html), and so does CDO (bilinear, bicubic, conservative, distance weighted, largest area...), all of them have use cases.

[edzer]

@adrfantini I agree, but note that raster always only had 2, and I doubt how many people ever chose for the non-default... stars provides gdal_utils, including warp, but that does file -> file, really, without having data in R; it does have all the options, obviously. Not sure how much value there is in reimplementing all of them.

@tim-salabim https://github.com/r-spatial/stars/commit/520d2a12f63c4ec3417e0c42289a781e564a5d54 speeds up st_tranform of images by a factor 10 (but now needs fresh install of sf from github, as I removed the proj dependency of stars).

[adrfantini]

@edzer yeah I understand that's not a priority right now. CDO does have some options iirc that GDAL does not tho.

[edzer]

x = st_stars("CORDEXRUN_4_RAD.1990010100.nc", sub = "cll")
plot(x, names = NULL, col = grey((5:10)/10))

[edzer]

Update: with more modern stars version you'd need

plot(x, main = NULL, col = grey((5:10)/10), breaks = "equal")

for the plot above.