How can I generate a TEC map by interpolating a 3D matrix in a .mat file?

Question

Özge 2020-1-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/499008-how-can-i-generate-a-tec-map-by-interpolating-a-3d-matrix-in-a-mat-file

评论： Özge 2020-1-9

IGS-2-hour-TEC-map-taken-from-ftp-cdawebgsfcnasagov-pub-data-gps-tec2hr-igs.png

I have hourly data and a TECU value for each grid point in the matte file. And I want to produce a TEC map using this data.

long=-180:5:180;

lat=87.5:-2.5:-87.5;

[long2, lat2]=meshgrid(long, lat);

My data 73x71x13 in such a 3 dimensional matrix. I got the error how I did the interpolation. How can i produce a TEC map like in the image?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Meg Noah 2020-1-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/499008-how-can-i-generate-a-tec-map-by-interpolating-a-3d-matrix-in-a-mat-file#answer_408784

在 MATLAB Online 中打开

clc
close all
clear all
filename = 'gps_tec2hr_igs_20150317_v01.cdf';
[data, info] = cdfread(filename);
latIGS = cdfread(filename,'Variable','lat'); latIGS = latIGS{1};
lonIGS = cdfread(filename,'Variable','lon'); lonIGS = lonIGS{1};
igs = cdfread(filename,'Variable','tecIGS'); 
Epoch = cdfread(filename,'Variable','Epoch'); 
figure(1)
for iplot = 1:12
subplot(3,4,iplot)
imagesc(lonIGS,latIGS,igs{iplot},[0 104]);
colormap(jet);
colorbar
title([datestr(todatenum(Epoch{iplot}),'dd-mmm-yyyy HH:MM:SS')]);
set(gca,'ydir','normal');
end
figure(2)
iplot= 10; % get the data for 18:00
imagesc(lonIGS,latIGS,igs{iplot});
colormap(jet);
colorbar
title({'IGS 2h TEC Map';'TEC/TECU'; ...
    [datestr(todatenum(Epoch{iplot}),'dd-mmm-yyyy HH:MM:SS')]});
set(gca,'ydir','normal');
axis equal
axis tight
% georeference the data
latIGS = double(latIGS);
lonIGS = double(lonIGS);
R1 = makerefmat('RasterSize',size(igs{iplot}), ...
    'Latlim', [min(latIGS(:)) max(latIGS(:))], ...
    'Lonlim', [min(lonIGS(:)) max(lonIGS(:))] );
dataTEC = flipud(igs{iplot});
figure(3);
hold on;
mapshow(zeros(size(dataTEC)),R1,'CData',dataTEC,'DisplayType','surface');
cW = mapshow(dataTEC,R1,'DisplayType','contour','linecolor','black', ...
    'showtext','on');
colormap(jet)
colorbar('location','southoutside');
axis equal
axis tight
title({'IGS 2h TEC Map';'TEC/TECU'; ...
    [datestr(todatenum(Epoch{iplot}),'dd-mmm-yyyy HH:MM:SS')]});
figure(4);
gcm = worldmap([min(latIGS) max(latIGS)], [min(lonIGS) max(lonIGS)]);
hold on;
load('coast');
geoshow(gcm,dataTEC,R1,'DisplayType','texturemap');
plotm(lat,long,'color','black');
colormap(jet)
colorbar('location','southoutside');
axis equal
axis tight
title({'IGS 2h TEC Map';'TEC/TECU'; ...
    [datestr(todatenum(Epoch{iplot}),'dd-mmm-yyyy HH:MM:SS')]});

4 个评论
显示 2更早的评论隐藏 2更早的评论

Meg Noah 2020-1-7

在 MATLAB Online 中打开

Here's a solution for the 1-hour data I found online. You can modify it for your data, but I don't know what hours you have.

clc
close all
clear all
% https://cdaweb.gsfc.nasa.gov/pub/data/gps/tec1hr_igs/2015/
filename = 'gps_tec1hr_igs_20150317_v01.cdf';
[data, info] = cdfread(filename);
latIGS = cdfread(filename,'Variable','lat'); latIGS = latIGS{1};
lonIGS = cdfread(filename,'Variable','lon'); lonIGS = lonIGS{1};
tecUHR = cdfread(filename,'Variable','tecUHR');
Epoch = cdfread(filename,'Variable','Epoch');
%% 4 ways to visualize a tecUHR
%  in this case - average integrated over 1 hour
figure(1)
for iplot = 1:24
    subplot(6,4,iplot)
    imagesc(lonIGS,latIGS,tecUHR{iplot},[0 104]);
    colormap(jet);
    colorbar
    title([datestr(todatenum(Epoch{iplot}),'dd-mmm-yyyy HH:MM:SS')]);
    set(gca,'ydir','normal');
