[gps-developers] Re: Matlab plotting of gps data

William Rideout brideout at haystack.mit.edu
Wed Sep 3 14:11:38 EDT 2003


John and Anthea,

I've installed the script plot_gps_cart on hyperion.  Just type

help plot_gps_cart

from Matlab on hyperion to get a description of how to use it.

Also attached are the two *.m files, which you can use on your local 
version of Matlab.  These files are in CVS under 
Millstone/gps/source/scripts.

I'll let you know when I incorporate some of your suggestions.

Bill


John Foster wrote:
> Neat  
> 
> I want a copy of this script!
> 
> My version ('gpsort') is probably the grandparent of what progressed
> to Phil/Anthea/Bill and still takes ~120 sec (on my laptop) to process
> 20 min of data.
> 
> Not great for making my own movies!!!
> 
> JF
> 
> Comments:
> 
> having lots of user-accessible switches to set the plot characteristics
> is VG.
> 
> Having an polar-projection version would be VG, as well
> 
> my options include overlaying a magnetic coord contour grid (developed from
> your 2-D model output - usualy run for 350 km altitude) - useful at times
> 
> I also like adding F-region the shadow line (SZA=100) which can be
> obtained from an SZA grid for the date/time of the frame (more complexities)
> 
> Final data-manipulation suggestion:
> 
> being able to fill in some of the missing bins by interpolation is very
> helpful when preparing final/publication versions of the maps - alternately,
> we could set up a few hundred more GPS sites or launch a bunch more GPS
> satellites.
> 
> _______________________________________________
> gps-developers mailing list
> gps-developers at openmadrigal.org
> http://www.openmadrigal.org/mailman/listinfo/gps-developers
> 
> 


-- 
Bill Rideout
MIT Haystack Observatory
Email: brideout at haystack.mit.edu
Phone: 781 981-5624
-------------- next part --------------
function [] = plot_gps_cart(directory, ...
              start_time, end_time, ...
              step_lat, start_lat, end_lat, ...
              step_lon, start_lon, end_lon, ...
              min_tec, max_tec, ...
              filename, whiteBackground)
% plot_gps_cart plots gps data in cartesian coordinates
% 
%     directory - path to directory with tec* files
%     start_time - Matlab date string or date number to set starting time
%     end_time - Matlab date string or date number to set ending time
%     step_lat - degrees latitude per bin
%     start_lat - lower limit of latitude to show in map (-90 to 90)
%     end_lat - upper limit of latitude to show in map (-90 to 90)
%     step_lon - degrees longitude per bin
%     start_lon - lower limit of longitude to show in map (-180 to 180)
%     end_lon - upper limit of longitude to show in map (-180 to 180)
%     min_tec - lower limit of TEC color scale
%     max_tec - upper limit of TEC color scale
%     filename - filename to save image.  If blank, just display, and let user save
%     whiteBackground - 1 for white background, 0 for black
%
%   Notes:
%    Uses the median of the binned data, not the mean.
%    The bins are not overlapping, so that each data point shows up in one and only one bin.
%
%   Example:
%    plot_gps_cart('/usr/algo/DATA_04_2002',731321.0,731321.0139,3,0,90,3,-90,0,0,100,'p1.jpeg',1) 
%
%   Written by B. Rideout 9/2/2003

% Check for input mistakes
if (nargin ~= 13)
    error('Not enough arguments')
end

% cd to directory
presentDir = cd;
cd(directory);
filelist = dir('tec-*-*-*.dat');
fileSize = size(filelist);
numFiles = fileSize(1,1);


% throw an error if not files found
if (numFiles == 0)
    cd(presentDir)
    error('No tec-???-??-??.dat files found in directory')
end


% convert start_time and end time into Matlab Datenum's
startTime = datenum(start_time);
endTime = datenum(end_time);

% make sure both times have the same year, since files don't identify year
startvec = datevec(startTime);
startyear = startvec(1,1);
endvec = datevec(endTime);
endyear = endvec(1,1);
if (startyear ~= endyear);
    cd(presentDir)
    error('Plot cannot cross year boundry, since files do not include year');
end

% load all files within those times
count = 0;
for i=1:numFiles
  thisfile = filelist(i).name;
  % skip files of wrong length
  if (length(thisfile) == 17)
    % compute its time
    day = double(str2num(thisfile(5:7)));
    hour = double(str2num(thisfile(9:10)));
    minute = double(str2num(thisfile(12:13)));
    thistime = datenum(startyear,0,0,0,0,0) + day + hour/24.0 + minute/1440.0;
    % skip files outside the time range
    if (thistime >= startTime & thistime <= endTime)
      fid = fopen(thisfile);
      if (count == 0)
          this_data = fscanf(fid, '%g', [4 inf]).';
          count = 1;
      else
          A = fscanf(fid, '%g', [4 inf]).';
          this_data = [this_data; A];
      end
      fclose(fid);
    end
  end
end

