% ld.m
% script to do the following:
% 0. read in the binary SST data file
% 1. apply latitude weights prior to covariance estimation
% 2. remove land locations from the global grid
% 3. reshape the data into lon,lat,time data structure
% 4. plot a selected time window 
%
% data source: Bottomley et al. (1990)
% retrieved from markov.ldeo.columbia.edu
% with the following grid structure:
% Time
%    grid: /T (months since 1960-01-01) ordered Jan 1961 to Dec 1990 by 1. N= 360 pts :grid 
%Longitude
%    grid: /X (degree_east) periodic 177.5W to 177.5E by 5. N= 72 pts :grid 
%Latitude
%    grid: /Y (degree_north) ordered 87.5S to 87.5N by 5. N= 36 pts :grid 
%
% missing_value: NaN

clear

% useful definitions
% by setting these values in one place, I can change them
% globally here as well, should I apply this script
% to a different dataset.
nlat=36; % number of latitude gridboxes per time slice2
nlon=72; % number of longitude gridboxes per time slice
nxny=36*72; % number of spatial points in each timeseries
ny=30; % number of years in timeseries
nm=30*12; % number of months in timeseries
nlats=-87.5;
nlatn=87.5;
nlonw=-177.5;
nlone=177.5;
gridstep=5;

% 0. load the binary data
% the following commands read in the binary dataset 
[fid,m1]=fopen('m6190.r4','r','ieee-be');
[sst,c]=fread(fid,[nlon*nlat,nm],'float');
[nxny,nt]=size(sst)
% rows are points in space on a 5x5 grid;
% longitude varies the fastest, then latitude
% columns are time (jan 1961-dec 1990)

% 1. apply area weights
% grid area is proportional to cos(lat)
% consider the distortion in high latitudes 
% in a mercator projection
% we need to weight the data appropriately
% prior to EOF analysis
% otherwise we give inordinate weight to relatively
% small areas in the covariance and EOF calculations

latwt=cos(pi*[-87.5:5:87.5]./180);
lw=diag(latwt)*ones(nlat,nlon);
lwnxny=reshape(lw',1,nxny);
% apply weights
for i=1:nm
sstw(:,i)=sst(:,i).*lwnxny';
end
% after analysis, we'll reverse this weighting to plot the field in space:
% e.g. sstuw(:,i)=sstw(:,i)./lwnxny;

% 2. Remove the land area gridpoints to get the marine data locations:
[fid,m1]=fopen('m5mask.r4','r','ieee-be');
[land,c]=fread(fid,[nlon*nlat,1],'float');
[nxny]=length(land)
     ogrid=find(~isnan(land));
     sstwo=sstw(ogrid,:);
% put back to grid: data(nxny,nm): data(ogrid,n)=sstwo(:,n),n=1,nm

% 3. Put the data into (lon,lat,monthtime) 3-D matrix form:
sstf=reshape(sstw,nlon,nlat,nm); % see help for reshape

% 4. Plot using m_pcolor
% note the use of transpose and flipud to get data in form to plot
% using m_pcolor
% see Rich Pawlowicz' m_map routines for more information
% http://www2.ocgy.ubc.ca/~rich/private/mapug.html
% make sure the m_map functions are in your matlab path
addpath /usr/share/MATLAB/toolbox/m_map

% the next four lines will give you a map of a given month's data:
% for instance, let test= Jan 1961:
test=squeeze(sstf(:,:,1));  % see help squeeze

 [lon,lat]=meshgrid([nlonw:gridstep:nlone],nlats:gridstep:nlatn);
 m_proj('gall-peters','lon',[-180 180],'lat',[-90 90])
 m_pcolor(lon,lat,flipud(test')); 
 m_coast;
 m_grid;


save lddat.mat