% TEC mapper script for global TEC maps.
% Written by Tibor Durgonics.
% Modified by Hans-Henrik von Benzon (Now runing on Matlab2015 too).

states = shaperead('landareas.shp', 'UseGeoCoords', true); % Read vector features and attributes from the landareas shapefile.
   
j = 8; % Which epoch should be displayed (1..12) 

figure(j)

[XI,YI]= meshgrid(-180:5:180,87.5:-2.5:-87.5); % Original extent with different Latitude and longitude resolution.
XCol = XI(1,:); XCol = transpose(XCol);
YCol = YI(:,1);
[XI,YI]= meshgrid(-180:2.5:180,87.5:-2.5:-87.5); % The new extent which has now equal reslution in both directions.

TEC=griddata(XCol,YCol,TECMatrixDay(:,:,j),XI,YI,'v4');

worldmap world
meshm(TEC,[0.402,87,-180],[90 180]); 
title(['CODG-',num2str(EPOCH_OF_FIRST_MAP(1,1)),'/',...
num2str(EPOCH_OF_FIRST_MAP(2,1)),'/',num2str(EPOCH_OF_FIRST_MAP(3,1))...
    ,' : ',num2str((j-1)*2),' hour(s) [UT-19 sec]',' (#St:',num2str(N_OF_STATIONS),' #Sat:',num2str(N_OF_SATELLITES),')']);
mapColor=hsv(256);
myMap=mapColor(1:180,:);
colormap(myMap(180:-1:1,:));
colorbar;
h=colorbar;
ylabel(h,'TECU','FontSize', 11);
maxTEC = ceil(max(max(max(TECMatrixDay(:,:,j)))));
caxis([0 maxTEC]);

NrOfStates=length(states);
for ii=1:NrOfStates
    allstates=states(ii);
    geoshow(allstates.Lat,allstates.Lon, 'Color','black','LineWidth',1);
end