% if no files in given time range, throw error
if (count == 0)
    error('No files found in given time range');
end

% convert long to -180 to 180
for row=1:length(this_data)
    if (this_data(row,3) > 180.0)
        this_data(row,3) = this_data(row,3) - 360.0;
    end
end

% bin the data
[bin_dat] = gps_bin(step_lat, step_lon, this_data);

% create lon_grid and lat_grid
lat_grid = -90:step_lat:90;
lon_grid = -180:step_lon:180;

% find start and end indexes of grids
startLatIndex = 1;
endLatIndex = 1;
for i=1:length(lat_grid)
    if (start_lat >= lat_grid(i))
        startLatIndex = i;
    end
    if (end_lat >= lat_grid(i))
        endLatIndex = i;
    end
end

startLonIndex = 1;
endLonIndex = 1;
for i=1:length(lon_grid)
    if (start_lon >= lon_grid(i))
        startLonIndex = i;
    end
    if (end_lon >= lon_grid(i))
        endLonIndex = i;
    end
end

% shrink the grids and data to desired range

bin_dat = bin_dat(startLatIndex:endLatIndex , startLonIndex:endLonIndex);
lat_grid = lat_grid(startLatIndex:endLatIndex);
lon_grid = lon_grid(startLonIndex:endLonIndex);



% plot data
clf;
axis([min(lon_grid) max(lon_grid) min(lat_grid) max(lat_grid)]);
set(gcf,'position',[232 258 560 420]);
pcolor(lon_grid, lat_grid, bin_dat);
caxis([min_tec, max_tec]);
cb = colorbar;
set(get(cb,'Title'), 'String', 'TEC');

title(sprintf('GPS TEC Map  from %s to %s', ...
      datestr(startTime), datestr(endTime)));
xlabel('Geodetic Longitude, Deg');
ylabel('Geodetic Latitude, Deg');
shading flat;
  
if (whiteBackground == 0)
  % turn background black
  set(gca,'Color', [0 0 0]);
end
  
% add coastline
load coast.mat;
hold on;
if (whiteBackground == 0)
  plot(long, lat, 'w-');
else
  plot(long, lat, 'k-');
end

% save if filename given
if (length(filename) > 1)
    eval(['print -djpeg -zbuffer ' filename]);
end
  
  
% return to original dir
cd(presentDir)
-------------- next part --------------
function [bin_dat] = ...
    gps_bin(lat_step, lon_step, gdata)
% [bin_dat] = ...
%       GPS_BIN(lat_step, lon_step, gdata)
%
% Bin-sort GPS TEC files (from Anthea Coster) and return a bin-sorted array.
%
%  Modified by B. Rideout 8/29/2003 to improve speed, return median rather
%     than mean of binned data.
%
% INPUTS:
%   lat_step    step size of latitude bins, in degrees.
%   lon_step    step size longitude bins, in degrees.
%   gdata       N x 4 array containing TEC values, in the format:
%               time gdlat gdlon tec
%               ...
%               time = decimal day of the year (with fraction)
%               gdlat = geodetic latitude, deg
%               gdlon = geodetic longitude, deg
%               tec   = TEC, std units
%
% RETURNS:
%   bin_dat     K x L array [K = 180/lat_step = 1  L = 360/lon_step + 1]
%               with binned TEC data.  Bins with no data are set to NaN.
%

% To make things faster, set maximum size of bin
maxBinSize = 50;

% sp_bin_dat is 3D array holding binned data 
sp_bin_dat = zeros(ceil(180.0/lat_step)+1, ceil(360.0/lon_step)+1, maxBinSize);

% bin_dat is 2D array holding mean of binned data - to be returned
bin_dat = zeros(ceil(180.0/lat_step)+1, ceil(360.0/lon_step)+1);

% bin_index is 2D array holding index of each binned data cell
bin_index = zeros(ceil(180.0/lat_step)+1, ceil(360.0/lon_step)+1);


% loop through each data point
for i=1:length(gdata)
    % get index into sp_bin_dat
    indexLat = ceil(((gdata(i,2) + 90.0)+(lat_step/2.0))/lat_step);
    indexLon = ceil(((gdata(i,3) + 180.0)+(lon_step/2.0))/lon_step);
    % check whether that bin is full
    newIndex = bin_index(indexLat, indexLon);
    if (newIndex < maxBinSize)
        % append this point
	sp_bin_dat(indexLat, indexLon, newIndex + 1) = gdata(i,4);
	bin_index(indexLat, indexLon) = newIndex + 1;
    end
end

% now that data is loaded in sp_bin_dat, populate bin_dat
for i=1:ceil(180.0/lat_step)+1
    for j=1:ceil(360.0/lon_step)+1
        if (bin_index(i,j) > 0)
	   bin_dat(i,j) = median(sp_bin_dat(i,j,:));
	else
	   bin_dat(i,j) = NaN;
	end
    end
end



More information about the gps-developers mailing list