Data analysis with client tools: Matlab¶
A. Sea level from tide gauges¶
There is an overview at http://www.soest.hawaii.edu/coasts/sealevel/. The sea level from tide gauges can be found here at http://uhslc.soest.hawaii.edu/data/. At this link you can choose the tide-gauge of preference and download as .dat, .csv and .netcdf, alternately as hourly or daily files. Note the file naming convention is a “d” for daily and “h” for hourly, then the station ID.
The .csv data are most easily used, since these can be edited (if needed) with a text editor and can be read directly into Matlab using read. For example,
>> load d001.csv
Will create a variable called d001 with four columns. The first column is year, followed by month, day and sea level value (usually in mm). Looking at the .dat equivalent, the file contains a lot more metadata, but also has a mixture of strings and numbers that are hard to import into Matlab. Of course the netCDF file has both the metadata and ease of import.
As an example, let’s try download data from a couple Pacific tide gauges to look at tsunami propagation. The largest earthquake (magnitude 9.5) of the 20th century occurred on May 22, 1960 off the coast of South Central Chile. It generated one of the most destructive tsunamis in the Pacific. Near the generating area, both the earthquake and tsunami were extremely destructive, particularly in the area from Concepcion to Isla Chiloe. Huge tsunami waves measuring as high as 25 m arrived within 10 to 15 minutes after the earthquake, killing at least 200 people. It then took the tsunami about 15 hours to travel the 10,000 miles to Hawaii.
Let’s see how long it took to get to Midway (about 2130 kilometers). Since it’s a fast-moving wave (relatively), we’ll need hourly data. The data are available from the UH Sea Level Center (UHSLC), but it’s already been downloaded and put in our class directory. First, we access the data using the absolute path, and assign the data to the variables H (for Honolulu) and M (for Midway).
>> H = load (‘/home/iniki/ges/ocn363-data/ssh/h057.csv’);
>> M = load (‘/home/iniki/ges/ocn363-data/ssh/h050.csv’);
The data are presented in columns, with the first four being time (year, month, day, hour) and the last (fifth) containing the sea level in millimeters. We can save this final column as a variable:
>> hono_clev = H(:,5);
>> midw_clev = M(:,5);
Next, we can look for the dates of the tsunami. To do this, we’ll use the Matlab “find” function. Note that this returns the index (row,column) of the constraint, not the number itself. Here we look for the times around May 22-25, 1960:
>> I = find ( H(:,1) == 1960 & H(:,2) == 5 & H(:,3) > 21 & H(:,3) < 26 );
>> J = find ( M(:,1) == 1960 & M(:,2) == 5 & M(:,3) > 21 & M(:,3) < 26 );
% 4. Make plot, Honolulu as red circles and line; Midway as green
% squares and line
% Note: we only plot the range given from above (May 22-25 1960)
plot(hono_clev(I),’ro-‘)
hold
plot(midw_clev(J),’gs-‘)
% 5. Pimp the plot; add legend and labels
title(‘Sea level following May 1960 Chilean Earthquake’)
xlabel(‘hour’)
ylabel(‘sea level [mm]’)
legend(‘Honolulu’,’Midway’)
% 6. find when sea level max, assume this is arrival time of tsunami
% Note: here we look for the sea level when it is maximum at each
% site
time_hono = find ( hono_clev(I) == max(hono_clev(I) ));
time_midw = find ( midw_clev(J) == max(midw_clev(J) ));
% 7. compute speed
2180 / ( time_midw - time_hono )
% 8. define time to be linearly increasing
% Note: it’s difficult to represent day/month/year as a linearly
% increasing number. Matlab (like other s/w) proides a
% continuously increasing time defined as “days since some
% reference date” (or “hours since…”). For Matlab this is
% days since Jan 1, 0000.
% The two commands used to convert between the two are
% “datenum” and “datevec”. If you have day/month/year and
% want “days since”, use datenum. That’s exactly what we
% have here.
% Use datenum to get “days since Jan 1, 0000”; pass it year
% month, day and hour (minutes and seconds are 0, 0)
hono_time = datenum(H(:,1),H(:,2),H(:,3),H(:,4),0,0);
midw_time = datenum(M(:,1),M(:,2),M(:,3),M(:,4),0,0);
% 9. we can now replot and take advantage of linear time to make our
% x-axis much better:
figure
plot(hono_time(I),hono_clev(I),’bo-‘)
hold
plot(midw_time(J),midw_clev(J),’gs-‘)
datetick(‘x’,’mmm-dd HH-MM’)
ylabel(‘Sea level (mm)’)
xlabel(‘May 1960’)
title(‘Tide gauge sea level’)
legend(‘Honolulu’,’Midway’)
print -dpng tide_plot.png
% 9. we can now replot and take advantage of linear time to make our
% x-axis much better:
figure
plot(hono_time(I),hono_clev(I),’bo-‘)
hold
plot(midw_time(J),midw_clev(J),’gs-‘)
datetick(‘x’,’mmm-dd HH-MM’)
ylabel(‘Sea level (mm)’)
xlabel(‘May 1960’)
title(‘Tide gauge sea level’)
legend(‘Honolulu’,’Midway’)
print -dpng tide_plot.png
Let’s continue with an example using the Honolulu tide gauge (station ID = 057). Download the Honolulu tide gauge daily dataset (name=d057.csv), and start Matlab. To load the dataset, be aware of where you are running Matlab and where the dataset is (recall our discussion about paths). If the dataset is in your Downloads directory, and you are running Matlab from your home directory, you’ll have to specify the path:
>> load Downloads/d057.csv
>> h=d057;
>> h;
>> figure()
>> plot(h(:,4)) % creates a figure plotting the sea level in mm (all rows, column 4)
>> min(h(:,4)) % finds the bad values as the minimum values
>> i=find(h(:,4)==-32767); % finds the indexes with the bad values -32767
>> whos
>> h(i,4)=NaN; % “masks” the bad values with NaN (Not a number)
>> mean(d057(:,4)); % calculating mean not-excluding NaN
>> nanmean(h(:,4)); % calculating mean excluding NaN
>> mymean=nanmean(h(:,4));
Create a date vector to plot versus the sea level
>> t = datenum(h(:, 1:3));
>> whos t % to see what is t in your file
>> figure()
>> plot(t, h(:,4))
>> datetick(‘x’,’yyyy’) % plots
>> datestr(t); % transforms matlab format to time format
>> g = [t h(:,4)]; % concatenate time and sea level in one matrix
>> find(isnan(g(:, 2))) % find the NaNs in the second column, like we did previously
>> help sin % gives you information on the function sin
>> i = find(isnan(g(:,2)));
>> g(i, :) = []; %gets rids of all NaN data
>> whos h
>> whos g
>> figure()
>> plot(g(:,1), g(:,2))
Now we will do a least-square fit of the data
>> A = polyfit(g(:,1), g(:,2), 1); %use the matlab function polyfit
>> W = polyval(A, g(:,1));
>> hold on % plot the trend on top of our data
>> plot(g(:,1), W, ‘r’) % plot the line on top of the line
>> title(‘HNL sea level, blue=gauge, red=linear fit (1.48mm/yr)’)
>> xlabel(‘date’)
>> ylabel(‘level (mm)’)
>> datetick(‘x’, ‘yyyy’)
>> g = [t h(:,4)]; % concatenate time and sea level in one matrix
>> find(isnan(g(:, 2))) % find the NaNs in the second column, like we did previously
>> help sin % gives you information on the function sin
>> i = find(isnan(g(:,2)));
>> g(i, :) = []; %gets rids of all NaN data
>> figure()
>> plot(g(:,1), g(:,2))
Now we will do a least-square fit of the data
>> A = polyfit(g(:,1), g(:,2), 1); %use the matlab function polyfit
>> W = polyval(A, g(:,1));
>> hold on % plot the trend on top of our data
>> plot(g(:,1), W, ‘r’) % plot the line on top of the line
>> title(‘HNL sea level, blue=gauge, red=linear fit (1.48mm/yr)’)
>> xlabel(‘date’)
>> ylabel(‘level (mm)’)
>> datetick(‘x’, ‘yyyy’)
If you want to have a 5th and 20 degree polynomial
>> A = polyfit(g(:,1), g(:,2), 5); %use the matlab function polyfit
>> W = polyval(A, g(:,1));
>> hold on % plot the trend on top of our data
>> plot(g(:,1), W, ‘g’) % plot the line on top of the line
>> A = polyfit(g(:,1), g(:,2), 20); %use the matlab function polyfit
>> W = polyval(A, g(:,1));
>> hold on % plot the trend on top of our data
>> plot(g(:,1), W, ‘y’) % plot the line on top of the line
Now let’s check fluctuations at higher frequencies
>> figure()
>> plot(g(:,1), g(:,2))
>> title(‘HNL sea level, blue=gauge, red=linear fit (1.48mm/yr)’)
>> xlabel(‘date’)
>> ylabel(‘level (mm)’)
>> datetick(‘x’, ‘yyyy’)
Let’s plot 4 years of data instead of the whole time series In Matlab type:
>> help reshape
>> 41303/1461 % the length of the data divided by four years gives you 28
>> k = reshape(h(1:28*1461,4),1461, 28); % rearrange the time series every four years, you will >> get a matrix of 28 rows, 1461 columns
>> i=find(k(:)==-32767); % finds the indexes with the bad values -32767
>> k(i)=NaN; %
>> figure()
>> plot(k(:, 1)) % plots the first four years of the data
>> hold on
>> plot(k(:,2)) % plots from year 8 to year 12 of data
>> plot(k(:,3))
>> plot(k(:,4))
We can plot the average every 4 years so we can get the yearly seasonality
>> avk = nanmean(k,2);
>> figure()
>> plot(avk)
In the homework, you must get hourly Honolulu sea level, everything is similar as previously done except the reshape command
>> k = reshape(h(1:24*41303), 24, 41303);
>> avk = nanmean(k,2);
B. Acoustic data¶
In the previous homework, we learned how to plot time series of tide gauge data. We found a rich content, at different periods: tidal oscillations at 12 and 24 hours, yearly variations corresponding to seasons, and if you plotted multi-year data, some gauges will show variations due to El Nino and global warming.
Clearly, tools are needed to separate such processes happening at vastly different time scales. This is what spectral analysis is about. In short, it provides a way to convert from the time domain to frequency.
We will begin with the most intuitive time series: sound files (click on score to listen; save file in your matlab directory).
This “audio” file (downloaded from youtube), obtained by feeding the signal of a microphone to an Analog-to-Digital converter (ADC), is perhaps the most common time series, used in music CD. Each sample is 16 bits (from -32768 to +32767), and standard sampling is 44100 samples/seconds.
We can read this file in Matlab using the wavread function:
>> y=wavread(‘twinkle.wav’); % read .wav file scaled to [-1:1]
>> audioinfo(‘twinkle.wav’) %reads the .wav file and give you the file information
>> !file twinkle.wav %reads a .wav file outside matlab (!) allows you to get out of matlab for one line
>> FS=44100; % set the sampling frequency
>> NBITS=16; % set de bit depth of the samples
>> whos
Name Size Bytes Class Attributes
FS 1x1 8 double
NBITS 1x1 8 double
y 1754568x1 14036544 double
Note that matlab converted the “16 bits” to an internal format “double”, which allows computations. Let’s plot the time series:
>> plot((1:1754568)/FS,y) % create x-axis in units of s
>> title(‘twinkle time series’)
>> xlabel(‘time (seconds)’)
>> ylabel(‘normalized amplitude’)
Your plot should look like this:
Let’s zoom on it (1 second):
>> axis([3 4 -0.5 0.5])
and again (1/100 second)
>> axis([3 3.01 -0.5 0.5])
Visually, it is very similar to a plots of tide gauge time series: variations over long periods, short period oscillations. However, when listening, our brain is doing more than decoding the time series: it identifies easily the contents of the score, telling an A from an F. This is the objective of spectrum analysis.
We will skip the theory, and take an empirical approach using Matlab functions. First, we will plot what is called the spectrogram:
>> figure
>> specgram(y,FS/2,FS) % for each 1/2 sec estimate spectrum
>> axis([0 40 400 800]) % music spans A440 to A880
>> xlabel(‘time (seconds)’)
>> ylabel(‘frequency (Hz)’)
>> caxis([-20 20]) % scale the color axis
>> colorbar % put colorbar (scale in dB = 10log10)
>> title(‘twinkle spectrogram’)
The score is simply a transcription of the spectrum for the musician. Music for humans is recorded through its spectrogram, not its time series, because the human ear is a spectrum analyzer, not a time-series sampler. Applying spectrum techniques to any time series (like sea level) is therefore a powerful tool to enhance our understanding of their content.
The spectrogram tells us more than the score, however. Its frequency ranges from 0 to FS/2, half the sampling frequency (the CD sampling frequency was chosen at >40 kHz because the human ear goes to 20 kHz). Change the scale:
>> axis([0 40 0 1e4])
>> caxis([-20 40])
The “fundamental” sounds between 400 and 800 Hz, are replicated as part of the melody, to much higher frequencies, at 2*f, 3*f, … n*f. These are called the “harmonics” of the fundamentals, and their relative strength is what makes a sound specific to an instrument.
There is an important compromise to make when using spectrogram: the spectral resolution dF (the vertical axis) and the temporal resolution dT (the horizontal axis) are linked. Require a higher spectral resolution, you will lose in temporal resolution, and vice-versa. The spectral and temporal resolutions are linked with dT*dF=1. Let’s experiment with this. In the original spectrogram, we chose to have a temporal resolution using FS/2 samples, thus of 1/2 seconds. The frequency resolution was thus 2 Hz. Convince yourself by doing:
>> axis([0 40 400 410])
Let’s try a shorter averaging time of 0.05 seconds using FS/20 samples each:
>> specgram(y,FS/20,FS)
>> axis([0 40 400 800])
The temporal resolution is now very fine, but the gain is mitigated by a much coarser frequency resolution of 20 Hz, as you can see with:
>> axis([0 40 400 500])
Likewise, if we take a longer averaging time of 5 seconds using 5*FS samples each:
>> specgram(y,5*FS,FS)
>> axis([0 40 400 800])
The frequency resolution is now very fine, but the temporal resolution almost useless for the process that we are studying - we cannot distinguish individual notes:
>> axis([0 40 400 402])
The choice of temporal vs. spectral resolution depends on what you want to analyze. Sometimes for process with fast temporal changes, you may need to sacrifice spectral resolution. At the other end, stationary processed that do not change in time (except for the sinusoidal variations of their components) do not require temporal resolution. At this extreme, the interest is in the highest spectral resolution possible. So, let’s compute the spectrogram using the entire length of the record (saving the “spectral data” in array x, don’t forget the “;”):
>> x=specgram(y,1754568,FS);
There is now ultra-fine spectral resolution of about 0.025 Hz which is FS/1754568, as seen with:
>> axis([0 40 400 402])
Since in such plot there is no information in the x-axis, it is customary to plot the spectral data in a different way. “x” is called the Fourier transform of the time series; it is a complex number. To plot it, we need to take its absolute value, and its logarithm (because the scale of variations if very large):
>> plot(0:FS/1754568:FS/2,20*log10(abs(x)))
>> xlabel(‘frequency (Hz)’)
>> ylabel(‘variance (dB)’)
>> title(‘twinkle power spectrum’)
As a forensic analyst, you notice that while the file was uploaded to youtube at the standard 44,100 Hz sampling rate, its author actually initially sampled it at 32,000 Hz, because there is no spectral content above 16 kHz! let’s scale to your useful fundamentals frequency band:
>> axis([400 800 -20 80])
Now you see the exact notes of the A-major scale played by the violonist: A-B-Csharp-D-E-F-sharp. But you don’t have any information on the temporal evolution of the score.
Such techniques are extremely useful for analyzing natural sounds in the ocean, and especially to track marine mammals. An excellent collection of examples including time series and spectrograms is given at http://ocr.org/sounds/marine-mammal-audiographs/. Another cool database of deep ocean sounds is at http://www.listentothedeep.net/acoustics/ (this needs a flash player). SOEST has their own observatory 100 km north of Oahu, called the ALOHA Cabled Observatory (ACO), with hydrophones (underwater microphones) connected in real-time to the shore station by fiber optics cable.
Time series of Doppler signal
Spectrum and spectrogram are important in navigation, both underwater acoustic navigation and radio satellite navigation at the surface and for aircrafts (Global Positioning System - GPS). We now analyze a simple case of navigation data, the Doppler signal of the horn of a passing car, which is identical in concept with the Doppler radio signal of a passing satellite.
>> y=wavread(‘dopplercar.wav’);
>> FS=44100;
>> plot((1:228111)/FS,y) % create x-axis in units of s
>> title(‘Doppler time series’)
>> xlabel(‘time (seconds)’)
>> ylabel(‘normalized amplitude’)
>> figure
>> specgram(y,FS/20,FS)
>> xlabel(‘time (seconds)’)
>> ylabel(‘frequency (Hz)’)
>> caxis([-20 20]) % scale the color axis
>> colorbar % put colorbar (scale in dB = 10log10)
>> title(‘Doppler spectrogram’)
>> caxis([-20 20])
The fundamental frequencies are 425 and 475 Hz. Note how complex the spectrum is: two fundamental frequencies, full of harmonics. A horn is designed to be heard over noise. A strong Doppler shift between incoming and outgoing horn on the car is observed. To estimate the Doppler shift:
>> axis([2.5 5 2000 4000])
>> ginput(2)
then click on a spectrum line before and after the passage of the car, for example the one at 3000 Hz:
ans =
3.08 3049.71
4.62 2821.64
That shows that the mean frequency (of that harmonic) was 2935 Hz, and the frequency shift was +-115 Hz. There is an approximate relationship between the Doppler shift Df of a moving sound or radio source of frequency f, and its velocity V, as long as V << C (much smaller):
V = C * Df / f
where C is the propagation speed of the wave (300,000 km/s for light/radio, 340 m/s for sound in air, 1500 m/s for sound in water). We can thus infer:
V = 13.3 m/s = 48 km/h = 30 mph.
The time at which the Doppler shift goes through 0 is called the “point of closest approach”. Estimating this point is the base for satellite navigation. (Side exercise: can you estimate the Doppler shift of a transmitter of 1.572 GHz on a low altitude satellite (GPS signal)?).
C. Waves¶
D. CTD from HOT¶
We will continue to look at time-series but now expand to variables that are now functions of depth as well. In addition, we will look at multiple measurements made at the same location. The data we will use come from the Hawaii Ocean Time-series (HOT) program. This is an on-going effort to monitor the North Pacific Subtropical Gyre via repeat ship-based measurements. An overview of HOT is available at http://hahana.soest.hawaii.edu/hot/ and at http://www.soest.hawaii.edu/HOT_WOCE/.
For our next exercise, we will make a plot of temperature vs. depth over time (two independent variables). To do this, we will make use of color-shaded contours (x-axis is time, y-axis depth, colors correspond to values). First, let’s run through an example of contour plotting in Matlab. There are a few different ways of doing this, including “contour”, “contourf”, “pcolor” and “imagesc”. The basic syntax is the same for all, e.g., contour(variable), or contour(x,y,variable).
Let’s start with a simple example. We define a variable, A, that is a function of y and z defined as follows:
y = sin(x)
z = cos(x)
To construct these in Matlab:
>> x = ( 0 : 0.01 : 2.0 * pi );
>> y = sin(x);
>> z = cos(x);
>> A = y’ * z;
We now have several ways to plot this: contour(x,x,A) contourf(x,x,A)
pcolor(x,x,A) imagesc(x,x,A)as shown below: 
We will do something similar, except we will have depth along the y-axis, time along the x-axis, and the contours will be of temperature. For the HOT data, we will use the cruise-averaged variables in the file /home/iniki/ges/ocn363-data/hot/HOTyears01-29.ctd. This file contains the temperature, salinity and oxygen values as measured by the CTD casts averaged over the entire cruise. So, there is essentially one depth-profile per cruise. The file is laid out where each line is a different depth, with each subsequent cruise measurements appended. We want to create a color-shaded plot of a variable (temperature) as a function of depth (y-axis) and time (x-axis). The task is complicated by the fact that there are, potentially, different number of depths measured during each cruise.
To deal with this, the process will be:
create an empty matrix that has different times in the columns and different depths in the rows
fill in each “cell” or “row/column” with a value if present, otherwise leave it blank
contour the matrix
Below is a schematic of our matrix:
To add some detail to our list:
find the maximum depth, and create rows from 0 to the maximum incremented by 2
find the total number of cruises and create columns from 1 to the total (note finding when a new cruise starts use diff)
loop through each row in the data
Here is the script:
% Now we try make our HOT plot
% 1. Understand the data set.
% a. The data are in /home/iniki/ges/ocn363-data/data/hot/HOTyears01_29ctd.dat
% b. This data set represents mean observations from CTD casts for each
% cruise
% c. Note (via “more”) the data set is layed out in columns, with the
% first being day since January 1, 1988 (start of program); second
% column is depth; third is density; fourth is temperature; fifth
% is salinity; sixth is oxygen. (Note that depth is in negative
% decibars, approximately meters, to make plotting easier; density
% is also negative, but I’m not sure why)
% d. The data get written after each cruise, so the rows are depth
% for a particular cruise, then the next cruise data is appended
% at the end.
% 2. Goal is to make a color-shaded plot of temp (or salt or O2) vs. time
% and depth
% 3. The process is we will make an empty matrix where each row is a depth
% bin, and each column is a time bin. The value in each particular
% row/column cell is the value measured. For example, row ten column 5
% will contain temperature at the 10th depth level and 5th time.
% load the data
hot = load (‘/home/iniki/ges/ocn363-data/hot/HOTyears01_29ctd.dat’);
% convert depths to be positive numbers (read all rows, second column)
hot(:,2) = -1.0 * hot(:,2);
% note that for our plot, we want to go to some maximum depth, but we
% don’t know what that is; each cruise will go to some different maximum
% depth. so, first find the maximum depth, and create a new depth
% index from the surface (0) to the maximum depth (max_depth) with an
% increment of 2.
max_depth = max ( hot(:,2) );
depth = [0:2:max_depth];
% now find the number of depth levels
ndep = length(depth);
% So, we now know our matrix will have “ndep” rows. Next we need to
% find the number of columns, or number of times. This will be equal
% to the number of cruises. To find the number of cruises, we search
% for when the day (column one) changes value.
icr = find ( diff ( hot(:,1) ) );
% now find the index to the start of each cruise (ic1)
% the index to the end of each cruise (ic2) and the
% total number of cruises (ncru)
ic1 = [1; icr+1];
ic2 = [icr; length(hot)];
ncru = length(ic1);
day = hot(ic1,1);
time = 1988 + day / 365;
% make a matrix (temp) that has “ndep” rows and “ncru”
% columns (depth will be rows and cruise/date is the
% columns)
temp = ones(ndep, ncru) + NaN; % Make a matrix w/ cruise
% finally, loop over all cruises, and replace the
% NaN elements in temp with real values (temperature is the 4th column)
for ic = 1:ncru,
j = hot ( ic1(ic) : ic2(ic), 2 ) / 2 + 1;
temp(j,ic) = hot ( ic1(ic) : ic2(ic), 4 )’;
end
% now plot
contourf(time,depth,temp)
set(gca,’YDir’,’Reverse’)
title (‘HOT temperature’)
colorbar
print -dpng hot_temp.png
E. ADCP¶
Time series of currents
Another application of Doppler shift is the measurement of currents, using an Acoustic Doppler Current Meter (ADCP). An ultrasound is broadcast through the water column, and the spectrum of the echoes are analyzed. If there is a relative water motion between the instrument and the water, the echo will be subject to a Doppler shift, directly related to the water velocity.
Data from a 300 kHz ADCP deployed off Kaena point with format description are posted at [].
The ADCP was upward looking on a mooring similar to this:
Right-click to download the data, save on your workstation, and load into matlab:
>> load a2_adcp.mat
>> whos
Name Size Bytes Class
d 1x1 68914870 struct
We have here a new matlab concept, analog to the concept of folders for files. Instead of giving all the variables separately, they are bundled into a “structure”. To look at the content of the structure, simple type the structure name:
>> d
d =
vel1: [34x15083 double]
vel2: [34x15083 double]
etc…
There are a lot of parameters recorded, but we’ll be mostly interested in velocity (vel1,2,3,4), the depth of the bins, the amplitude of the signal of each beam (amp1,2,3,4), the temperature, and the heading’ of the mooring. Each parameter is over 34 depth bins of 4 m each, and there are 15,083 time samples, one every 20 min (72 per day). We will explore the content of the data, plot time series, and spectrum.
Each variable can be accessed by pre-pending its name with d. such as
>> d.depth
First, look at the amplitude (strength) of the returned signal, in the mean:
>> plot(d.depth,mean(d.amp1’))
>> xlabel(‘range (m)’)
>> ylabel(‘signal strength (dB)’)
>> title(‘beam signal strength’)
The apostrophe after d.amp1 is the “transpose operator”. It changes the matrix d.amp1 from being 34x15083 to 15083x34. This is done so that the mean function operates over time, not depth. The figure shows that the signal decays with distance, from +130dB to 60dB, then suddenly at a range of 100 m, increases to a very high value. This is the echo of the surface. It limits the range at which an ADCP can sample the water.
Next, look at the time series of eastward and northward velocity, at a range of say 38-42 m (bin 10):
>> plot(d.dday-d.dday(1),d.vel1(10,:))
>> title(‘velocity at 40 m’)
>> xlabel(‘days since start’)
>> ylabel(‘velocity (m/s); east=blue, north=red’)
>> hold on
>> plot(d.dday-d.dday(1),d.vel2(10,:),’r’) % use different color for north component
Zoom in:
>> axis([50 55 -0.4 0.4])
We see a rich spectral content of the time series, with slow variations at 30-100 scale, and fast variations due to tides and internal waves. Let’s attempt a spectrogram:
>> FS=72 % unit is in days
>> figure
>> specgram(d.vel1(10,:),7*FS,FS) % spectrum window 7 days
>> xlabel(‘days since start’)
>> ylabel(‘frequency (cpd)’)
>> caxis([-20 20])
>> axis([0 200 0 6])
Signals at 1 cycle per day (cpd; called diurnal tide) and 2 cpd (called semidiurnal tide) dominate. There is a hint of intermittent signal at 3 cpd and 4 cpd, harmonics of the main tides. Clearly, we need more frequency resolution. Let’s do a spectrum over the whole time series (i.e. renounce to temporal resolution, get maximum frequency resolution):
>> specgram(d.vel1(10,5000:end),100*FS,FS)
Time series of currents
>> load a2_adcp.mat
>> FS=72 % unit is in samples/day
Let’s try to look at the velocity data
>> more on % this will enable scoll/stop with the space bar; type letter q to exit more
>> d.vel1(10,:) % look at the data of depth bin 10 of the east velocity
You see NaN, which stands for Not-a-Number - the start of the file has bad data. Keep hitting the space bar until you see real numbers coming, something like:
Columns 92 through 98
NaN NaN NaN NaN NaN NaN -0.0970
So good data starts at time step 98. Likewise, if you type:
>> d.vel1(10,15001:end)
You see bad data NaN starting again at 15063:
Columns 57 through 63
0.0520 0.0510 0.0450 0.0490 0.0680 0.1350 NaN
Columns 64 through 70
NaN NaN NaN NaN NaN NaN NaN
This happens often with real data. I left them in the file, so you can see how to get rid of the bad data in an array and not be taken by surprise if it happens to you. In all subsequent uses of d.vel1, refer to it as
>> ngood=98:15062;
>> n=size(ngood,2);
>> d.vel1(10,ngood)
to simply skip the bad data at the start and end, with the total length of the record called n.
Now let’s attempt a spectrum averaged over the whole period (i.e., not a time dependent spectrogram, a line spectrum like we did at the end of the lecture on spectra):
>> specgram(d.vel1(10,ngood),n,FS);
>> ylabel(‘frequency (cpd)’)
>> caxis([10 40])
>> axis([0 200 1.7 2.3])
Remember the parameters of the function specgram: TIME_SERIES,NUMBER_SAMPLES_AVERAGED,SAMPLING_FREQUENCY. If you do help specgram, matlab will tell you that there is a new function spectrogram. If you are adventurous, you can read about it; but it is more complicated to use. As a starter, stick with specgram.
Now that we have ensured maximum frequency resolution (at the expense of time resolution), we notice that the signal at 2-cpd called semidiurnal tide is clearly split into several components. Just like we did for the violin sound, we may as well plot this spectrum which has lost time resolution as a line drawing. That requires capturing the spectrum (also called Fourier transform) in an array, take the square of its absolute value to get the spectral variance, and plot it:
>> x=specgram(d.vel1(10,ngood),n,FS);
>> y=abs(x).^2;
>> whos y
Name Size Bytes Class Attributes
y 7483x1 59864 double
>> fscale=0:FS/n:FS/2; % frequency range from 0 to half the sampling frequency
>> plot(fscale,y)
>> xlabel(‘frequency (cpd)’)
>> axis([0 4 0 1e4])
You clearly see two distinct peaks near 2-cpd, and peaks at 1 cpd, 3 and
4 cpd. Let’s try to find their exact frequencies of the 2 cpd peaks:
>> axis([1.7 2.3 0 15000])
>> ginput(2) % click on the top of the two peaks near 2 cpd
ans =
1.9343 9984.4
2.0021 5102.3
>> 24./ans(:,1) % how many hours is this
ans =
12.4074
11.9876
>> (ans-12)*60 % how many minutes is this
ans =
24.4431
-0.7458
The theoretical values, based on astronomy, are 1.93236 and 2.00000. The second is exactly at 12 hours, that is called the solar semi-diurnal tide, because it is synchronous with the sun. The first one is more interesting: it is 12 h 25 min. This is half a lunar day: 28 lunar days take 29 solar days, thus 29/28*12 = 12 h 25 min. It is called the lunar semi-diurnal tide. These two frequencies in the spectrum “beat”: they together create spring and neap tides. In days, the period of the spring/neap cycle is simply 1/(1.9343-2.0021) which comes out to a about 14 days, the fortnightly cycle, as expected.
Spectra like this are noisy: any noise in the data translates into noise in the spectrum. Noise happens at all frequencies. To make spectra easier to analyze, it is customary to “smooth” them.
Reducing noise can be done two different ways: computing a succession of variance spectra over smaller lengths of the time series, and then averaging them, or directly smoothing a single spectrum in the frequency domain. The two methods are equivalent and are readily implemented.
Since we have already the spectrum, we will choose the second one. First, redo the un-smoothed plot in a sub-window, using a log scale:
>> subplot(1,2,1)
>> plot(fscale,10*log10(y))
>> xlabel(‘frequency (cpd)’)
>> axis([0 6 10 50])
We simply average the spectral variance over k consecutive bins, choosing k odd. The following code will do this (make a sketch to understand how it works):
>> k=7;
>> m=size(y,1);
>> l=floor(k/2); % this is half the smoothing window so that k=2*l+1
% for k=3 l=1, k=5 l=2, k=7 l=3 …..
>> for i=l+1:m-l-1
>> z(i)=mean(y(i-l:i+l)); % average over k elements i-l,…,i-1,i,i+1,…,i+l
>> end
>> for i=1:l-1
>> z(i)=mean(y(1:i+l)); % average over less bins near the start
>> end
>> for i=m-l+1:m
>> z(i)=mean(y(i-l:m)); % average over less bins near the end
>> end
>> subplot(1,2,2)
>> plot(fscale,10*log10(z),’r’)
>> xlabel(‘frequency (cpd)’)
>> axis([0 6 10 50])
The spectrum is much cleaner, but there is less resolution. Any additional smoothing would smear the difference between the lunar and solar semi-diurnal peaks, as you can see with zooming on the frequency axis:
>> subplot(1,2,1)
>> axis([1.7 2.3 10 50])
>> subplot(1,2,2)
>> axis([1.7 2.3 10 50])
The art of spectrum analysis, just like the art of making spectrograms, is to find the optimum averaging of the spectrum to make desired processes stand out of the random noise, without overly smearing them. There is no single answer; it is part of experience analyzing data.
Now let’s play with a ship-board ADCP
>> load coast % you will load the coastline
>> xlon = long; % renaming the variables in the file coast
>> ylat = lat;
>> plot(xlon, ylat, ‘k’) % you will see a map of the world
>> hold on % so you can plot more things on top of this figure
>> load /home/iniki/ges/ocn363-data/adcp/ship/mat/00175ave.mat
>> plot(lon, lat, ‘*m’) % you will see the ship track
>> axis([-120 -80 -20 40])
>> axis square % makes the distance between x and y coordinates the same
>> figure()
>> plot(xlon, ylat, ‘k’) % plot the coast
>> hold on
>> quiver(lon, lat, u(:,1), v(:,1), ‘r’) % plots the velocities at surface depth
>> axis equal
>> axis([-120 -90 -10 35])
>> axis([-112 -108 -5 5]) % do a closer zoom to see the currents
>> hold on
>> quiver(lon, lat, u(:,10), v(:,10), ‘g’) % lets plot the currents at 100m depth
>> load coast
>> ylat = lat;
>> xlon = long;
>> load /home/iniki/ges/ocn363-data/drifter/buoy.mat % load the buoy data from iniki server
>> whos % to see what variables are available
>> figure()
>> plot(xlon, ylat, ‘k’)
>> hold on
>> plot(lon, lat, ‘bo’) % you will see buoys floating around the Hawaiian Islands
>> axis square
ARGO buoys
>> load /home/iniki/ges/ocn363-data/drifter/drift.mat
>> whos % to see what variables are available
>> figure()
>> plot(xlon, ylat, ‘k’)
>> hold on
>> plot(Lon, Lat, ‘r.’) % you will see buoys floating around the Hawaiian Islands
>> axis square
>> axis([-230 -100 0 40])
F. HFR¶
Loops. Useful for repetition, instead of repeating a command use a loop. An example of the Fibonacci series with the Matlab command while.
>> N = ?; % choose the last number of the series, Instead of ?, write the number you wish
>> fib=zeros(1,N);
>> fib(1)=1;
>> fib(2)=1;
>> k=3;
>> while k <= N
>> fib(k)=fib(k-2)+fib(k-1);
>> k=k+1;
>> end
Now with the for command:
>> fib=zeros(1,N);
>> fib(1)=1;
>> fib(2)=1;
>> k=3;
>> for k=3:N;
>> fib(k)=fib(k-2)+fib(k-1);
>> k=k+1;
>> end
How to plot vectors in matlab!! Let’s start by downloading the data []. Start the terminal and load the data in the appropriate directory
>> load hfr_2011.mat
>>%This file includes the variables X and Y for the longitude and latitude, U and V for the zonal and meridional currents and T for the time % Plot the mean velocity for the whole time period
>> figure()
>> s=.05;
>> quiver(X,Y,nanmedian(U,3)*s,nanmedian(V,3)*s,0)
>> hold on
>> axis square
>> text(-157.6,20.82,’30 cm/s’)
>> quiver(-157.6,20.8,.30*s,0)
>> grid
>>%Plot just certain months of the year
>> [year,month,day]=datevec(T);
>> figure()
>> subplot(221)
>> fi = find(month==2);
>> s=0.04;
>> quiver(X,Y,nanmedian(U(:,:,fi),3)*s,nanmedian(V(:,:,fi),3)*s,0)
>> axis square
>> grid
>> fi = find(month==5);
>> subplot(222)
>> quiver(X,Y,nanmedian(U(:,:,fi),3)*s,nanmedian(V(:,:,fi),3)*s,0)
>> axis square
>> grid
>> fi = find(month==7);
>> subplot(223)
>> grid
>> quiver(X,Y,nanmedian(U(:,:,fi),3)*s,nanmedian(V(:,:,fi),3)*s,0)
>> axis square
>> grid
>> fi = find(month==11);
>> subplot(224)
>> quiver(X,Y,nanmedian(U(:,:,fi),3)*s,nanmedian(V(:,:,fi),3)*s,0)
>> axis square
>> grid
>>%lets add the titles to each subplot
>>% now let’s make a movie using for loops
>> figure()
>> for i=1:12
>> fi =find(month==i);
>> quiver(X(1:3:end, 1:3:end), Y(1:3:end, 1:3:end), nanmedian(U(1:3:end,1:3:end,fi),3)*s,nanmedian(V(1:3:end,1:3:end,fi),3)*s,0)
>> axis equal
>> axis([-158.4 -157.5 20.65 21.3])
>> hold on
>> text(-157.6,20.82,’30 cm/s’)
>> quiver(-157.6,20.8,.30*s,0)
>> title(sprintf(‘Currents for Month %d, 2011’,i),’FontSize’,14)
>> grid
>> drawnow
>> pause(1)
>> hold off
>> end
Let’s make a color map of the two-year mean currents
>> figure()
>> contourf(X, Y, nanmedian(U,3))
>> axis square
>> axis([-158.4 -157.5 20.65 21.3])
>> hold on
>> title(‘Two year median of surface currents from HFR’)
The Oahu coast is in the following directory that can be accessed from your computer. On your terminal enter ls /home/iniki/ges/ocn363-data/hfr/; you should be able to see three files where oahu.dat is the coast. Copy to the directory of your preference in your computer (cp /home/iniki/ges/ocn363-data/hfr/oahu.dat .).
You can also plot the data from different dimensions
>> figure()
>> pcolor(squeeze(V(:, 20, :)))
>> shading flat
>> colorbar
>> caxis([-6 6])
>> figure()
>> feather(squeeze(U(10,20,:)), squeeze(V(10,20,:)))
>> figure()
>> pcolor(T, Y(:,20), squeeze(U(:,20,:)))
>> shading flat
>> datetick(‘x’,’ mmyy’)
>> figure()
>> pcolor(Y(:,20), T, squeeze(U(:,20,:))’)
>> shading flat
>> datetick(‘y’,’ mmyy’)
>> colorbar()
>> caxis([-6 6])
>> figure()
>> pcolor(X(10,:), T, squeeze(U(10,:,:))’)
>> shading flat
>> datetick(‘y’,’ mmyy’)
>> colorbar()
>> caxis([-6 6])
>>%Let’s make a zoom where we see more energy and find the propagation speed of the perturbation observed.
>> axis ([-158.17 -157.66 734721.13 734750.09])
>> figure()
>> pcolor(X(10,:), T, squeeze(V(10,:,:))’)
>> shading flat
>> datetick(‘y’,’ mmyy’)
>> colorbar()
>> caxis([-6 6])
>>%Let’s make a zoom where we see more energy and find the propagation speed of the perturbation observed.
>> axis ([-158.17 -157.66 734721.13 734750.09])
Select two points in the graph to find the two latitudes and the two times
>> ginput(2)
>> speed = diff(ans)
The first value will give you the degrees of longitude the perturbation moved and the second value will give you the time it took to move. You should get something about 60 km/day or 0.8 m/s.
>> vprog = speed(1) * cosd(21) * 10e7/90/(speed(2) *86400) % in m/s
Histograms and normal distribution
Selecting one point in the HFR spatial domain
>> up = squeeze(U(15, 15, :));
>> figure()
>> subplot(211)
>> hist(up, 32) % you can select the number of “little bars” to separate your distribution
>> title(‘u’)
>> vp = squeeze(V(15, 15, :));
>> subplot(212)
>> hist(vp, 32) % you can select the number of “little bars” to separate your distribution
>> title(‘v’)
Selecting all the data in the HFR spatial domain
>> figure()
>> u0 = U(:);
>> hist(u0, 512)
>> figure()
>> v0 = V(:);
>> hist(v0, 512)
G. Drifters (surface and argo)¶
In this lecture, we will use a large data set of satellite-tracked buoys, to compute the mean circulation around the Hawaiian islands.
World Ocean Circulation Experiment drifting buoy principle
A WOCE buoy before deployment
The WOCE satellite-tracked buoys data for the Hawaiian area is []. It covers more than 10 years. Download it and import it into matlab. Let’s first look at all the buoys:
>> load drift.mat
>> whos
U and V have been pre-computed from the randomly-space ARGOS satellite fixes, and co-located to the interpolated latitude-longitude. The variables of interest to us are:
Name Size Bytes Class Attributes
Lat 418660x1 3349280 double
Lon 418660x1 3349280 double
U 418660x1 3349280 double
V 418660x1 3349280 double
id 418660x1 3349280 double
>> figure(1)
>> plot(Lon,Lat,’.’)
The data set covers 10 N to 30 N and 130 W to 180 W, or 100 deg^2. On average, there are about
>> size(id,1)/(20*50)
ans =
418.6600
points/deg^2. To have a graphic idea of the density of data, let’s look at a 2x2 deg^2 box centered on (-160,20):
>> axis([-161 -159 19 21])
Let’s compute the median currents in this box. First, we need to find the indices (in the 418660x1 matrices) of all the buoys that crossed the box, and compute the median currents in this box:
>> box=find(Lat<20+1 & Lat>20-1 & Lon>-160-1 & Lon<-160+1);
>> size(box,1)
ans = 2773
>> median(U(box))
ans = 0.0242
>> median(V(box))
ans = 0.0701
We use median rather than mean, because the median is robust to outliers (unflagged bad data), while the mean is not. The mean currents in this box are:
>> mean(U(box))
ans = 0.0007
>> mean(V(box))
ans = 0.0580
The means are different from the medians, and you can easily convince yourself that the distributions, especially of V, are not symmetric, skewed, by plotting:
>> hist(U(box),63)
>> hist(V(box),63)
Likewise, the standard deviations are
>> std(U(box))
ans = 0.2418
>> std(V(box))
ans = 0.2187
reflecting the energetic eddies present in the area, much larger than the mean currents.
Of course, our interest is to construct a map of mean currents, not just estimate them at a single point. For this we must define a grid. The grid resolution is a matter of compromise between precision and resolution: finer resolution will yield less buoys crossing each box, thus larger standard errors on the estimation of the mean, while coarser resolution will yield more precise estimation of the mean velocities. We chose a grid resolution of 2x2 deg., which is a good compromise; see (-160,20) above. We define two matrices cenlon and cenlat containing the (lon,lat) of each grid point (center of each box). cenlon is a matrix in which all lines are identical to each other, and cenlat is a matrix in which all columns are identical to each other. We populate the cenlon matrix by creating the first line with longitudes [-179:2:-131], and then replicating it for all latitudes. Likewise, we populate the cenlat matrix by populating the first column with latitudes [11:2:29], and the replicating it for all longitudes. In practice:
>> cenlat=zeros(10,25);
>> cenlat(:,1)=[11:2:29];
>> for j=2:25
>> cenlat(:,j)=cenlat(:,1);
>> end
>> cenlon=zeros(10,25);
>> cenlon(1,:)=[-179:2:-131];
>> for i=2:10
>> cenlon(i,:)=cenlon(1,:);
>> end
Here they are, filled:
cenlon =
-179 -177 -175 -173 -171 -169 -167 -165 -163 -161 -159 -157 -155 -153 -151 -149 -147 -145 -143 -141 -139 -137 -135 -133 -131
-179 -177 -175 -173 -171 -169 -167 -165 -163 -161 -159 -157 -155 -153 -151 -149 -147 -145 -143 -141 -139 -137 -135 -133 -131
-179 -177 -175 -173 -171 -169 -167 -165 -163 -161 -159 -157 -155 -153 -151 -149 -147 -145 -143 -141 -139 -137 -135 -133 -131
-179 -177 -175 -173 -171 -169 -167 -165 -163 -161 -159 -157 -155 -153 -151 -149 -147 -145 -143 -141 -139 -137 -135 -133 -131
-179 -177 -175 -173 -171 -169 -167 -165 -163 -161 -159 -157 -155 -153 -151 -149 -147 -145 -143 -141 -139 -137 -135 -133 -131
-179 -177 -175 -173 -171 -169 -167 -165 -163 -161 -159 -157 -155 -153 -151 -149 -147 -145 -143 -141 -139 -137 -135 -133 -131
-179 -177 -175 -173 -171 -169 -167 -165 -163 -161 -159 -157 -155 -153 -151 -149 -147 -145 -143 -141 -139 -137 -135 -133 -131
-179 -177 -175 -173 -171 -169 -167 -165 -163 -161 -159 -157 -155 -153 -151 -149 -147 -145 -143 -141 -139 -137 -135 -133 -131
-179 -177 -175 -173 -171 -169 -167 -165 -163 -161 -159 -157 -155 -153 -151 -149 -147 -145 -143 -141 -139 -137 -135 -133 -131
-179 -177 -175 -173 -171 -169 -167 -165 -163 -161 -159 -157 -155 -153 -151 -149 -147 -145 -143 -141 -139 -137 -135 -133 -131
cenlat =
11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13
15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15
17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17 17
19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19
21 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21
23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23
25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25
27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27
29 29 29 29 29 29 29 29 29 29 29 29 29 29 29 29 29 29 29 29 29 29 29 29 29
We are now ready to compute the mean currents:
>> meanU=zeros(10,25);
>> meanV=zeros(10,25);
>> for i=1:10
>> for j=1:25
>> box=find(Lat<cenlat(i,j)+1 & Lat>cenlat(i,j)-1 & Lon>cenlon(i,j)-1 & Lon<cenlon(i,j)+1);
>> meanU(i,j)=mean(U(box));
>> meanV(i,j)=mean(V(box));
>> end
>> end
It takes a little while to compute, but then it is trivial to plot the mean currents; we add Oahu coastline for reference.
>> quiver(cenlon,cenlat,meanU,meanV)
>> axis([-180 -130 10 30])
>> hold on
>> plot(coast(:,1),coast(:,2),’.r’)
>> title(‘mean currents Tropical Central Pacific from WOCE drifters’)
>> xlabel(‘Longitude’)
>> ylabel(‘Latitude’)
The mean currents show the well-defined southern part of the North Pacific gyre, with the westward North Equatorial Current (NEC) intensifying from the latitude of Hawaii, down to a maximum at 12N.
H. Satellite measurements¶
We will now investigate global-scale data measured by satellite. Earth-observing satellites have been measuring various parameters such as sea level, sea surface temperature, ocean color, etc. The benefits of satellite-derived data are that they are typically global in extent. The downsides are spatial and/or temporal resolution is sometimes coarse (e.g., 100 kilometers and weekly). Actually, most satellites have a tradeoff between the spatial and temporal resolution. For example, a geostationary satellite can get very high-resolution images, but maybe only once a month. Other limitations include cloud contamination, near-shore effects, and that typically satellites can only “see” the ocean surface.
Nevertheless, satellite data can be very powerful for studying global phenomena. We have a set of data on local disk for these exercises, but more data are available at the various NASA data centers, for example the Physical Oceanography Distributed Active Archive Center (PODAAC; https://podaac.jpl.nasa.gov).
First, let’s have a look at our local data. As in previous exercises, these are all available at /home/iniki/ges/ocn363-data in the subdirectories clims, ssh, misc, ersscat, etc.
Subdirectory |
Filename |
Description |
|---|---|---|
clims |
ers_clim.mat |
Monthly mean seasonal cycle of wind speed (m/s) ERS-1 |
rain_clim.mat |
Monthly mean seasonal cycle of rain (mm/day) CMAP |
|
ssh_clim.mat |
Monthly mean seasonal cycle of sea level (cm) TOPEX |
|
sst_clim.mat |
Monthly mean seasonal cycle of SST (degC) from AVHRR |
|
ersscat |
m####.mat |
Monthly wind speed (m/s) in 21 files (1/97 - 9/98) |
ssh/mat |
dev96_99.mat |
4 years of monthly mean sea level from TOPEX |
global_sl.mat |
11 years of annual mean sea level from TOPEX (1993-2003) |
|
aviso_ann.mat |
20 years of annual mean sea level from AVISO (1993-2012) |
|
misc |
avhrr_veg.mat |
6 years of monthly mean vegetation index (1995-2000) |
CMAP_precip |
10 years of monthly mean rainfall (1990-1999) |
|
SPac_bathy.mat |
High-res bathymetry for the South Pacific |
|
TOPEX_mm.mat |
9 years of monthly mean sea level (1993-2001) |
Since these data have large spatial extents, we will first make color-shaded “maps” of various properties. To help, we will also plot coastlines. Note that matlab includes a database of coastlines, called “coast” in a easy to format. One caution is that the variables included in “coast” are long and lat, and if you later load a data set with these variable names, they will be overwritten. Therefore, it is usually helpful to immediately redefine these.
>> load coast;
>> xlon = long;
>> ylat = lat;
>> plot(xlon,ylat,’k’)
For the first example, we will make a color-shaded plot of sea surface temperature (SST) from the Advanced Very-High Resolution Radiometer (AVHRR).
>> load /home/iniki/ges/ocn363-data/misc/avhrr_sst
>> contourf(lon,lat,squeeze(sst(1,:,:))’)
>> colorbar
>> hold
>> plot(xlon,ylat,’k’)
>> plot(xlon+360,ylat,’k’)
We can also select individual points from the satellite record. In this case, let’s have a look at sea level from the AVISO product. We can load in the annual mean data and make a time-series plot at a specific location.
load /home/iniki/ges/ocn363-data/ssh/mat/aviso_ann
% plot time series for some grid point
plot(time,ssh(:,500,400),time,ssh(:,500,400),’o’)
% make a global map of mean sea level anomalies
h = reshape ( ssh, 20, 1080 * 915 );
ssh_mean = nanmean ( h’ );
plot ( time, ssh_mean, ‘o’ );
% fit a straight line
b = polyfit ( time, ssh_mean, 1 );
% slope (sea level rise in cm/year)
b(1)
y = polyval ( b, time );
hold
plot ( time, y, ‘r’ )
% can also do y = m * x + b, where m=b(1) and b=b(2)
% now compute the rate of sea level rise over the whole globe
% note we just use a subset since the whole data set takes a
% long time to compute
for i = 301 : 841
for j = 327 : 590
x = squeeze ( ssh(:,i,j) )’;
b = polyfit ( time, x, 1 );
htrend(i-300,j-326) = b(1);
end
end
contourf(lon(301:841),lat(327:590),htrend’)
hold
plot(xlon,ylat,’k’)
colorbar
caxis([-1 1])
Now do a comparison with tide gauges:
% load same annual mean sea level data set
load /home/iniki/ges/ocn363-data/ssh/mat/aviso_ann
% find a point
J = find ( lat > 21.2 & lat < 21.4 );
I = find ( lon > 202.2 & lon < 202.4);
h = squeeze(ssh(:,I,J));
newh = h - nanmean(h);
plot(time,newh,’o’)
hold
b=polyfit(time,newh’,1)
y = polyval(b,time)
plot(time,y,’r’)
% compare to tide gauge
load /home/iniki/ges/ocn363-data/ssh/dat/honolulu.dat
time2 = honolulu(:,1);
ssh = honolulu(:,2);
J = find ( time2 >= time(1) & time2 < time(end) );
ssh2 = ssh(J);
I = find(ssh2==-999);
ssh2(I) = NaN;
avessh = nanmean(ssh2);
newh=0.1*(ssh2-nanmean(ssh2));
newh(I)=0.0;
plot(time2(J),newh)
b = polyfit(time2(J),newh,1)
y = polyval(b,time2(J));
plot(time2(J),y,’g’)
% now look at entire record
ssh2 = 0.1 * ( ssh - nanmean(ssh));
plot(time2,ssh2)
b = polyfit(time2,ssh2,1);
y = polyval(b,time2);
plot(time2,y,’k’)
Now plot some climatologies.
% load climatologies
load /home/hokulea/ges/ocn363/data/clims/sst_clim
A=squeeze(sst(1,:,:));
pcolor(lon,lat,A’)
shading interp
colorbar
% 3. plot animation of sst anomalies
clear all
load coast;
xlon = long;
ylat = lat;
load /home/hokulea/ges/ocn363/data/ssh/mat/dev96_99;
format bank
year = floor(t);
month = floor ( ( t - year ) * 12 ) + 1;
cmon = [‘Jan’;’Feb’;’Mar’;’Apr’;’May’;’Jun’;’Jul’;’Aug’;’Sep’;’Oct’;’Nov’;’Dec’];
for i = 1 : 48
pcolor ( x, y, squeeze ( dev ( i,:,:) ) );
shading interp;
caxis ( [-50 50 ] )
colorbar
axis image
%date=[num2str(month(i)),’/’,num2str(year(i))];
date = [ cmon(month(i),:),’,’, num2str(year(i))];
text(240,50,date)
hold on
plot(xlon,ylat,’k’)
plot(360+xlon,ylat,’k’)
hold off
command=[‘print -djpeg fig’,num2str(i+10)];
eval(command);
pause
% the command pause waits for you to hit the “enter” key before going on
% you can also specify the amount of time pause waits before automatically
% going on:
% pause(1)
end
% now load another data set that contains monthly mean ssh anomalies from Topex/Poseidon
% for the 97-98 El Nino/La Nina period:
clear
load /home/hokulea/ges/ocn363/data/ssh/mat/dev96_99
whos
README
% make a movie of maps of monthly mean ssh anomalies
% we can store in a matlab variable the figures created during the loop,
% so that we can play back the movie as many times as we want without having to redo
% the loop each time:
help movie
help getframe
figure
for i=1:48
pcolor(x,y,squeeze(dev(i,:,:)));
shading interp;
caxis([-50 50])
colorbar
axis image
M(:,i)=getframe;
end
whos
% the movie is stored in the “struct array” M.
% we can now play it back:
movie(M)
% we can set up the number of times the movie is played and the speed in frames per second:
movie(M,2,3)
% now we have to display also information about the date
% so we now when we are during the movie.
% let’s see what the time format is:
format bank
t
% it is in decimal years, so it would be nice to extract the month information:
year=floor(t);
month=floor((t-year)*12)+1;
% we could put that information in the title, but there is a bug with the command
% movie, which does not refresh the title.
% so let’s put it in the figure itself, by using the command “text”
help text
text(240,50,’hello’)
% now the question is how to make the text change during the loop so that it displays
% the corresponding month and year ?
% we have to convert the numerical values of the month and year variables into characters:
i=1;
year(i) % is a number
num2str(year(i)) % is a character
text(240,40,num2str(year(i)))
% now we can display the year and month together by creating a character vector:
date=[num2str(month(i)),’/’,num2str(year(i))]
text(240,60,date)
% all is left to do is to plug these commands into the loop:
clf % clear figure
for i=1:48
pcolor(x,y,squeeze(dev(i,:,:)));
shading interp;
caxis([-50 50])
colorbar
axis image
date=[num2str(month(i)),’/’,num2str(year(i))];
text(240,50,date)
M(:,i)=getframe;
end
whos
movie(M,2,3)
% finally, what if we want to save our movie in our home directory to put it on our web site ?
% the best way is to print a series of figures in our home directory, and then to transform
% them into a movie, by using the linux command “convert”.
% the question is how to print each figure in different figure names without having to
% write the plot and print commands 48 times ?
% we have to enter the command we would write for each figure into a character vector,
% as we did previously for the date:
command=[‘print -djpeg fig’,num2str(i)]
% and then we have to tell matlab to execute this command, using the “eval” function:
% (first go to a directory where you have the right to write, i.e. not the one I go here,
% but your home directory for exemple)
cd /home/hokulea/ges/ocn363/students/jimp
eval(command)
end
ls
% you should see all the 48 different figures in your directory !
% note that there is a trick: I have added 10 to the variable “i” value, so that
% the figure names begin at 11 instead of 1.
% this is to keep the files in order when we list them in Linux,
% due to the way linux lists the files alphabetically.
% In class we tried to print them in jpeg format, but there seems to be a bug in matlab
% for printing figures in a loop in this format, so instead I have printed them in postscript here.
% Now quit matlab, go to the directory where you have saved the figures, and type
% in the terminal window:
% convert fig*.eps sshmovie.gif
% this will make an animated gif file, that you can visualize in your web browser !
% you can now erase all the figures (except the movie file !) to save disk space.
Opendap data
>> lon = ncread(‘http://apdrc.soest.hawaii.edu:80/dods/public_data/satellite_product/Rain-CMAP/clima/enh‘, ‘lon’);
>> lat = ncread(‘http://apdrc.soest.hawaii.edu:80/dods/public_data/satellite_product/Rain-CMAP/clima/enh‘, ‘lat’);
>> rain = ncread(‘http://apdrc.soest.hawaii.edu:80/dods/public_data/satellite_product/Rain-CMAP/clima/enh‘, ‘precip’);
Or, to save on typing…
>> URL = ‘http://apdrc.soest.hawaii.edu:80/dods/public_data/satellite_product/Rain-CMAP/clima/enh’;
>> lon = ncread(URL, ‘lon’);
>> lat = ncread(URL, ‘lat’);
>> rain = ncread(URL, ‘precip’);
>> rain_dec = squeeze(rain(:, :, 12));
>> figure()
>> pcolor(lon, lat, rain_dec’) % mean precipiation in december
>> colorbar
>>
Now lets play with a different data set and we will extract the long term mean and annual cycle
>> load /home/iniki/ges/ocn363-data/misc/avhrr_sst
>> time(1) % to see the beginning of the time set
>> I = find(lon > 201 & lon > 202);
>> J = find(lat > 20 & lat < 25);
>> ssh_hi = sst(:, I(1), J(1)); % select only the area around Hawaii
>> figure()
>> plot(time, ssh_hi)
>> X = ones(408,5);
>> w = 2.0 * pi / 12;
>> for K = 1:408 % this generates cycle with a combination of two cosines and two sines
>> wdt = w *K;
>> wdt2 = K * w * 2;
>> X(K,1) = 1.0;
>> X(K, 2) = cos(wdt);
>> X(K,4) = cos(wdt2);
>> X(K,5) = sin(wdt2);
>> end
>> figure()
>> plot(X(:, 1));
>> plot(X(:, 2));
>> plot(X(:, 3));
>> plot(X(:, 4));
>> [b, int] = regress(ssh_hi, X, 0.05); % this regresses the ssh to the annual-semiannual cycle created in the loop
>> ann_sst = X(:,1) * b(1) + X(:,2) * b(2) + X(:,3) * b(3) + X(:,4) * b(4) + X(:,5) * b(5);
>> figure()
>> plot(time, ssh_hi) % this is the raw time series of ssh
>> hold on
>> plot(time, ann_sst, ‘r’) % plotting the annual cycle on top of the raw time series
>> figure()
>> plot(time, ssh_hi - ann_sst) % this is the difference between raw time series and the mean annual-semiannual harmonic created in the for loop
script1:
load coast
xlon = long;
ylat = lat;
URL = ‘http://oos.soest.hawaii.edu/thredds/dodsC/hioos/satellite/dhw_5km‘;
%ncdisp(URL)
lat = ncread(URL,’latitude’);
lon = ncread(URL,’longitude’);
time = ncread(URL,’time’);
I = find ( lon > 133.5 & lon < 135.5 );
J = find ( lat > 6.5 & lat < 8.5 );
%datestr(double(time)/86400+datenum(1981,1,1,0,0,0))
K = find ( time == 86400 * (datenum(2015,7,15,12,0,0)-datenum(1981,1,1,0,0,0)))
sst = ncread(URL,’CRW_SST’,[I(1) J(1) K], [ length(I) length(J) K+30]);
anom = ncread(URL,’CRW_SSTANOMALY’,[I(1) J(1) K], [ length(I) length(J) K+30]);
dhw = ncread(URL,’CRW_DHW’,[I(1) J(1) K], [ length(I) length(J) K+30]);
pcolor(lon(I),lat(J),squeeze(anom(:,:,1))’)
shading interp
for L = 1:30
pcolor(lon(I),lat(J),squeeze(anom(:,:,L))’)
shading interp
caxis([-0.5 1])
colorbar
pause
end
for L = 1:30
pcolor(lon(I),lat(J),squeeze(anom(:,:,L))’)
date = datestr(double(time(K-1+L))/86400+datenum(1981,1,1,0,0,0))
shading interp
caxis([-0.5 1])
colorbar
text(133.7,8.3,date)
pause
end
script 2
clear all
load coast
xlon = long;
ylat = lat;
cmon = [‘Jan’; ‘Feb’; ‘Mar’; ‘Apr’; ‘May’; ‘Jun’; ‘Jul’; ‘Aug’; ‘Sep’; ‘Oct’; ‘Nov’; ‘Dec’ ];
load /home/iniki/ges/ocn363-data/misc/seawifs
%for i = 1: length(time)
for i = 1: 500
pcolor(lon,lat,log(squeeze(chl(i,:,:))’))
shading interp
caxis([-4 4])
colorbar
hold
plot(xlon,ylat,’k’)
t = datevec(time(i));
date = [ cmon(t(2),:), ‘ ‘, num2str(t(3)), ‘, ‘, num2str(t(1)) ];
text ( 102, 13, date )
M(:,i) = getframe;
% command = [‘print -dpng fig’, num2str(i+100)];
% eval(command)
hold
end
movie(M)
script3
%
% Plot climatologies
%
% load climatologies
load /home/iniki/ges/ocn363-data/clims/sst_clim
DJF = ( sst(12,:,:) + sst(1,:,:) + sst(2,:,:) ) / 3.0;
MAM = ( sst(3,:,:) + sst(4,:,:) + sst(5,:,:) ) / 3.0;
JJA = ( sst(6,:,:) + sst(7,:,:) + sst(8,:,:) ) / 3.0;
SON = ( sst(9,:,:) + sst(10,:,:) + sst(11,:,:) ) / 3.0;
subplot(2,2,1)
pcolor(lon,lat,squeeze(DJF)’)
shading interp
colorbar
subplot(2,2,2)
pcolor(lon,lat,squeeze(MAM)’)
shading interp
colorbar
subplot(2,2,3)
pcolor(lon,lat,squeeze(JJA)’)
shading interp
colorbar
subplot(2,2,4)
pcolor(lon,lat,squeeze(SON)’)
shading interp
colorbar
%load /home/hokulea/ges/ocn363/data/clims/ssh_clim
%A=squeeze(ssh(1,:,:));
%load /home/hokulea/ges/ocn363/data/clims/rain_clim
%A=squeeze(rain(1,:,:));
%load /home/hokulea/ges/ocn363/data/clims/ers_clim
%taux=squeeze(U(1,:,:));
%tauy=squeeze(V(1,:,:));
%quiver(lon,lat,taux’,tauy’)
We will continue our discussion around satellite data and expand our data analysis “toolkit” by looking at correlations. From wikipedia:
In statistics, dependence or association is any statistical relationship, whether causal or not, between two random variables or two sets of data. Correlation is any of a broad class of statistical relationships involving dependence, though in common usage it most often refers to the extent to which two variables have a linear relationship with each other.
Correlations are useful because they can indicate a predictive relationship that can be exploited in practice. For example, an electrical utility may produce less power on a mild day based on the correlation between electricity demand and weather. In this example there is a causal relationship, because extreme weather causes people to use more electricity for heating or cooling; however, correlation is not sufficient to demonstrate the presence of such a causal relationship (i.e., correlation does not imply causation).
Formally, dependence refers to any situation in which random variables do not satisfy a mathematical condition of probabilistic independence. In loose usage, correlation can refer to any departure of two or more random variables from independence, but technically it refers to any of several more specialized types of relationship between mean values. There are several correlation coefficients, often denoted ρ or r, measuring the degree of correlation. The most common of these is the Pearson correlation coefficient, which is sensitive only to a linear relationship between two variables (which may exist even if one is a nonlinear function of the other). Other correlation coefficients have been developed to be more robust than the Pearson correlation – that is, more sensitive to nonlinear relationships. Mutual information can also be applied to measure dependence between two variables.
The most familiar measure of dependence between two quantities is the Pearson product-moment correlation coefficient, or “Pearson’s correlation coefficient”, commonly called simply “the correlation coefficient”. It is obtained by dividing the covariance of the two variables by the product of their standard deviations.
The population correlation coefficient ρX,Y between two random variables X and Y with expected values μX and μY and standard deviations σX and σY is defined as:
ρX,Y = corr(X,Y) = cov(X,Y) / (σ:sub:X σY) = E[(X − μX) (Y − μY)] / (σ:sub:X σY)
where E is the expected value operator, cov means covariance, and corr is a widely used alternative notation for the correlation coefficient.
The Pearson correlation is +1 in the case of a perfect direct (increasing) linear relationship (correlation), −1 in the case of a perfect decreasing (inverse) linear relationship (anticorrelation), and some value in the open interval (−1, 1) in all other cases, indicating the degree of linear dependence between the variables. As it approaches zero there is less of a relationship (closer to uncorrelated). The closer the coefficient is to either −1 or 1, the stronger the correlation between the variables.
We will create a function in Matlab to compute the cross-correlation between two time-series. A function in Matlab is just a script, and it takes (optionally) input and output variables. For our function, we will want to supply two time-series (e.g., X and Y) and a range of lead/lags. For example, let’s say we want to correlate monthly mean temperature (SST) to sea level (SSH) for 10 month leads/lags, our function would be like corr(SST,SSH,10).
The return from the function should be the correlation coefficient as a function of lead/lag, and perhaps the lead/lag itself. So, a sample call to the function would be [r,lag] = corr(SST,SSH,10). Here is an example:
function [z,lag] = conn(x,y,n)% function [z, lag] = conn(x,y,n)% computes cross-correlation Cxy for lags -n to n%% convention x leads for +ve lag% y leads for -ve lag%% like smc/xcorfor.m except it does not plot.x = x-mean(x);y = y-mean(y);sx = std(x);sy = std(y);lx = length(x);i= 1;for lag = -n:-1;iy = 1:(lx+lag);
ix = (1-lag):lx;
z(i) = sum(x(ix).*y(iy))/length(ix);
i = i+1;
end
for lag = 0:n;
ix = 1:(lx-lag);
iy = (lag+1):lx;
z(i) = sum(x(ix).*y(iy))/length(ix);
i = i+1;
end
lag = -n:n;
z = z/(sx*sy);
Now we can run it on some data:
>> clear all>> load /home/iniki/ges/ocn363-data/misc/avhrr_sst>> J = find ( lat > 21.4 & lat < 21.8 );>> I = find ( lon > 201 & lon < 202 );>> sst_hi = sst(:,I,J);>> load /home/iniki/ges/ocn363-data/misc/CMAP_precip>> J = find ( lat == 21.25 );>> I = find ( lon == 201.25 );>> rain_hi = squeeze(rain(:,I,J));>> [z,lag] = xcorrel(rain_hi,sst_hi,12);>> plot(lag,z)
I. Numerical models¶
Thus far we have looked at a variety of observational data, both in situ and remote. Each has its own benefits and limitations. Some of the limitations include the necessary restrictions on data availability. For example, for the most part satellites can only “see” the surface. In situ measurements can be made of the sub-surface, but these are limited in geographic extent. Another constraint is that it is difficult to perform multiple experiments in situ. To overcome these, many people use numerical models. Numerical models are approximations of physical systems and use numerical methods to solve geophysical equations. There are many different types of models, including ocean-only models that require inputs of surface forcing (e.g., surface winds and heat fluxes), atmospheric models that require information at the land/ocean surface (e.g., surface temperature), and coupled models that are entirely “self contained”. For this class we will look at coupled model results from the IPCC AR-4 experiments. A few examples have been placed in our data directory, /home/iniki/ges/ocn363-data/IPCC:
With these runs, just a very small subset, we can do things like model inter-comparisons, scenario comparisons, etc. First, let’s just make some plots. Note that these model output are in netCDF format (see notes from earlier this semester, so we need to use “ncdisp” and “ncread” to access them. % check out the contents of an example file ncdisp(‘/home/iniki/ges/ocn363-data/IPCC/SresA1FI_Z1_atmos.2001-2100.precip.nc’) % notice that we have variables lon (144), lat (90), time (100), and precip (lon,lat,time) [kg/m2/s] % let’s save these three (recall the syntax on ncread is file name, variable name, index start, number, stride. % here we read the whole variable, so no need the last part (index start, number, stride). First, since we % may need it, let’s load the coastline file
|
|---|