Decorrelation Scales of Sea Surface Temperature
Imagery in the Gulf of Maine
Final Report to the
Maine Space Grant Consortium Fellowship Program
September 1, 1999
Peter Brickley
Andrew Thomas
School of Marine Sciences
University of Maine
Orono, ME 04469
peter@grampus.umeoce.maine.edu
Abstract
The purpose of this research was to investigate the regional differences in decorrelation spatial scales of sea surface temperature in the Gulf of Maine, which reflect variability in ocean physics. Dominant patterns of spatial scales were extracted from three years of high resolution sea surface temperature images using robust statistical measurement of two-dimensional autocorrelation. Several major regional variations in spatial scales, and preferred orientation of the correlation ellipses were found, segregating major circulation features in the Gulf of Maine. Comparison between several years showed that mean spatial scale maps exhibit interannual variation. Estimates of the spatial scales of SST pattern evolution, their degree of regional variation, and statistical significance will provide more confidence in using satellite-observed SST as a diagnostic tool, and in turn, suggest realistic spatial scales for properties which are less easily sensed remotely.
Introduction
Satellite images of sea surface temperature (SST) in the Gulf of Maine (GOM) commonly exhibit dramatic temperature changes over spatial scales of a few km associated with both mesoscale or gulf-scale circulation features. Examples of this include a plume of relatively cold SST erupting from the eastern Maine coastal current (EMCC), "streamers" in the terminal end of EMCC, and cold SST bands between the Jordan and Georges Basin Gyres. These thermal gradients have a finite spatial extent and temporal persistence. To date, efforts to quantify this variability using shipboard transects and moored sensors have suffered from an acknowledged degree of under sampling which restricts our ability to interpret these features in their proper oceanographic context. The purpose of this research was to investigate regional differences in decorrelation spatial scales of sea surface temperature in the Gulf of Maine. These differences indicate fundamental spatial scales of variability in the underlying ocean physics.
Several studies suggest that spatial scales vary between sub-regions of the Gulf of Maine, and in particular between eastern GOM and western GOM. Such a east-west "split" was hypothesized by Bigelow (1927) and recent estimates of air-sea heat flux (Mountain, et al, 1996) and observation of thermal fronts (Ullman and Cornillon, 1999) support this idea. Several studies (e.g. Vermersch, et al. (1989) in the western GOM; Hetland (1996) in the eastern coastal GOM, and Brickley (1994) in Massachusetts Bays) have surmised the existence of relatively energetic eddy fields from observations of low coherence between currents and pressure observations separated by distances as little as 10-30 km. Little is known about the hydrographic heterogeneity of the Gulf’s interior and extent to which this effects SST fields. Complex SST features, such as streamers, eddies, fronts, and other small scale, large gradient SST variations are frequently noted in SST imagery. However their small spatial and temporal scales result in their absence from even short term averages of SST. Here we address the issue of spatial variability. By calculating decorrelation scales from individual images and then averaging these, information about extent and position of small scale features is preserved. Robust statistical measurements of dominant patterns of spatial scales can be extracted from high-resolution sea surface temperature images using techniques of cross-correlation and spectral analysis. In addition the orientation of these spatial correlation scales provides an estimate of "variance isotropy", or dominant direction of hydrographic variability.
Estimates of the spatial scales of SST pattern evolution, their degree of regional variation, and statistical significance will provide more confidence in using satellite-observed SST as a diagnostic tool, and in turn, suggest realistic spatial scales for other surface oceanographic properties which cannot be measured remotely. This may also guide coordination of future sampling efforts to capture synoptic features by suggesting spatial limits on how often and at what scales to sample to resolve spatial variability in the region. Quantifying the dominant spatial scales of patterns in SST data will assist interpretation of satellite data, in situ moored and ship-borne measurements. This is especially critical to improving our understanding of circulation features when it is not possible to collect data at sea.
Data and Methods
Image database
The image database used was the NOAA/AVHRR Pathfinder image catalog from 1990-1994. This data set was derived from the 5 channel Advanced Very High Resolution Radiometer (AVHRR) carried on board the NOAA satellites. The AVHRR measures radiance in five spectral bands ranging from 0.55 to 12.5 m m. Raw data are converted to a surface temperature (in degrees Celsius) using a multi-channel sea surface temperature algorithm (MCSST) algorithm and the subjected to the Pathfinder processing (Smith et al., 1996; McClain, et al., 1985). Thermal precision, or resolution, in satellite derived surface temperature from these images is on the order of 0.15° C. Each satellite pass has a spatial resolution of 1.2 km at nadir. A composite daily image was formed from several images for each day. During development of the processing algorithms this archive was supplemented by nearly 3 years of SST imagery from the School of Marine Science TeraScan satellite receiving station (Sea Space Corp., San Diego, Ca.). Both the Univ. of Maine and the NOAA Pathfinder images were navigated to the same region at an accuracy of approximately one pixel and processed to remove cloud features. Image sizes were 604 x 684 pixels with lower left-hand and upper right-hand corners located at approximately (-72, 38.5) and (-64,45), respectively. This region included large portions of the Middle Atlantic Bight and Gulf Stream. Calculations used a 400 x 400 pixel sub-image to reduce computation time.
GOM summer images from June-August of three consecutive years (1990-1992) were subjected to spatial scale analysis. The multi-year array of images was inspected to assess quality. We chose to focus on the summer period (June-August) when thermal signatures of regional circulation features are enhanced by solar heating and seasonal stratification. The higher temperature contrast also emphasize regional variations across the GOM. Because of prevailingly mild weather patterns, a high percentage of cloud-free images were obtained during summer months of these three years. This resulted in approximately 8 nearly cloud-free images per month and approximately 12 images per month where important sub-regions within the GOM were visible. Each image was subjectively rated on a scale of 1-10 (10 being cloud free) and grouped as "acceptable" if the subjective rating was greater than 6 (Table 1).
The results of Vermersch, et al. (1989) and Hetland (1996) suggest that spatial decoherence of current fields occurs rapidly between 10 and 30 km. This may be reflected in spatial decorrelation of surface temperature patterns. In all calculations a 30 x 30 pixel subimage is the basic computational element of the correlation analysis. This sub-image size is a compromise between excessive computation time required by an incremental number of smaller areas and larger dimensions which result in too few spatial points to resolve spatial pattern. A sub-region overlap between adjacent regions of 50% (15 pixels) was used to increase spatial resolution.
Results
1-dimensional algorithm
The autocorrelation function is defined as
(1)
where k is lag, k = 0,1,2,...,m (m<N) and where the lag .
In practice the procedure is to calculate the autocovariance function by subtracting a mean value. Failure to subtract the mean can change the lag value at which the correlation approaches zero. The autocovariance is defined as
(2)
where and is the lag time for k sampling increments (Emery and Thomson, 1997). It is seen that the autocovariance function is symmetric with respect to lag. Division by (N-k) for each lag k gives the unbiased estimate, however, in practice this function is often normalized by the maximum zero lag covariance. Normalization results in an autocorrelation function with limits between 1 and -1. The correlation function value of zero represents the maximum lag beyond which points are no longer statistically dependent, and thus determines the "decorrelation scale". Computing an autocorrelation function for various lags determines the dominant spatial scales in SST variability within a selected range of lags.
2-dimension spatial correlation algorithm
Extension of the one dimension covariance function to two dimensions is straightforward and follows by extending equation (2) to two dimensions. The resulting two-dimensional autocorrelation function is sensitive to feature shape, size and orientation and provides a measure decorrelation length scale of image features and their orientation
Algorithms for this method are included in Appendix A. The algorithm implements the autocorrelation calculation for 30 x 30 sub-images of SST. A least-squares planar trend is removed from each 30x30 (or nxn) sub-image prior to calculation of the autocovariance to remove influence trends in SST at scales larger than the sub-image. The autocovariance is calculated from a two-dimensional fast Fourier transform (FFT), except where bad pixel data are flagged in the sub-image. In the latter case, a two-dimensional summation routine is used (see Appendix A, correl2.m and xcorr_s.m). The resulting autocovariance function is normalized by the maximum zero-lag covariance to yield the autocorrelation function. An example is shown in Figures 1-3 of the autocorrelation (the autocovariane normalized by the maximum zero-lag autocovariance). Figure 1 shows a 22 x 22 pixel sub-image of July SST from the central Gulf of Maine. Figure 2 shows a surface plot of the two-dimensional normalized autocorrelation function. The correlation function zero level frames an ellipse about the central maximum. Figure 3 shows that the orthogonal major and minor axis length scales of the correlation ellipse have a length scale of approximately 12 and 4.5 pixels (zero crossing of the autocorrelation).
Correlation ellipse parameters are extracted from the normalized covariance output by contouring the value of zero correlation. Parameters of major and minor axis length, ellipse orientation, and ellipticity are measured from the origin to the zero correlation contour. The length scale of maximum correlation is defined by the major axis length a and is oriented in a xy plane at an angle defined by , where a and b are the major and minor axes lengths. Ellipticity is defined and is used as a measure of isotropy. Note that is 1 if b=0, and 0 if a=b. Sub-image features that do not have a dominant direction of orientation have high isotropy and =0.
Figure 4 is July 21, 1994 SST image of the Gulf of Maine. The boxed area is a 50 x 50 pixel image of a thermal front over Truxton Swell between cold Scotian Shelf water and warmer waters over Jordan Basin. Figure 5 shows this sub-image with overlying vectors indicating the major axis length scale and orientation of the autocorrelation ellipse (called a correlation map). Each arrow is a vector representation of the maximum decorrelation length scale. The vector originates at the center of each sub-image and the arrowhead has no significance, except to indicate the orientation of correlation ellipse. In the region of the frontal boundary, the major axis length scales increase (reaching about 22 km) and become aligned with the front. Far from the front, the orientation of the ellipse major axis is less organized. The mean ellipticity also decreases, indicating greater isotropy in thermal features (not shown).
Figure 6 shows a Gulf wide SST image from July 9, 1991 (yearday 190). Figures 7 and 8 show the corresponding map with vectors of the ellipse major axis length and orientation (Figure 7) and the major axis length (Figure 8). In spite of its complexity, close inspection reveals several regions with common length scales and orientation. For example, major axis length scales increase along the thermal boundary between the eastern and western Gulf and tend or orient along this feature. In the eastern Gulf of Maine length scales are generally longer (Figure 8) and average between 22 and 30 pixels in length. Longer spatial scales are reflective of more coherent, spatially persistent physical processes modulating SST pattern. To the west, major axis length scales show greater spatial variability across subimages, a general decrease in length, and although this is not shown, have greater isotropy (lower ellipticity). Spatial scales are also extended along the frontal region around Georges Bank.
An average of all available subimages for June, July, and August were compiled into monthly correlation maps. Missing or bad values were excluded in the mean. An example of these results is shown in Figures 9 and 10. Figure 9 shows the major and minor axes lengths of the mean correlation ellipses for July 1991. A reference scale appears in the upper right corner. Figure 10 presents the number of subimage correlation values used to form the mean correlation map (N). Several dominant regional variations in spatial scales are seen. In particular, a shift in spatial scale and orientation along a curving boundary between eastern and western GOM is apparent. This appears in the mean correlation maps of June-August in all three years. This region is coincident with the dominant pattern of SST evident in most summer imagery of the interior Gulf of Maine (as can be seen in Figure 6). The calculations performed here quantify the along and cross length scales of this feature. Recent work by Ullman and Cornillon (1999) shows the presence of a Gulf-scale frontal system that extends in summer across the Gulf of Maine from the shelf west of Penobscot Bay to the northern flank of Georges Bank. The broad area shown in Figure 9 where spatial scales become elongated recurs with great fidelity along a similar corridor crossing the central Gulf of Maine, and appears to confirm these observations. The frontal region around Georges Bank shows a thin but permanent region of increased spatial scale, oriented parallel to the front and averaging 16-25 pixels in length. The similarity of Figure 9 to Figure 6 indicates semi-permanence to regional scales throughout the month of July. Comparison to other years shows that these mean spatial scale maps exhibit some interannual variation.
Conclusions
This research is an attempt to limit guesswork when interpreting spatial variability in continental shelf SST imagery. This research uses a simple autocorrelation technique to quantify the spatial variation of length scales and dominant direction of variability without discarding fine-scale features. Tests of this method on several years of summer SST images have shown that regional variation in correlation ellipse orientation and major axis length is complex, but shows a general fidelity to major Gulf-scale SST features. Future improvements of this technique should seek to establish confidence limits on derived spatial scales, extend the image database across more seasons and years, and compare results with numerical model simulations of SST and in situ moored data.
References
Bigelow, H. B., Physical Oceanography of the Gulf of Maine, Department of Commerce, Washington, DC, 1927
Bisagni, J. J., Gifford, D. J. and C. M. Ruhsam, The spatial and temporal distribution of the Maine Coastal Current during 1982, Cont. Shelf. Res. 16, 1, 1996.
Bisagni, J. J., R. C. Beardsley, C. M. Ruhsam, J. P. Manning, and M. J. Williams, Historical and recent evidence of Scotian Shelf water on Georges Bank, Deep Sea Res., 43, 1439, 1996.
Brickley, P. J., Wind stress and subtidal ciruculation in Massachusetts and Cape Cod Bay, M.S. thesis, University of Maine, Orono, Me, 130 pg., 1994.
Brooks, D. A., The influence of warm-core rings on slope water entering the Gulf of Maine, J. Geophys. Res., 92 (C8), 8183, 1987.
Chelton, D. B. and M. G. Schlax, Estimation of time-averaged chlorophyll concentrations from irregularly-spaced satellite observations, J. Geophys. Res., 95, 1990.
Denman, K. L. and M. R. Abbott, Time scales of pattern evolution from cross-spectrum analysis of AVHRR and CZCS imagery, J. Geophys. Res., 99 (C4), 7433, 1994.
Deschamps, P. Y., Frouin, R. and L. Wald, Satellite determination of the mesoscale variability of the sea surface temperature, J. Phys. Oceanogr., 11, 864, 1981.
Emery, W. J. and M. Ikeda, A comparison of geometric correction methods for AVHRR imagery, Can. J. Rem. Sens., 10,1, 46, 1984.
Emery, W. J., and R. E. Thomson, Data analysis methods in physical oceanography. 1st ed. Elsevier Science, New York, New York, 634 pp., 1997.
Hetland, R. D. The evolution of the vernal circulation in the eastern Gulf of Maine, M.S. thesis, University of Maine, Orono, Me, 98 pg., 1997.
McClain, E. P., W. G. Pichel, and C. C. Walton, Comparative performance of AVHRR-based sea surface temperatures, J. Geophys. Res., 90, 11,587-11,601, 1985.
Mountain, D. G., G. A. Strout, and R. C. Beardsley, Surface heat flux in the Gulf of Maine, Deep-Sea Res. II, 7-8, 1533-1546, 1996.
Pettigrew, N. R. and R. D. Hetland, The cyclonic gyre of the eastern Gulf of Maine. Abstract GLO-5 Proceedings of the Oceanography Society, 1995.
Pettigrew, N. R., D. W. Townsend, H. Xue, J. P. Wallinga, P. J. Brickley and R. D. Hetland. Observation of the Eastern Maine Coastal Current and its offshore extensions in 1994. J. Geophys. Res. 103 (C13), 30, 623-30,639, 1998.
Townsend, D. W. and N. R. Pettigrew, The role of frontal currents in larval fish transport on Georges Bank, Deep Sea Res., 43, 1773, 1996.
Ullman, D. A, and P. C. Cornillon, Satellite-derived sea surface temperature fronts on the continental shelf off the northeast U.S. coast, J. Geophys. Res., 104 (C10), 23,459-23,478, 1999.
Vermersch, J. A., Beardsley, R. C. and W. S. Brown, Winter circulation in the western Gulf of Maine. Part 2: Current and pressure observations. J. Phys. Oceanogr. 9, 768, 1979.
Tables
Table 1
Subjective image quality estimation of GOM MCSST images for 1990-1992. The correlation scale is computed for each 30 x 30 pixel sub-image if N_{good} > 70%. Missing column entries indicate a score less than 5.
Yearday |
Month |
Q-1990 |
Q-1991 |
Q-1992 |
Yearday |
Month |
Q-1990 |
Q-1991 |
Q-1992 |
152 |
June |
8 |
8 |
198 |
8 |
5 |
6 |
||
153 |
8 |
199 |
8 |
8 |
6 |
||||
154 |
6 |
5 |
200 |
7 |
|||||
155 |
6 |
5 |
201 |
5 |
|||||
156 |
8 |
7 |
202 |
5 |
|||||
157 |
8 |
203 |
6 |
||||||
158 |
9 |
204 |
|||||||
159 |
9 |
205 |
7 |
||||||
160 |
8 |
206 |
5 |
||||||
161 |
6 |
207 |
6 |
||||||
162 |
8 |
208 |
5 |
||||||
163 |
5 |
7 |
209 |
7 |
|||||
164 |
6 |
8 |
210 |
7 |
|||||
165 |
7 |
5 |
7 |
211 |
9 |
||||
166 |
8 |
6 |
212 |
||||||
167 |
5 |
213 |
August |
8 |
8 |
||||
168 |
5 |
214 |
5 |
8 |
8 |
||||
169 |
6 |
215 |
7 |
9 |
|||||
170 |
216 |
8 |
8 |
||||||
171 |
5 |
217 |
8 |
||||||
172 |
6 |
218 |
6 |
8 |
5 |
||||
173 |
7 |
219 |
9 |
||||||
174 |
220 |
6 |
9 |
||||||
175 |
7 |
5 |
221 |
5 |
6 |
||||
176 |
8 |
222 |
5 |
||||||
177 |
8 |
223 |
6 |
7 |
|||||
178 |
5 |
5 |
224 |
8 |
|||||
179 |
225 |
6 |
9 |
||||||
180 |
7 |
226 |
7 |
||||||
181 |
9 |
227 |
8 |
||||||
182 |
July |
8 |
228 |
8 |
|||||
183 |
8 |
7 |
229 |
8 |
6 |
||||
184 |
230 |
8 |
|||||||
185 |
8 |
5 |
231 |
||||||
186 |
8 |
232 |
|||||||
187 |
8 |
233 |
7 |
||||||
188 |
8 |
5 |
234 |
9 |
|||||
189 |
8 |
235 |
8 |
8 |
|||||
190 |
9 |
8 |
236 |
6 |
8 |
||||
191 |
7 |
9 |
237 |
8 |
7 |
||||
192 |
9 |
238 |
6 |
5 |
|||||
193 |
7 |
8 |
239 |
6 |
7 |
||||
194 |
7 |
240 |
8 |
||||||
195 |
241 |
||||||||
196 |
7 |
6 |
242 |
6 |
|||||
197 |
6 |
243 |
9 |
9 |
Figures
Figure 1. 22 x 22 pixel SST sub-image with a narrow, elongated thermal feature
Figure 2. Surface plot of the two-dimensional normalized autocovariance function. The covariance function zero level frames an ellipse about the central maximum.
Figure 3. Orthogonal major and minor axis length scales of the image in Figure 1.
Figure 4. July 21, 1994 SST image of the Gulf of Maine. A 50 x 50 pixel box indicates the sub-image used in the example autocovariance map.
Figure 5. Vectors indicating the major axis length scale and orientation of the autocorrelation ellipses for the image in Figure 4. Arrows originate on the sub-image center and indicate the orientation of the ellipse, not direction of flow.
Figure 6. SST image from July 9, 1991 (yearday 190).
Figure 7. Map of major and minor axis length scales of the autocorrelation ellipses for the image in Figure 6. Arrowheads indicate the orientation of the ellipse not direction of flow
Figure 8. Map of major axis length scales of the autocovariance ellipses for the image in Figure 6.
Figure 9. Mean autocorrelation map for the month of July, 1991. The markers indicate the major and minor ellipse axes lengths and the orientation of each ellipse. A reference length scale is in the upper right hand corner. The ellipse scales are plotted over the July 9, 1991 image (from Figure 6) to give a sense of location.
Figure 10. Number of values used to compute the mean autocorrelation ellipse length scales in Figure 9. Sub-image autocorrelation was computed for subimages with fewer than 30% bad pixels in any composite daily image.
Appendix A
MALTA Functions
function [xl,yl,xi,yi,eci,eca,ecb]=corr_compile(direc,ydlim)
%function [xl,yl,xi,yi,eci,eca,ecb]=corr_compile(direc,ydlim)
% compile stack of correlation scale and isotropy
% in m,n,p data matrix for each variable (xi;yi;eci)
% direc - directory of files (string)
% ydlim - yearday limits for filename compilation
% decompress and load in each image file
[lst]=corr_list(direc,ydlim) % compile the list of filenames
% main loop
% first uncompress each file and then pass it into im_corr2
for ii=1:length(lst(:,1))
eval(['!uncompress ',lst(ii,:)]);
I=im_path(lst(ii,1:11));
[xl,yl,xi(:,:,ii),yi(:,:,ii),eci(:,:,ii),eca(:,:,ii),ecb(:,:,ii)]=im_corr2(I(50:450,150:550)); % insert into 3d array
h=corr_plot(lst(ii,1:11),I,xl,yl,xi(:,:,ii),yi(:,:,ii),eci(:,:,ii),eca(:,:,ii),ecb(:,:,ii)); % output
% determine which are good images
eval(['cd ',direc]);
eval(['!compress ',lst(ii,1:11)]);
end;
function I=im_path(filename)
%function I=im_path(filename)
%function to fread a Pathfinder MCSST image given a input filename
% calls the following
% fopen ('filename','r','b')
% fread(fid,[684,604],'uint8')
if (isstr(filename)==0) then
disp ('Error, must be a string'), break,
end;
a=pwd
fid=fopen(filename,'r','b')
I=fread(fid,[684,604],'uint8');
I=I'.*0.125 - 2.0;
fclose(fid);
function [xl,yl,xi,yi]=im_corr2(im)
%function [xl,yl,xi,yi]=im_corr2(im)
% auto corrlation of image
% set boundary pad to sub-array size
% nb=mb=19
[m,n]=size(im); % image dim
nb=9;mb=9; % boundary
ns=m-(2*mb+1); % sqrt(number of samples)
% preallocate arrays
i_e_ij=zeros(2,2,ns^2);
xy_ij=zeros(1,2,ns); % allocate 3 d arrays
% generate indices
isec=[[1:ns]' [1:ns]'+(2*mb+1)]
jsec=isec;
pause
k_1=1;
for ii=1:ns
for jj=1:ns
i_e_ij(:,:,k_1)=[[isec(ii,:);jsec(jj,:)]];
k_1=k_1+1;
end;
end;
% main loop
for kk=1:ns^2
kk
[mu,mv]=dist_orient(correl2(im(i_e_ij(1,1,kk):i_e_ij(1,2,kk),i_e_ij(2,1,kk):i_e_ij(2,2,kk))));
xy_ij(:,:,kk)=[mu mv];
end;
xi=[reshape(xy_ij(1,1,:),ns,ns)]'; % (ns x ns) matrix of x scales
yi=[reshape(xy_ij(1,2,:),ns,ns)]'; % (ns x ns) of y scales
% offset axis
xl=[1:ns]-(2*mb-1);
yl=[1:ns]-(2*nb-1);
function [mu,mv,ec,eca,ecb]=dist_orient(arr)
%function [mu,mv,ec,eca,ecb]=dist_orient(arr)
% arr - input 2d autocorrelation function
% mu,mv - output is x,y coordinates of maximum zero level contour
% ec is eccentricity
% eca, ecb are ellipse major and minor axis lengths
lev=0.05;
c=contourc(arr,[lev lev]); %contour 0.05 level
[mm,nn]=size(arr);
m1=length(c(1,:));
%form vectors, first search segments for the longest (central)
a1=find(c(1,:)==lev);
n_seg=length(a1);
[q,i]=sort(c(2,a1)); % search for longest segment
i_len=a1(i(n_seg)); % index of first
s_len=q(n_seg); % segment length
u=c(1,i_len+1:i_len+s_len);
v=c(2,i_len+1:i_len+s_len);
w=u+sqrt(-1).*v; % omplex coords
%
maxw=max(demean(w)); %tricky multiple contour lines
minw=min(demean(w));
theta=180./pi.*angle(maxw);
mu=real(maxw);
mv=imag(maxw);
ec=sqrt((abs(maxw)).^2 - (abs(minw)).^2)./abs(maxw); % eccentricity
eca=abs(maxw);
ecb=abs(minw);
%ec2=(abs(minw))./(abs(maxw)) % ratio of min to max axis
if (mu > mm)
mu=mm;
end;
if (abs(mv) > mm)
mv=mm;
end;
%mv=sign(theta).*mv;
% orientation in cartesian coords is
%aa=2.*mean(demean(u).*demean(v));
%bb=mean(demean(u).^2)-mean(demean(v).^2); % 0 deg = (x,y=0)
%theta=atan2(aa,bb).*180./(2.*pi);
function c_norm=correl2(arr)
%function c_norm=correl2(arr)
% input array and output autocorrelation function estimate
% farray - nxn array
% c_norm - m x m autocorrelation function
% m = 2n-1
% implements both summation and fft autocorrelation routines
% normalized by zeroth lag autocorrelatin value Rxx(0), or mean square.
% removes the mean (if not removed already), so really a autocovariance
% function
% 90% significance level depends on number of degrees of freedom
% Use mean autocorrelation length scale. edf=(n*n)./length
% see Bendat and Piersol
%
% 6/10/98 pb
tol=10*eps;
[m,n]=size(arr);
%flag any nan values
nflg=any(any(isnan(arr)));
%replace bad vals with mean values for purposes of fitting a trend
if nflg
ind=find(isnan(arr));
m_arr=mean(mean(arr(~isnan(arr)))); %replace with denaned mean val
arr(ind)=mean(mean(arr(~isnan(arr))));
end;
% remove a linear trend (a best-fit plane)
% form regression matrix
X=[ones(m,1) [1:m]'];
a1=X\mean(arr)'; % solve equation set.
a2=X\mean(arr')';
if (a1(2) <tol | a2(2) < tol) % machine precision error
a1(2)=0;
a2(2)=0;
end;
B1=ones(m,1)*(X*a1)';
B2=ones(m,1)*(X*a2)';
C=(B1+B2')./2;
ARR=arr-C;
clear B1 B2 a1 a2 arr C X
% remove mean, averaging over good points
p=ARR-mean(gmean(ARR));
%raw autocorrelation
if nflg
c_arr=xcorr2_s(p,p); % via summation
else
c_arr=xcorr2(p,p); % via fft
end;
% normalize by maximum zero lag correlation Rxx(0)
if max(max(c_arr))==0
c_norm=c_arr;
else,
c_norm=c_arr./max(max(c_arr));
end;
function r=xcorr2_s(a,b)
%function r=xcorr2_s(a,b)
% 2d cross-correlation routine, using summation,
% rather than FFT. Handles NaN's without propagation
% Returns unnormalized cross-correlation function
% 6/24/98 pb
%
if nargin == 1
b = a;
end
[ma,na] = size(a);
[mb,nb] = size(b);
b = conj(b(mb:-1:1,:));
apad = [a ; zeros(mb-1,na)];
r = zeros(ma+mb-1,na+nb-1);
for k=1:(na+nb-1)
count = k*(k<min(na,nb))+min(na,nb)*(k>=min(na,nb))*(k<=max(na,nb))+(na+nb-k)*(k>max(na,nb));
starta = 1*(k<=nb)+(k-nb+1)*(k>nb);
startb = (nb-k+1)*(k<=nb)+1*(k>nb);
for i=0:(count-1)
b1=b(:,startb+i);
a1=apad(:,starta+i);
na1=length(apad(:,starta+i));
nb1=length(b(:,startb+i));
for kk=1:(na+nb-1)
indb=fliplr([1:kk-((kk-nb1)*(kk>nb1))]);
inda=[1+((kk-nb1)*(kk>nb1)):kk*(kk<=nb1)+kk*(kk>nb1)];
P=a1(inda).*b1(indb);
r(kk,k) = r(kk,k)+sum(P(~isnan(P)));
end;
end;
end;
function c = xcorr2(a,b)
%XCORR2 Two-dimensional cross-correlation.
% XCORR2(A,B) computes the crosscorrelation of matrices A and B.
if nargin == 1
b = a;
end
[ma,na] = size(a);
[mb,nb] = size(b);
b = conj(b(mb:-1:1,:));
apad = [a ; zeros(mb-1,na)];
c = zeros(ma+mb-1,na+nb-1);
for k=1:(na+nb-1)
count = k *(k<min(na,nb)) ...
+min(na,nb) *(k>=min(na,nb))*(k<=max(na,nb)) ...
+(na+nb-k) *(k>max(na,nb));
starta = 1 *(k<=nb) ...
+(k-nb+1) *(k>nb);
startb = (nb-k+1) *(k<=nb) ...
+1 *(k>nb);
for i=0:(count-1)
c(:,k) = c(:,k) + filter( b(:,startb+i), 1, apad(:,starta+i) );
end
end
function h=corr_plot(lst,I,xl,yl,xi,yi,eci,eca,ecb);
%function h=corr_plot(lst,I,xl,yl,xi,yi,eci,eca,ecb);
% function to display pathfinder images overlaid with
% correlation scale vectors and indices of isotropy
% plot each image, overlay the major and minor axes.
% lst - image filename
% I - image variable
% xl, yl - scales
% xi, yi, - coordinates of major axes vector
% eca, ecb - major and minor axes scales
% eci - eccentricity of correlation ellipses
% complementary set of minor axes coordinates
theta=(180./pi).*asin(yi./eca);
x2i=ecb.*cos((theta+90).*pi./180);
y2i=ecb.*sin((theta+90).*pi./180);
%plot image array
figure(1);clf;
h=imagesc([1:400],[1:300],I(50:450,150:550));
%h=imagesc(I(50:150,150:250));
hold on
% scale axes by 1/2 grid distance
p=diff(xl(1:2))./2; % 1/2 grid dimension
q1=quiver(xl,yl,xi./p,yi./p,0,'w-');
q1a=quiver(xl,yl,-xi./p,-yi./p,0,'w-');
q2=quiver(xl,yl,x2i./p,y2i./p,0,'w-');
q2a=quiver(xl,yl,-x2i./p,-y2i./p,0,'w-');
xlabel('Distance [pixels]');
ylabel('Distance [pixels]');
title(['Daily composite SST: ',lst])
axis xy
axis('square')
posit=get(gca,'position')
cb=colorbar;
cb_posit=get(cb,'position');
set(cb,'position',[cb_posit(1) posit(2) cb_posit(3)-0.025 posit(4)]);
axes(cb);
ylabel('SST [\circC]')