Error using netcdflib NC_EINVALCOORD

3 次查看（过去 30 天）
Steffen Eisele
Steffen Eisele 2020-8-3
Hi,
I'm trying to get some wind direction timelines for different heights. The dataset includes wind data for 37 pressure levels. And I'm getting the following error:
h =
1
Error using netcdflib
The NetCDF library encountered an error during execution of 'getVarsShort' function - 'Index exceeds dimension bound (NC_EINVALCOORDS)'.
Error in netcdf.getVar (line 140)
data = netcdflib(funcstr,ncid,varid,varargin{:});
Error in internal.matlab.imagesci.nc/read (line 635)
data = netcdf.getVar(gid, varid, ...
Error in ncread (line 66)
vardata = ncObj.read(varName, varargin{:});
Error in DirectionTimelines_from_ERAInt6hrly_SJ (line 44)
uwnd_all_temp = ncread(Wu(f).name,'u', [(int16((VolcLong/0.75)+1)) (int16(241-((90+VolcLat)/0.75))) h 1],[1 1 1 (length(ncread(Wu(f).name,
'time')))],[1 1 1 1]); %Some years have a leap year and thus 1464 time periods, but leave at 1460 so doesn't need adapting for leap years.
Error in run (line 91)
evalin('caller', strcat(script, ';'));
The code is the following:
cd ~/Documents/MATLAB/Wind
Wu = dir('u*.nc');
Wv = dir('v*.nc');
YearRange = '2009-2018';
NumYears = 10;
Levels = length(ncread(Wu(1).name, 'level'));
for f = 1:length(Wu)
if (f==1)
date(1:16072) = ncread(Wu(f).name, 'time');
else
date((length(date)+1):((length(date)+1)+length(ncread(Wu(f).name, 'time'))-1)) = ncread(Wu(f).name, 'time');
end
end
Time = length(date);
%Change for each volcano run...
Volcano = '5S'; VolcLat = -5; VolcLong = 95;
mkdir(strcat('../../Output-Wind/',Volcano,'_WindInfo_2019-2018'));
if (VolcLong < 0) %i.e. in Western Hem
VolcLong = 360 + VolcLong;
end
uMatrix = zeros(Time, 37); %Matrix of n rows (1460 time records per year for 6 hourly daily data) by 37 (heights)
%To convert mb or hPa (unit for pressure level data) to altitude, try: www.csgnetwork.com/pressurealtcalc.html
for h = 1:37 %Number of levels
h %So can track progress when running
count = 0;
for f = 1:length(Wu)
uwnd_all_temp = ncread(Wu(f).name,'u', [(int16((VolcLong/0.75)+1)) (int16(241-((90+VolcLat)/0.75))) h 1],[1 1 1 (length(ncread(Wu(f).name, 'time')))],[1 1 1 1]); %Some years have a leap year and thus 1464 time periods, but leave at 1460 so doesn't need adapting for leap years.
uwnd_all_temp = squeeze(uwnd_all_temp);
uMatrix((count+1):(count+length(uwnd_all_temp)),h) = uwnd_all_temp;
count = count+length(uwnd_all_temp);
end
end
vMatrix = zeros(Time, 37); %Matrix of rows (no. years * no. days * no. records/day) by 37 cols (heights)
for h = 1:37 %Number of levels
h
count = 0;
for f = 1:length(Wv)
vwnd_all_temp = ncread(Wv(f).name,'v', [(int16((VolcLong/0.75)+1)) (int16(241-((90+VolcLat)/0.75))) h 1],[1 1 1 (length(ncread(Wu(f).name, 'time')))],[1 1 1 1]); %Some years have a leap year and thus 1464 time periods, but leave at 1460 so doesn't need adapting for leap years.
vwnd_all_temp = squeeze(vwnd_all_temp);
vMatrix((count+1):(count+length(vwnd_all_temp)),h) = vwnd_all_temp;
count = count+length(vwnd_all_temp);
end
end
Dmatrix = zeros(Time,Levels); % matrix of direction (rows) at each height (cols)
Imatrix = zeros(Time,Levels); % matrix of intensity (rows) at each height (cols)
AltitudeLabels = {'32.5 km' '31 km' '29.5 km' '28 km' '27 km' '26 km' '23.5 km' '21.5 km' '19.5 km' '17.5 km' '16 km' '14.5 km' '13.5 km' '12.5 km' '12 km' '11 km' '10.5 km' '9 km' '8 km' '7 km' '6.5 km' '5.5 km' '5 km' '4.2 km' '3.6 km' '3 km' '2.5 km' '2.2 km' '1.9 km' '1.7 km' '1.5 km' '1.2 km' '1 km' '800 m' '500 m' '300 m' '100 m'};
%To calculate and save direction and intensity per record with height:
for h = 1:Levels
Dmatrix(:, h) = mod(90-atan2d(vMatrix(:,h),uMatrix(:,h)), 360); %calculation from u and vmatrix... for Direction.
Imatrix(:, h) = sqrt((uMatrix(:,h).^2)+(vMatrix(:,h).^2)); %calculation from u and vmatrix... for Speed.
end
save(strcat('../../Output-Wind-Info/',Volcano,'_WindInfo_2009-2018/Dir-Int_Matrix.mat'),'Dmatrix','Imatrix');
%From: http://www.csgnetwork.com/pressurealtcalc.html
Altitude = {32435,30760,29675,28180,27115,25905,23315,21630,19315,17660,15790,14555,13505,12585,11770,11030,10360,9160,8115,7180,6340,5570,4865,4205,3590,3010,2465,2205,1950,1700,1455,1220,990,760,540,325,110};
AltAsNumber = cell2mat(Altitude);
AltAsNumber = reshape(AltAsNumber,37,1);
Has anyone an idea how to solve the problem?
Thanks in advance
Steffen

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Signal Analysis 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by