% SAMPLE_SCRIPT

% This sample script was provided by Neil Gordon with some modification by
% Joel Norris.  It reads in and plots meteorological data from 1993 March 14
% 00Z for a region off the east coast of the U.S. (20-50N, 270-310E).

% The parameters we will be reading in are Z1000, Z500, T1000, u1000, v1000,
% u500, and v500.

% Note that the header line ("yr  mn  dy  hr  lat  lon  temperature...") must
% be deleted from the data files before the data can be read in correctly

% Start up matlab by typing "matlab" in the directory where the script and the
% data files are located.

% "matlab -nojvm" will run matlab without the fancy popup menu.

% If a semicolon is not put after a command, then Matlab will print the results
% to the screen.  You can print to the screen the values of any variable simply
% by typing the name of the variable.

% The commands in this script can be entered interactively into matlab using
% copy and paste.

% The entire script can run at once by typing "sample_script" in the matlab
% window.

% Now let's begin.

% First we specify the file we want to read Z1000 data from.
infile = 'Z1000.1993031400.1x1.txt';

% Now we read in the data and store it in matrix X.
X = load(infile);

% Each column of matrix X corresponds to a different variable.
year_raw = X(:,1);
month_raw = X(:,2);
day_raw = X(:,3);
hour_raw = X(:,4);
latitude_raw = X(:,5);
longitude_raw = X(:,6);
Z1000_raw = X(:,7);

% An alternative method is to read in all data using a single command.
% We designate the type of numbers to be read in for each column using
% '%d%d%d%d%f%f%f', where %d means an integer number and %f means a floating
% point number (4 integers for yr mn dy hr followed by 3 floating points for
% lat lon data).
% This command will be too long for a single line, so we use "..." to indicate
% the command extends to the following line
[year_raw,month_raw,day_raw,hour_raw,latitude_raw,longitude_raw,Z1000_raw] =...
  textread('Z1000.1993031400.1x1.txt','%d%d%d%d%f%f%f');

% Now we read in the rest of the parameters.  We only need to extract the last
% column from matrix X since the other variables will be the same.
infile = 'T1000.1993031400.1x1.txt';
X = load(infile);
T1000_raw = X(:,7);

infile = 'u1000.1993031400.1x1.txt';
X = load(infile);
u1000_raw = X(:,7);

infile = 'v1000.1993031400.1x1.txt';
X = load(infile);
v1000_raw = X(:,7);

infile = 'Z500.1993031400.1x1.txt';
X = load(infile);
Z500_raw = X(:,7);

infile = 'u500.1993031400.1x1.txt';
X = load(infile);
u500_raw = X(:,7);

infile = 'v500.1993031400.1x1.txt';
X = load(infile);
v500_raw = X(:,7);

% Now lets check to see what variables we have in matlab (not necessary for the
% plot).
whos

% The raw data are loaded in as one long series of numbers, so we need to
% reshape the 1-dimensional arrays into 2-dimensional arrays of 360 longitude
% by 181 latitude points
Z1000_glob = reshape(Z1000_raw,[360 181]);
T1000_glob = reshape(T1000_raw,[360 181]);
u1000_glob = reshape(u1000_raw,[360 181]);
v1000_glob = reshape(v1000_raw,[360 181]);
Z500_glob = reshape(Z500_raw,[360 181]);
u500_glob = reshape(u500_raw,[360 181]);
v500_glob = reshape(v500_raw,[360 181]);
longitude_glob = reshape(longitude_raw,[360 181]);
latitude_glob = reshape(latitude_raw,[360 181]);

% It will be helpful to have longitude in a single 360 array and latitude in a
% single 181 array rather than 360x181 2-d arrays.
lon1d = longitude_glob(:,1);
lat1d = latitude_glob(1,:);

