光栅砖中的交换轴

Question

使用 R 的 raster 包，我有一个从文件中获取的brick，带有以下ncdump 标头（我展示了一个小示例文件，实际文件要大得多）：

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 ;

在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

但是，为了更快的访问和更高的压缩率，两个文件尺寸已被交换，因此分块更适合我们需要的用途。所以文件会是这样的（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 ;

在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

如您所见，文件打开时好像列数为 125000。我想将列数与层数交换，而不读取所有数据。我认为从光栅手册中我应该使用layer 或lvar，因为：

层：整数。在多层文件中使用的层（变量）， 或从 RasterStack/Brick 中提取的图层或 SpatialPixelsDataFrame 或 SpatialGridDataFrame。一个空 如果“layer=0”，则返回 RasterLayer（无关联值）

.......

‘lvar’：整数 > 0（默认值=3）。选择“级别变量” （第 3 维变量）使用，如果文件有 4 维 （例如深度而不是时间）

但这似乎不起作用，例如设置layer="time"，因为它没有任何改变。

我该怎么做？

Answer 1

如果您不介意在打开/阅读后进行整形，我认为您可以使用 ncdf4library 读取变量中的数据，然后转置它。比如：

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

一个可能的问题是data_new 不再是raster* 对象，但您可以轻松地从中重新创建一个对象。

HTH，

洛伦佐