Stars: read_stars on irregular netcdf converts to a regular grid

[kendonB]

File is:

mrso_Lmon_CMCC-CESM_rcp85_r1i1p1_200001-200412.zip

library(tidyverse)
library(stars)
climate_model_output <- read_stars("mrso_Lmon_CMCC-CESM_rcp85_r1i1p1_200001-200412.nc")  
#> Warning in CPL_read_gdal(as.character(x), as.character(options),
#> as.character(driver), : GDAL Message 1: Latitude grid not spaced evenly.
#> Setting projection for grid spacing is within 0.1 degrees threshold.

Doesn't look like it comes out as irregular

climate_model_output %>% 
  st_dimensions() %>% 
  .$y
#> $from
#> [1] 1
#> 
#> $to
#> [1] 48
#> 
#> $offset
#> [1] 89.01354
#> 
#> $delta
#> [1] -3.708898
#> 
#> $refsys
#> [1] NA
#> 
#> $point
#> [1] NA
#> 
#> $values
#> NULL
#> 
#> attr(,"class")
#> [1] "dimension"

Showing that the raw file is indeed irregular:

ncdf_file <- ncdf4::nc_open("mrso_Lmon_CMCC-CESM_rcp85_r1i1p1_200001-200412.nc")
diff(ncdf_file$dim$lat$vals)
#>  [1] 3.680158 3.701891 3.706801 3.708663 3.709564 3.710067 3.710376
#>  [8] 3.710579 3.710720 3.710821 3.710895 3.710952 3.710996 3.711030
#> [15] 3.711057 3.711079 3.711096 3.711110 3.711121 3.711129 3.711135
#> [22] 3.711140 3.711142 3.711143 3.711142 3.711140 3.711135 3.711129
#> [29] 3.711121 3.711110 3.711096 3.711079 3.711057 3.711030 3.710996
#> [36] 3.710952 3.710895 3.710821 3.710720 3.710579 3.710376 3.710067
#> [43] 3.709564 3.708663 3.706801 3.701891 3.680158

^{Created on 2019-02-21 by the reprex package (v0.2.1.9000)}

[mdsumner]

That's GDAL behaviour, to "assume regular" within constrains and is possibly close enough. With read_ncdf you should get the lat as a "values" dimension.

A note about the context, things I've learned the hard way:

Whether or not that was intended by the creators or it should be treated as regular is social history problem - can you find out from the source experts? (I'll have a closer look soon, but often I find it's not possible to answer and the right way forward simply depends on what you need to do - I see incorrectly mapped northern hemisphere regions in southern-focussed climate output, and get this: it turns out you are supposed to ignore the northern parts of the grid - only Greenland really shows that's it obviously out, but there's no way anyone to be able to tell that - you have to glean it from the culture/community. My data-experience and inclination leads me to hunt down why the entire grid doesn't map properly, but that was the wrong approach!).

[kendonB]

Thanks so much for the insights @mdsumner. My immediate goal is to pull out some basic changes over the century for a single location but for all models. But it seems that some of the models shouldn't be used for this goal!

[kendonB]

Is the long term intention for read_stars to use read_ncdf? And read_ncdf is currently "experimental"?

[mdsumner]

Long term I'd want to be able to choose the GDAL way or the pure NetCDF way - I do think they can be wrapped, but the idioms for reading partial-arrays is quite different and there's considerable unknown work to do that. Any feedback you have about usage with real world data can help, and see https://github.com/r-spatial/stars/issues/62

(GDAL is 2D-layer-bound and offset/srcsize/outsize based - higher dimensions present as "bands", which can be thought of as 2D layers, NetCDF is nD-arrays with offset-start/offset-end indexing, with no within-stride sub-sampling )

read_ncdf is currently much raw-er than read_stars (no metadata, crs interpretation is applied), but is closer to the native format in terms of the interface and will allow some flexibility not provided by the GDAL interface (though the GDAL one can be used in ever more more flexible ways ... this is not a simple topic to summarize neatly).

[adrfantini]

My immediate goal is to pull out some basic changes over the century for a single location but for all models. But it seems that some of the models shouldn't be used for this goal!

In general climate models should not be trusted to provide a reliable representation over a single grid point - I realise it might sometimes be tempting to do so, but it is discouraged. You should at least do an analysis in a small subset around your selected point. Due to topography smoothing and sometimes imprecise coastal boundaries this is especially true in coastal and mountain areas.

[edzer]

I would like to port the time parsing from read_stars to read_ncdf (should be trivial with the units pkg), and also have it read and use coordinate bounds, if they are supplied, since the stars data model caters for them.

[edzer]

In this case, the lat values are not the mid points from the lat_bnds pairs, and working your way back from lat to lat_bnds in stars would not give you back the lat_bnds in the .nc file.

[dblodgett-usgs]

Re time parsing, the RNetCDF function utcal.nc to convert NetCDF time to POSIX is quite good in my experience. Maybe your suggested approach using units to handle the time units would get there too though.

[edzer]

Thanks! RNetCDF also uses udunits2, so that should do the same. Tricky might be to handle 360- and 365-day years, now done through pkg PCICt in read_stars (converts into PCICt, rather than POSIXt).

[dblodgett-usgs]

Ahh... looking at some code now -- looks like about the same difference. I'd port the read_stars over to read_ncdf then.

[edzer]

This change tries to read dimension first from a bounds variable, if present. Using one of the files from @kendonB , we used to get

library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.7.0, GDAL 2.3.2, PROJ 5.2.0
(x = read_ncdf("mrso_Lmon_GFDL-CM2p1_rcp45_r1i1p1_200601-201012.nc"))
# stars object with 3 dimensions and 1 attribute
# attribute(s):
#      mrso           
#  Min.   :   0.0000  
#  1st Qu.:   0.0000  
#  Median :   0.0000  
#  Mean   :  16.3720  
#  3rd Qu.:   0.0275  
#  Max.   :2079.4417  
# dimension(s):
#      from  to      offset delta refsys point                    values    
# lon     1 144 2.22045e-16   2.5     NA    NA                      NULL [x]
# lat     1  90          NA    NA     NA    NA [90] -89.4944,...,89.4944 [y]
# time    1  60          NA    NA     NA    NA  [60] 52976.5,...,54771.5    
st_bbox(x)
#          xmin          ymin          xmax          ymax 
#  2.220446e-16 -8.949438e+01  3.600000e+02  8.949438e+01 

but now have

library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.7.0, GDAL 2.3.2, PROJ 5.2.0
(x = read_ncdf("mrso_Lmon_GFDL-CM2p1_rcp45_r1i1p1_200601-201012.nc"))
# stars object with 3 dimensions and 1 attribute
# attribute(s):
#      mrso           
#  Min.   :   0.0000  
#  1st Qu.:   0.0000  
#  Median :   0.0000  
#  Mean   :  16.3720  
#  3rd Qu.:   0.0275  
#  Max.   :2079.4417  
# dimension(s):
#      from  to offset delta  refsys point
# lon     1 144      0   2.5      NA    NA
# lat     1  90     NA    NA      NA    NA
# time    1  60     NA    NA POSIXct    NA
#                                           values    
# lon                                         NULL [x]
# lat            [-90,-88.98876),...,[88.98876,90) [y]
# time 2006-01-17 12:00:00,...,2010-12-17 12:00:00    
st_bbox(x)
# xmin ymin xmax ymax 
#    0  -90  360   90 

The problem with the grid centers of a rectilinear grid is that you can't recover the original boundaries.