% To avoid unneeded processing of global data, let's extract data for the
% region of interest (20-50N, 270-310E).
lon1d_reg = find(lon1d>=270 & lon1d<=310);
lat1d_reg = find(lat1d>=20 & lat1d<=50);
Z1000_reg = Z1000_glob(lon1d_reg,lat1d_reg);
T1000_reg = T1000_glob(lon1d_reg,lat1d_reg);
u1000_reg = u1000_glob(lon1d_reg,lat1d_reg);
v1000_reg = v1000_glob(lon1d_reg,lat1d_reg);
Z500_reg = Z500_glob(lon1d_reg,lat1d_reg);
u500_reg = u500_glob(lon1d_reg,lat1d_reg);
v500_reg = v500_glob(lon1d_reg,lat1d_reg);
longitude_reg = longitude_glob(lon1d_reg,lat1d_reg);
latitude_reg = latitude_glob(lon1d_reg,lat1d_reg);

% Let's load in coastline location data so we can plot it on the map.
load coast;

% Now let's make a contour plot.  Use "clf" to clear the figure.
clf 

% Plot the map.
plot(long,lat,'k');

% Use "hold on" to add one plot to another.
hold on

% Since the coastline data are -180 to 180 but our data are 0 to 360, we will
% add another map to right of the first map to get a 0-360 region with
% coastlines (Matlab is lame this way, but a more sophisticated user could 
% probably do better).
plot(long+360,lat,'k');

% Restrict the domain to 20-50N, 270-310E.
axis([270 310 20 50]);

% Plot contours of Z1000 over the map.
hold on
contour(longitude_reg,latitude_reg,Z1000_reg);

% Add a title.
title('Z1000 (gpm) [1993 March 14 00Z]');

% Let's redo the plot.  This time we will use only black contours, where 'k'
% indicates black.  For future reference, 'r', 'b', 'g', 'm', 'c', and 'y'
% correspond to red, blue, green, magenta, cyan, and yellow.
clf
plot(long,lat,'k');
hold on
plot(long+360,lat,'k');
axis([270 310 20 50]);
hold on
contour(longitude_reg,latitude_reg,Z1000_reg,'k');
title('Z1000 (gpm) [1993 March 14 00Z]');

% Let's redo the plot again to specify that contours should be drawn between 
% -400 and +400 at intervals of 40.  The sequence [A:B:C] means go from A to C
% at intervals of B.
clf
plot(long,lat,'k');
hold on
plot(long+360,lat,'k');
axis([270 310 20 50]);
hold on
contour(longitude_reg,latitude_reg,Z1000_reg,[-400:40:400],'k');
title('Z1000 (gpm) [1993 March 14 00Z]');

% Let's redo the plot again to label the contours.  Note the changes to the
% line with the "contour" commmand and the new "clabel" command.
clf
plot(long,lat,'k');
hold on
plot(long+360,lat,'k');
axis([270 310 20 50]);
hold on
[c,h] = contour(longitude_reg,latitude_reg,Z1000_reg,[-400:40:400],'k');
clabel(c,h,'color','k');
title('Z1000 (gpm) [1993 March 14 00Z]');

% Let's now add red dashed contours of T1000 to the plot with an interval of 4.
% The two hyphens in 'r--' tell Matlab to make dashes.  It may be difficult to
% see the dashes in the plot on the screen, but they will come out better when
% printed.
hold on
[c,h] = contour(longitude_reg,latitude_reg,T1000_reg,[-40:4:40],'r--');
clabel(c,h,'color','r');

% We should also change the title to include the new parameter.
title('Z1000 (gpm) and T1000 (C) [1993 March 14 00Z]');

% Let's redo the plot, this time without the labels.
clf
plot(long,lat,'k');
hold on
plot(long+360,lat,'k');
axis([270 310 20 50]);
hold on
contour(longitude_reg,latitude_reg,Z1000_reg,[-400:40:400],'k');
hold on
contour(longitude_reg,latitude_reg,T1000_reg,[-40:4:40],'r--');
title('Z1000 and T1000 [1993 March 14 00Z]');
                                                                                
% Now add 1000-hPa wind vectors colored blue and update title.
hold on
quiver(longitude_reg,latitude_reg,u1000_reg,v1000_reg,'b');
title('Z1000, T1000, and wind1000 [1993 March 14 00Z]');

