lagmoran

 PURPOSE: calculate Moran's I for a lag in x and/or y direction on
          a regular grid. 
 -------------------------------------------------------------------
 USAGE: m = lagmoran(raster, xlag, ylag);
 where: [raster] is the grid to calculated Moran's I for
        [xlag],[ylag] is the lag in x and y directions. Postive
            numbers calculate Morans I for cells to the right (x) 
            and below (y) each cell, negative for cells to the 
            left (x) and above (y).
 -------------------------------------------------------------------------
 OUTPUTS: [m] single global value [m] for the analysed raster

 See also: localmoran(), moransI()

 -------------------------------------------------------------------
 NOTES: Based on code from the dissertation of Jo Wood:
        http://www.soi.city.ac.uk/~jwo/phd/04param.php

 Felix Hebeler, Geography Dept., University Zurich, March 2006.

0001 function m = lagmoran(raster, xlag, ylag);
0002 % PURPOSE: calculate Moran's I for a lag in x and/or y direction on
0003 %          a regular grid.
0004 % -------------------------------------------------------------------
0005 % USAGE: m = lagmoran(raster, xlag, ylag);
0006 % where: [raster] is the grid to calculated Moran's I for
0007 %        [xlag],[ylag] is the lag in x and y directions. Postive
0008 %            numbers calculate Morans I for cells to the right (x)
0009 %            and below (y) each cell, negative for cells to the
0010 %            left (x) and above (y).
0011 % -------------------------------------------------------------------------
0012 % OUTPUTS: [m] single global value [m] for the analysed raster
0013 %
0014 % See also: localmoran(), moransI()
0015 %
0016 % -------------------------------------------------------------------
0017 % NOTES: Based on code from the dissertation of Jo Wood:
0018 %        http://www.soi.city.ac.uk/~jwo/phd/04param.php
0019 %
0020 % Felix Hebeler, Geography Dept., University Zurich, March 2006.
0021 if nargin < 3
0022     error('This function needs 3 input arguments (inputraster, xlag, ylag)');
0023 end
0024 
0025 
0026 var=0;vari=0;varj=0;zi=0;zj=0;covar=0;
0027 ncols= size(raster,2);
0028 nrows= size(raster,1);
0029 
0030 % calculate mean of union of cells i U j
0031    % determine cells to include
0032 if (xlag>=0) istartcol =1;iendcol = ncols-xlag; jstartcol =1+xlag;jendcol = ncols;
0033 else istartcol = 1-xlag; iendcol = ncols; jstartcol =1;jendcol = ncols+xlag;
0034 end
0035 if (ylag>=0) istartrow =1;iendrow = nrows-ylag; jstartrow =1+ylag;jendrow = nrows;
0036 else istartrow = 1-ylag; iendrow = nrows; jstartrow =1;jendrow = nrows+ylag;
0037 end
0038 
0039 % average over union of cells
0040 zbar = mean(nanmean(raster(istartrow:iendrow,istartcol:iendcol))); % mean over i
0041 zbar = zbar + mean(nanmean(raster(jstartrow:jendrow,jstartcol:jendcol))); % mean over j
0042 
0043 for row=1:1:nrows;
0044       if ((row+ylag <= nrows) && (row+ylag > 0))
0045             for col=1:1:ncols;
0046                  if ((col+xlag <= ncols) && (col+xlag > 0))
0047                         zi = raster(row,col);
0048                         zj = raster(row + ylag,col + xlag);
0049                         if (isnan(zi) || isnan(zj)) break; end; %don't include NaNs
0050                         covar = covar + ((zi - zbar) * (zj - zbar));
0051                         vari = vari + ((zi - zbar) * (zi - zbar));
0052                         if ( ((row + ylag*2) >= nrows) || ((row + ylag*2) < 0) || ((col + xlag*2) >=ncols) || ((col + xlag*2) < 0))
0053                               varj = varj + ((zj - zbar) * (zj - zbar));
0054                         end
0055                  end
0056             var = (vari + varj); % Note that varj is in fact the var(i U j) - var (i)
0057             end
0058       end
0059 end              
0060 if (var ~= 0.0)
0061     m = covar/var;
0062 else
0063     m = 1.0 ;
0064 end