morans_I

 PURPOSE: calculate global Moran's I for an input grid (matrix) by calculating all 
          local Moran's I for a given moving windows size using a weight matrix. 
 -------------------------------------------------------------------
 USAGE: M = moransI(grid, W, s);
 where: [grid] is the matrix to analyse
        [W] is the normalized weight matrix of the size the local Moran's
            I will be calculated for (uneven sized!)
        [s] is an optional flag to use zscores of input values for
        calculation. Set to 'true' if zscores of local grid should be
        calculated. Leave blank if not desired or input values are already
        standardized. 
 -------------------------------------------------------------------------
 OUTPUTS:
        [M] matrix of all local Moran's I 
 -------------------------------------------------------------------
 NOTES: Weight matrix needs to be 'moving window' style, not contiguity
        matrix: Moran's I is calculated and weighted for neighbours to center cell.
        Matrix needs to be normalized (weights sum to 1) and center cell weight 
        will be set to 0 if not already. Uses localmoran.m
        -> Use nanmean(M(:)) to get the average global Moran's I.

 See Anselin (1995, 'LISA.', Geogr. Analysis 27(2),p.93f) for details on 
 standardized variables in calculation of local Moran's I. 

 EXAMPLE:  M = moransI(rand(20,20),ones(5,5),'true')

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

0001 function M = morans_I(grid,W,s)
0002 % PURPOSE: calculate global Moran's I for an input grid (matrix) by calculating all
0003 %          local Moran's I for a given moving windows size using a weight matrix.
0004 % -------------------------------------------------------------------
0005 % USAGE: M = moransI(grid, W, s);
0006 % where: [grid] is the matrix to analyse
0007 %        [W] is the normalized weight matrix of the size the local Moran's
0008 %            I will be calculated for (uneven sized!)
0009 %        [s] is an optional flag to use zscores of input values for
0010 %        calculation. Set to 'true' if zscores of local grid should be
0011 %        calculated. Leave blank if not desired or input values are already
0012 %        standardized.
0013 % -------------------------------------------------------------------------
0014 % OUTPUTS:
0015 %        [M] matrix of all local Moran's I
0016 % -------------------------------------------------------------------
0017 % NOTES: Weight matrix needs to be 'moving window' style, not contiguity
0018 %        matrix: Moran's I is calculated and weighted for neighbours to center cell.
0019 %        Matrix needs to be normalized (weights sum to 1) and center cell weight
0020 %        will be set to 0 if not already. Uses localmoran.m
0021 %        -> Use nanmean(M(:)) to get the average global Moran's I.
0022 %
0023 % See Anselin (1995, 'LISA.', Geogr. Analysis 27(2),p.93f) for details on
0024 % standardized variables in calculation of local Moran's I.
0025 %
0026 % EXAMPLE:  M = moransI(rand(20,20),ones(5,5),'true')
0027 %
0028 % Felix Hebeler, Geography Dept.,de University Zurich, March 2006.
0029 
0030 %% Check if standardising should be done
0031 if exist('s','var')
0032     if strcmp(s,'true');
0033         grid=zscore(grid);
0034     elseif strcmp(s,'false')
0035         %do nothing
0036     else
0037         error('Invalid option for s: set [true] to calculated zscores to determine local Moran or leave blank if values are already standardized.');
0038     end
0039 end
0040 if (mod(size(W,1),2)| mod(size(W,2),2))~=1
0041    error('Weight matrix W needs to have uneven size (eg. 5x5)') 
0042 end
0043 %% Do local Morans I calc of the grid.
0044 M = NaN(size(grid,1),size(grid,2));
0045 wsx=floor(size(W,1)/2);
0046 wsy=floor(size(W,2)/2);
0047 % Do local morans I calc for moving window ws
0048 for row=1+wsy:1:size(grid,1)-wsy;
0049     for col=1+wsx:1:size(grid,2)-wsx;
0050         M(row,col) = get_moran(grid(row-wsx:row+wsx,col-wsy:col+wsy),W);
0051     end
0052 end
0053 
0054 %% calculate local Moran's I
0055 function m=get_moran(raster,W)
0056 ncols= size(raster,2);
0057 nrows= size(raster,1);
0058 zi = raster(ceil(nrows/2),ceil(ncols/2));%  value of center cell (note: no weight applied!)
0059 if (isnan(zi));
0060     m=NaN; 
0061     return; 
0062 end;
0063 raster=raster.* W; % Weight values in window
0064 raster(ceil(nrows/2),ceil(ncols/2))=0; %set center cell to zero to exclude zi from sum
0065 zj = nansum(raster(:)); % sum of weighted values excluding zi
0066 m = zi * zj; % calculate local Moran's I and return