% The plot looks busy because the vectors are too closely spaced.  Let's redo
% the whole plot, but this time only plot every other vector.  Note how we use
% the [A:B:C] notation to plot from 1 to 41 longitude points and 1 to 31 
% latitude points at intervals of 2.  There are 41 longitude points in 270-310E
% and 31 latitude points in 20-50N (inclusive).
clf
plot(long,lat,'k');
hold on
plot(long+360,lat,'k');
axis([270 310 20 50]);
hold on
contour(longitude_reg,latitude_reg,Z1000_reg,[-400:40:400],'k');
hold on
contour(longitude_reg,latitude_reg,T1000_reg,[-40:4:40],'r--');
hold on
quiver(longitude_reg([1:2:41],[1:2:31]),latitude_reg([1:2:41],[1:2:31]),...
  u1000_reg([1:2:41],[1:2:31]),v1000_reg([1:2:41],[1:2:31]),'b');
title('Z1000, T1000, and wind1000 [1993 March 14 00Z]');

% Let's add 500-hPa wind plotted in green and update the title.
hold on
quiver(longitude_reg([1:2:41],[1:2:31]),latitude_reg([1:2:41],[1:2:31]),...
  u500_reg([1:2:41],[1:2:31]),v500_reg([1:2:41],[1:2:31]),'g');
title('Z1000, T1000, wind1000, and wind500 [1993 March 14 00Z]');

% One disadvantage of the automatic scaling that Matlab applies to the vectors
% is that it is difficult to compare the magnitude of wind at 1000-hPa with
% wind at 500-hPa because they are scaled differently.  Putting a 0 near the
% end of the quiver command turns off the scaling.  Let's redo the plot using
% that.
clf
plot(long,lat,'k');
hold on
plot(long+360,lat,'k');
axis([270 310 20 50]);
hold on
contour(longitude_reg,latitude_reg,Z1000_reg,[-400:40:400],'k');
hold on
contour(longitude_reg,latitude_reg,T1000_reg,[-40:4:40],'r--');
hold on
quiver(longitude_reg([1:2:41],[1:2:31]),latitude_reg([1:2:41],[1:2:31]),...
  u1000_reg([1:2:41],[1:2:31]),v1000_reg([1:2:41],[1:2:31]),0,'b');
title('Z1000, T1000, and wind1000 [1993 March 14 00Z]');
hold on
quiver(longitude_reg([1:2:41],[1:2:31]),latitude_reg([1:2:41],[1:2:31]),...
  u500_reg([1:2:41],[1:2:31]),v500_reg([1:2:41],[1:2:31]),0,'g');
title('Z1000, T1000, wind1000, and wind500 [1993 March 14 00Z]');

% Turning off the scaling resulted in really long wind vectors.  We will need 
% to apply our own scale factor to the winds to reduce them (unlike Matlab, we
% will apply the same factor for all winds).  We will need to use trial and
% error to see what scale factor is best (0.03 seems to be best for this
% figure).  This is lame, but perhaps a more sophisticated user could do
% better.  At least we can see that the 500-hPa winds are clearly faster than
% those at 1000-hPa.
scale_factor = 0.03;
u1000_scale = u1000_reg * scale_factor;
v1000_scale = v1000_reg * scale_factor;
u500_scale = u500_reg * scale_factor;
v500_scale = v500_reg * scale_factor;
clf
plot(long,lat,'k');
hold on
plot(long+360,lat,'k');
axis([270 310 20 50]);
hold on
contour(longitude_reg,latitude_reg,Z1000_reg,[-400:40:400],'k');
hold on
contour(longitude_reg,latitude_reg,T1000_reg,[-40:4:40],'r--');
hold on
quiver(longitude_reg([1:2:41],[1:2:31]),latitude_reg([1:2:41],[1:2:31]),...
  u1000_scale([1:2:41],[1:2:31]),v1000_scale([1:2:41],[1:2:31]),0,'b');
title('Z1000, T1000, and wind1000 [1993 March 14 00Z]');
hold on
quiver(longitude_reg([1:2:41],[1:2:31]),latitude_reg([1:2:41],[1:2:31]),...
  u500_scale([1:2:41],[1:2:31]),v500_scale([1:2:41],[1:2:31]),0,'g');
title('Z1000, T1000, wind1000, and wind500 [1993 March 14 00Z]');

% Now that the plot is finished, let's print it to a postscript file named
% "sample_exercise1.ps" so we can print out a hardcopy version. "-dpsc" means
% print a ps file in color and "-dps" means print a ps file in black and white.
print -dpsc sample_exercise1.ps

