Terra freezes when opening netcdf file

[fabeit]

Hi,

I am having troubles using the 'rast' function to read certain netcdf files, I am not sure which part of these files makes them problematic but I am attaching a sample file.

It's definitely not the size as even a small file causes the rast function to freeze and forces me to terminate R. Some of these files are created using a program coded in fortran, not sure if there's some kind of formatting at the level of netcdf fortran libraries that create a netcdf file that's not compatible with terra.

Any idea what could be the problem? Thanks in advance.

output.zip

[pvanlaake]

There is nothing wrong with your file. It is only that the file is a tad odd and terra might not be the best package to use here (this is from looking at your sample file with a single lon and lat coordinate, i.e. a single point).

The data variables in the file have definitions like float anpop(time_counter, num_ani, lat, lon) with shape [54810, 24, 1, 1]. The GDAL library that terra uses for reading netCDF files tries to turn that into 1,315,440 layers (54,810 x 24) - this is likely the source of the freeze.

Since the data variable is just a temporal profile at a given location, package ncdfCF might be better suited here:

library(ncdfCF)

(ds <- open_ncdf("~/output.nc"))
#> <Dataset> output 
#> Resource   : ~/output.nc 
#> Format     : classic 
#> Collection : Generic netCDF data 
#> Conventions: (not indicated) 
#> Keep open  : FALSE 
#> 
#> Variables:
#>  name      units   data_type axes                           
#>  anpop     Kg/km2  NC_FLOAT  lon, lat, num_ani, time_counter
#>  ugb               NC_FLOAT  lon, lat, num_ani, time_counter
#>  int       Kg/ha   NC_FLOAT  lon, lat, num_ani, time_counter
#>  intlitter         NC_FLOAT  lon, lat, num_ani, time_counter
#>  fat       Kg      NC_FLOAT  lon, lat, num_ani, time_counter
#>  starv             NC_FLOAT  lon, lat, num_ani, time_counter
#>  expend            NC_FLOAT  lon, lat, num_ani, time_counter
#>  estd              NC_FLOAT  lon, lat, num_ani, time_counter
#>  esty              NC_FLOAT  lon, lat, num_ani, time_counter
#>  grass     KgDM/ha NC_FLOAT  lon, lat, time_counter         
#>  leaf      KgDM/ha NC_FLOAT  lon, lat, time_counter         
#>  fruit     KgDM/ha NC_FLOAT  lon, lat, time_counter         
#>  seed      KgDM/ha NC_FLOAT  lon, lat, time_counter

# The "time_counter" axis, an unusual name that might not be recognised by all netCDF readers
# Note the "noleap" calendar: all years have 365 days
ds[["time_counter"]]
#> <Time axis> [5] time_counter
#> Long name  : Time axis 
#> Length     : 54810 (unlimited)
#> Axis       : T 
#> Calendar   : noleap 
#> Range      : 1000-01-02 ... 1150-03-02 (days) 
#> Bounds     : (not set) 
#> 
#> Attributes:
#>  name          type     length value                         
#>  axis          NC_CHAR   1     T                             
#>  standard_name NC_CHAR   4     time                          
#>  long_name     NC_CHAR   9     Time axis                     
#>  calendar      NC_CHAR   6     noleap                        
#>  units         NC_CHAR  30     days since 1000-01-01 00:00:00
#>  actual_range  NC_FLOAT  2     1, 54810

# Open a data variable as an array
(anpop <- ds[["anpop"]]$data())
#> <Data array> anpop 
#> 
#> Values: [1e-04 ... 42448.57] Kg/km2
#>     NA: 548100 (41.7%)
#> 
#> Axes:
#>  axis name         long_name length values                      unit  
#>  X    lon                        1  [-113.029999]               degrees_east
#>  Y    lat                        1  [48.560001]                 degrees_north
#>       num_ani                   24  [1 ... 24]                 
#>  T    time_counter Time axis 54810  [1000-01-02 ... 1150-03-02] days since 1000-01-01 00:00:00
#> 
#> Attributes:
#>  name         type     length value              
#>  units        NC_CHAR  6      Kg/km2             
#>  coordinates  NC_CHAR  7      lon lat            
#>  actual_range NC_FLOAT 2      1e-04, 42448.566406

You can process the data variable as you require. Should you want to convert your processed data into a terra::SpatRaster for plotting or spatial analysis (after reducing the data set to a manageable number of layers), you can call something like r <- anpop$terra().