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

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?

Answer 1

[Meg Noah]

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')]});

[Özge]

Thank you for the answer, but my data file is a mat file. I also recorded the TEC data of the global ionospheric maps of the Code with a .mat extension to take up less space. I need to generate a TEC map from these .mat files. My data 73x71x13 in such a 3 dimensional matrix in a .mat file. I want to produce 15 min maps by interpolating from this .mat file.

[Meg Noah]

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);

One of the images in the movie:

I think I'll share this on matlab central!

Here's a link to the movie

TEC Movie

[Meg Noah]

Posted a project for it here:

https://www.mathworks.com/matlabcentral/fileexchange/73858-space-weather-global-tec-displayed-with-geoshow-mapshow

[Özge]

Thank you!!