Matlab读取NetCDF指定范围数据报错：horzcat维度不匹配

我需要从NetCDF文件中读取变量tref的指定范围数据：目标经纬度范围（已通过索引筛选）、第24个M（集合成员）、第4个L（预报时效）以及所有S（预报起始时间）。

使用的代码如下：

filename = 'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NCEP-CFSv2/.HINDCAST/.MONTHLY/tref/1+index/.tref/tref/tref/dods';
ncdisp(filename) % so the structure of tref is (X*Y*M*L*S which means lon*lat*ensemble*leadtime*date) 
lon = ncread(filename, 'X');
lat = ncread(filename, 'Y');
date = ncread(filename, 'S');
leadtimes = ncread(filename, 'L');
ensembels = ncread(filename, 'M');
ind_lon = find(lon>=44 & lon<=64);
ind_lat = find(lat>=24 & lat<=40);
tref = ncread(filename, 'tref', [ind_lon ind_lat 24 4 inf]);

执行代码时失败，报错信息：

Error using horzcat Dimensions of arrays being concatenated are not consistent.

以下是ncdisp()的输出结果：

Source:
           http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NCEP-CFSv2/.HINDCAST/.MONTHLY/tref/1+index/.tref/tref/tref/dods
Format:
           classic
Global Attributes:
           Conventions = 'IRIDL'
Dimensions:
           L = 10
           M = 24
           S = 348
           X = 360
           Y = 181
Variables:
    M   
           Size:       24x1
           Dimensions: M
           Datatype:   single
           Attributes:
                       standard_name = 'realization'
                       long_name     = 'Ensemble Member'
                       pointwidth    = 1
                       gridtype      = 0
                       units         = 'unitless'
    X   
           Size:       360x1
           Dimensions: X
           Datatype:   single
           Attributes:
                       standard_name = 'longitude'
                       pointwidth    = 1
                       gridtype      = 1
                       units         = 'degree_east'
    L   
           Size:       10x1
           Dimensions: L
           Datatype:   single
           Attributes:
                       standard_name = 'forecast_period'
                       long_name     = 'Lead'
                       pointwidth    = 1
                       gridtype      = 0
                       units         = 'months'
    S   
           Size:       348x1
           Dimensions: S
           Datatype:   single
           Attributes:
                       standard_name = 'forecast_reference_time'
                       long_name     = 'Forecast Start Time'
                       pointwidth    = 0
                       calendar      = '360'
                       gridtype      = 0
                       units         = 'months since 1960-01-01'
    Y   
           Size:       181x1
           Dimensions: Y
           Datatype:   single
           Attributes:
                       standard_name = 'latitude'
                       pointwidth    = 1
                       gridtype      = 0
                       units         = 'degree_north'
    tref
           Size:       360x181x24x10x348
           Dimensions: X,Y,M,L,S
           Datatype:   single
           Attributes:
                       pointwidth    = 0
                       standard_name = 'air_temperature'
                       PTVersion     = 2
                       long_name     = 'Reference Temperature'
                       process       = 'Spectral Statistical Interpolation (SSI) analysis from "Final" run.'
                       PDS_TimeRange = 3
                       gribNumBits   = 21
                       center        = 'US Weather Service - National Met. Center'
                       grib_name     = 'TMP'
                       gribcenter    = 7
                       scale_min     = 206.4713
                       gribparam     = 11
                       scale_max     = 315.4577
                       GRIBgridcode  = 3
                       gribleveltype = 105
                       gribfield     = 1
                       units         = 'Kelvin_scale'

问题原因 #

ncread的start参数要求是每个维度的起始索引组成的向量，直接把数组ind_lon、ind_lat和单个数值拼接会导致维度不匹配（数组和标量无法横向拼接）。此外，ncread需要通过start指定起始位置，count指定每个维度读取的数量，不能直接在start里传范围数组。

修正代码 #

filename = 'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NCEP-CFSv2/.HINDCAST/.MONTHLY/tref/1+index/.tref/tref/tref/dods';
ncdisp(filename) 
lon = ncread(filename, 'X');
lat = ncread(filename, 'Y');
date = ncread(filename, 'S');
leadtimes = ncread(filename, 'L');
ensembels = ncread(filename, 'M');

% 筛选经纬度索引
ind_lon = find(lon>=44 & lon<=64);
ind_lat = find(lat>=24 & lat<=40);

% 定义start和count参数
% start: [X起始索引, Y起始索引, M起始索引, L起始索引, S起始索引]
start = [ind_lon(1), ind_lat(1), 24, 4, 1];
% count: [X读取数量, Y读取数量, M读取数量, L读取数量, S读取数量]
count = [length(ind_lon), length(ind_lat), 1, 1, inf];

% 读取指定范围数据
tref = ncread(filename, 'tref', start, count);

说明 #

• start向量的每个元素对应tref的维度顺序X,Y,M,L,S，分别设置各维度的起始位置；
• count向量对应每个维度要读取的数量：经纬度取筛选后的索引长度，集合成员和预报时效各取1个（对应第24个M、第4个L），S维度用inf读取所有起始时间；
• 读取后的tref维度为[length(ind_lon), length(ind_lat), 1, 1, 348]，符合需求。

内容的提问来源于stack exchange，提问作者Behzad