% Instead of contour plots, let's now try a simpler x-y plot.  For example, we
% could plot how Z1000 and Z500 vary with longitude between 270-310E along the
% constant latitude circle of 40N.  To begin, we will need to extract the data
% along 40N.  We'll use the "find" command again.
lat_40 = find(lat1d==40);
longitude_40 = longitude_glob(lon1d_reg,lat_40);
Z1000_40 = Z1000_glob(lon1d_reg,lat_40);
Z500_40 = Z500_glob(lon1d_reg,lat_40);

% Now do the plot, with Z1000 in black and Z500 in green.
clf
plot(longitude_40,Z1000_40,'k');
hold on
plot(longitude_40,Z500_40,'g');

% Let's label the axes, provide a title, and print to a file.
xlabel('Longitude');
ylabel('Geopotential Height (gpm)');
title('Z1000 and Z500');
print -dpsc sample_exercise2.ps

% Let's do the plot again, this time changing the Y-axis using the "axis"
% command so that it goes to the approximate height of the tropopause.
clf
plot(longitude_40,Z1000_40,'k');
hold on
plot(longitude_40,Z500_40,'g');
axis([270 310 -1000 12000]);
xlabel('Longitude');
ylabel('Geopotential Height (gpm)');
title('Z1000 and Z500');

% In addition to plotting meteorological data, we will need to perform various
% calculations and manipulations.  Let's try repeatedly applying 1-2-1
% smoothing to the Z1000 field (actually, it will be 1-4-1 smoothing since we
% are working in two dimensions).  The equation is below.
%
% Xsmooth(i,j) = (4*X(i,j) + X(i+1,j) + X(i-1,j) + X(i,j+1) + X(i,j-1)) / 8
%
% We will need to go back to the global data to obtain the edge values needed
% for smoothing that are outside the region of interest.  Note how the array
% indices are arranged to carry this out (A:B means go from A to B using an
% interval of 1.
Zsmooth_glob = Z1000_glob;
Zsmooth_glob(2:358,2:180) = (4*Zsmooth_glob(2:358,2:180) +...
  Zsmooth_glob(1:357,2:180) + Zsmooth_glob(3:359,2:180) +...
  Zsmooth_glob(2:358,1:179) + Zsmooth_glob(2:358,3:181)) / 8;

% The above step will only provide light smoothing, so let's repeat it 100 more
% times to obtain something dramatic.
for n=1:100
  Zsmooth_glob(2:358,2:180) = (4*Zsmooth_glob(2:358,2:180) +...
    Zsmooth_glob(1:357,2:180) + Zsmooth_glob(3:359,2:180) +...
    Zsmooth_glob(2:358,1:179) + Zsmooth_glob(2:358,3:181)) / 8;
end

% Now extract the smoothed values in the region of interest.
Zsmooth_reg = Zsmooth_glob(lon1d_reg,lat1d_reg);

% Let's compare smoothed Z1000 (green) to original Z1000 (black) in a plot.
% Note how smoothing has partially filled in the low center and reduced the
% height of the high center.
clf
plot(long,lat,'k');
hold on
plot(long+360,lat,'k');
axis([270 310 20 50]);
hold on
[c,h] = contour(longitude_reg,latitude_reg,Z1000_reg,[-400:40:400],'k');
clabel(c,h,'color','k');
hold on
[c,h] = contour(longitude_reg,latitude_reg,Zsmooth_reg,[-400:40:400],'g');
clabel(c,h,'color','g');
title('Smoothed and Original Z1000 (gpm) [1993 March 14 00Z]');

% As another exercise, let's plot the variation of the Coriolis parameter in
% the region of interest.  First create a conversion factor from degrees to
% radians (and the value of Omega).
RAD = 3.141592654 / 180;
OMEGA = 7.292e-5;

% Now calculate the Coriolis parameter and plot it.
f_reg = 2 * OMEGA * sin(latitude_reg*RAD);
clf
plot(long,lat,'k');
hold on
plot(long+360,lat,'k');
axis([270 310 20 50]);
hold on
[c,h] = contour(longitude_reg,latitude_reg,f_reg,'k');
clabel(c,h,'color','k');
title('Coriolis Parameter (s-1) [1993 March 14 00Z]');