length(geotransform) == 6 is not TRUE when setting dimension as sfc column

[adrfantini]

I came upon this issue while trying to assign an sfc column dimension to a stars object derived from a point-timeseries netCDF file. Ideally sooner or later stars will understand these files directly (#30), but in the meantime I am trying to manually force the correct dimensions.
The object has two dimensions, station number and time, and I would like to assign lat-lon coordinates to each station.

Use this file (2.9MB) for testing.

library(ncdf4)
library(stars)

fin = '201001.nc'

nc = nc_open(fin)
lon = ncvar_get(nc, 'lon')
lat = ncvar_get(nc, 'lat')
nc_close(nc)

s = read_ncdf(fin)
sfc = st_as_sf(data.frame(lon = lon, lat = lat), coords = 1:2, crs = 4326)$geometry
s2 = st_set_dimensions(s, 'stat', sfc, point = TRUE)

Returns this error:

Error in xy_from_colrow(cbind(0, seq(x$from, x$to) - 1 + where), geotransform) : 
  length(geotransform) == 6 is not TRUE

[edzer]

Thanks! after

attr(attr(s, "dimensions"), "raster") = stars:::get_raster(dimensions = rep(NA_character_,2))

you can do

s2 = st_set_dimensions(s, 'stat', sfc, point = TRUE)

As you've noted earlier, we need to make the first step easier, but we also need to adapt read_ncdf such that it doesn't label the stat and time dimensions as x and y raster dimensions, as that is the cause of your issue here. I guess for the latter it should be enough to look at the refsys of time.

[edzer]

This allows you to unset [x] and [y] raster dimensions by

s2 = st_set_dimensions(s, xy = c(NA, NA))

Maybe @dblodgett-usgs or @mdsumner have an idea why it gets set in the first place, by read_ncdf, for this file.

[adrfantini]

Thanks, that's handy.

[mdsumner]

I just set the first two as the raster by default because I wanted something to plot. I am falling behind in what the stars vision is, I can see how this example is supposed to work but I'd probably start again with non raster stuff like this. The basis is there in ncmeta thanks to @dblodgett-usgs

ncmeta::nc_coord_var("~/201001.nc")
# A tibble: 2 x 6
  variable X     Y     Z     T     bounds
  <chr>    <chr> <chr> <chr> <chr> <chr> 
1 time     NA    NA    NA    time  NA    
2 pr       lon   lat   NA    time  NA 

I'll have to have a quiet look and a think about matching this up, it's a good example.

[mdsumner]

Quick update, I think I should not have set raster when coord_dim is not true for either of the first two dimensions, available via 'nc_meta()$dimension'

[dblodgett-usgs]

@mdsumner that would be a good fix for this issue. This is something I want to work on for #30. Trying to figure out how this will fit into the stars dimension data model is top of my list and you'll see me in the issues here and / or having a quiet look and think about it in August.

[mdsumner]

Oh cool, thanks! The thing I get stuck on is all the tail-chasing for what the geometries are, most of the time they will be long-lat and can assume WGS84 etc., but occasionally they might be 3D geometry, or perhaps lat-long, or maybe projected and incomplete metadata etc. But, this is also newly clarifying for me so I'll happily follow your lead ;)

[dblodgett-usgs]

Well, the CF1.8 geometries spec is live now, and it's all there! 😉 But seriously, it's going to be near impossible to figure out unless these kinds of data follow a pretty specific convention. We'll get that stuff in ncmeta and work it into stars and try to infer how things should wire up the best we can.

[mdsumner]

"in ncmeta" that reminds of of something I've wanted to bring up @edzer - there's a lot of time-handling in "read_ncdf" that I think belongs in a separate tool. While I was working on a tidync version of that function I found it sensible to get an object with all the variable and coordinate attributes with the units classes sorted out upfront.

https://github.com/mdsumner/stars/blob/netcdf-dev/R/ncdf.R#L57

That's what I think we need as a general background tool here - I used to think it was impossible, but it's interesting how much does get captured sensibly by churning through the units mill. I got a certain way into it but then got bogged down and torn away by other tasks.

It also might be best to wait until GDAL RFC75 comes through (seems to be going faster than I ever thought) and we might get a lot of stuff for free that way.

[edzer]

I agree that handling time, handling units and handling bounds could all go into ncmeta, but I guess we'd need to get them right first, or did we? I'm not so optimistic that RFC75 will help handling time, units and/or bounds, but one of us should really dive into what Even is doing!

[edzer]

See #555