end
iplot= 19; % get the data for 18:00
tecUHR_19UTC = tecUHR{iplot};
figure(2)
imagesc(lonIGS,latIGS,tecUHR_19UTC);
colormap(jet);
colorbar
strTitle = {'IGS 1h TEC Map';'TEC UHR'; ...
    [datestr(todatenum(Epoch{iplot}),'dd-mmm-yyyy HH:MM:SS')]};
title(strTitle);
set(gca,'ydir','normal');
axis equal
axis tight
% georeference the data
latIGS = double(latIGS);
lonIGS = double(lonIGS);
R1 = makerefmat('RasterSize',size(tecUHR_19UTC), ...
    'Latlim', [min(latIGS(:)) max(latIGS(:))], ...
    'Lonlim', [min(lonIGS(:)) max(lonIGS(:))], 'ColumnsStartFrom','North' );
[ny,nx] = size(tecUHR_19UTC);
dataTECUHR = tecUHR_19UTC; %flipud(tecUHR_19UTC);
load('coast');
figure(3);
hold on;
mapshow(zeros(size(dataTECUHR)),R1,'CData',dataTECUHR,'DisplayType','surface');
cW = mapshow(dataTECUHR,R1,'DisplayType','contour','linecolor','black', ...
    'showtext','on');
colormap(jet)
colorbar('location','southoutside');
axis equal
axis tight
title(strTitle);
mapshow(long,lat,'color','white');
figure(4);
gcm = worldmap([min(latIGS) max(latIGS)], [min(lonIGS) max(lonIGS)]);
hold on;
geoshow(gcm,dataTECUHR,R1,'DisplayType','texturemap');
plotm(lat,long,'color','black');
colormap(jet)
colorbar('location','southoutside');
axis equal
axis tight
title(strTitle);
%% interpolating to a new timegrid
% make into a block of one hour intervals
tecUHR3D = zeros(71,73,24);
for icell = 1:24
    tecUHR3D(:,:,icell) = tecUHR{icell};
end
fig = figure();
% create a placeholder image
gcm = worldmap([min(latIGS) max(latIGS)], [min(lonIGS) max(lonIGS)]);
setm(gcm, 'MlabelParallel', 'south');
hold on;
tecUHR_15min = squeeze(tecUHR3D(:,:,1));
geoshow(gcm,tecUHR_15min,R1,'DisplayType','texturemap');
plotm(lat,long,'color','black');
colormap(jet)
caxis([0 120]);
colorbar('location','eastoutside');
axis equal
axis tight
UTCHR = 0; alpha = 0;
title([datestr(todatenum(Epoch{1})+(UTCHR+alpha)/24,'dd-mmm-yyyy HH:MM:SS')]);
v = VideoWriter('TECUHR_20150317.mp4', 'MPEG-4'); % New
v.Quality = 95;
open(v);
% the most basic interpolation is linear between two points
% loop through hours 1 to 23
for UTCHR = 1:22
    
    iUTCHR = UTCHR + 1; % index into the NLATxNLONGxNHOUR array
    
    for itime = 0:15:45
        
        alpha = itime/60;  % first quarter of hour  = 0.25, etc.
        tecUHR_15min = squeeze((1-alpha).*tecUHR3D(:,:,iUTCHR) + alpha.*(tecUHR3D(:,:,iUTCHR+1)));
        geoshow(gcm,tecUHR_15min,R1,'DisplayType','texturemap');
        plotm(lat,long,'color','black');
        title([datestr(todatenum(Epoch{1})+(UTCHR+alpha)/24,'dd-mmm-yyyy HH:MM:SS')]);
        frame = getframe(gcf);
        writeVideo(v,frame);
        
    end
    
end
% very last frame
UTCHR = 23;
iUTCHR = UTCHR + 1; % index into the NLATxNLONGxNHOUR array
itime = 0;
alpha = itime/60;  % first quarter of hour  = 0.25, etc.
tecUHR_15min = squeeze((tecUHR3D(:,:,iUTCHR)));
geoshow(gcm,tecUHR_15min,R1,'DisplayType','texturemap');
plotm(lat,long,'color','black');
title([datestr(todatenum(Epoch{1})+(UTCHR+alpha)/24,'dd-mmm-yyyy HH:MM:SS')]);
frame = getframe(gcf);
writeVideo(v,frame);
% close the video
close(v);