Swap axis in raster brick

Question

Using R's raster package, I have a brick obtained from a file, with the following ncdump header (I show a small example file, the actual file is much bigger):

dimensions:
        lon = 2 ;
        lat = 3 ;
        time = UNLIMITED ; // (125000 currently)
variables:
        float lon(lon) ;
                lon:standard_name = "longitude" ;
                lon:long_name = "longitude" ;
                lon:units = "degrees_east" ;
                lon:axis = "X" ;
        float lat(lat) ;
                lat:standard_name = "latitude" ;
                lat:long_name = "latitude" ;
                lat:units = "degrees_north" ;
                lat:axis = "Y" ;
        double time(time) ;
                time:standard_name = "time" ;
                time:long_name = "Time" ;
                time:units = "seconds since 2001-1-1 00:00:00" ;
                time:calendar = "standard" ;
                time:axis = "T" ;
        short por(time, lat, lon) ;
                por:_FillValue = 0s ;
                por:missing_value = 0s ;

And in R:

class       : RasterBrick 
dimensions  : 3, 2, 6, 125000  (nrow, ncol, ncell, nlayers)
resolution  : 0.008999825, 0.009000778  (x, y)
extent      : 6.4955, 6.5135, 44.0955, 44.1225  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : /home/clima-archive/afantini/chym/chym_output/test.nc 
names       : X0, X3600, X7200, X10800, X14400, X18000, X21600, X25200, X28800, X32400, X36000, X39600, X43200, X46800, X50400, ... 
z-value     : 0, 449996400 (min, max)
varname     : por

However, for faster access and higher compression two of the file dimensions have been swapped, so that the chunking is better for the kind of use that we need. So the file will be like this (link to download the 1MB file):

dimensions:
        lon = UNLIMITED ; // (2 currently)
        lat = 3 ;
        time = 125000 ;
variables:                                                                                                                                                                                                                            
        float lon(lon) ;                                                                                                                                                                                                              
                lon:standard_name = "longitude" ;                                                                                                                                                                                     
                lon:long_name = "longitude" ;                                                                                                                                                                                         
                lon:units = "degrees_east" ;                                                                                                                                                                                          
                lon:axis = "X" ;                                                                                                                                                                                                      
        float lat(lat) ;                                                                                                                                                                                                              
                lat:standard_name = "latitude" ;                                                                                                                                                                                      
                lat:long_name = "latitude" ;                                                                                                                                                                                          
                lat:units = "degrees_north" ;                                                                                                                                                                                         
                lat:axis = "Y" ;                                                                                                                                                                                                      
        double time(time) ;                                                                                                                                                                                                           
                time:standard_name = "time" ;                                                                                                                                                                                         
                time:long_name = "Time" ;                                                                                                                                                                                             
                time:units = "seconds since 2001-1-1 00:00:00" ;                                                                                                                                                                      
                time:calendar = "standard" ;                                                                                                                                                                                          
                time:axis = "T" ;                                                                                                                                                                                                     
        short por(lon, lat, time) ;
                por:_FillValue = 0s ;
                por:missing_value = 0s ;

And in R:

class       : RasterBrick 
dimensions  : 3, 125000, 375000, 2  (nrow, ncol, ncell, nlayers)
resolution  : 3600, 0.009000778  (x, y)
extent      : -1800, 449998200, 44.0955, 44.1225  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /home/clima-archive/afantini/chym/chym_output/test_swapped.nc 
names       : X6.5, X6.50899982452393 
degrees_east: 6.5, 6.50899982452393 
varname     : por

As you can see the file is opened as if the number of columns is 125000. I want to swap the number of columns with the number of layers, without reading all the data. I think from the raster manual I should use layer or lvar, since:

layer: integer. The layer (variable) to use in a multi-layer file, or the layer to extract from a RasterStack/Brick or SpatialPixelsDataFrame or SpatialGridDataFrame. An empty RasterLayer (no associated values) is returned if ‘layer=0’

.......

‘lvar’: integer > 0 (default=3). To select the 'level variable' (3rd dimension variable) to use, if the file has 4 dimensions (e.g. depth instead of time)

But this does not seem to work, setting e.g. layer="time", as it changes nothing.

How can I do this?

Answer 1

If you don't mind reshaping AFTER opening/reading, I think you could just read the data in a variable using ncdf4library and then transpose it. Something like:

nc   <- nc_open(*your_nc_file*)
data <- ncvar_get(nc, por)     # "por" is the name of your variable, right ? 
data_new <- aperm(data, c(1,3,2)) # "transpose" the matrix

a possible problem could be that data_new is no longer a raster* object, but you can recreate one from it easily enough.

HTH,

Lorenzo