Neighborhood covariance for map comparison - Map Algebra help needed

[DarrenBoss]

I had posted in the Imagery and Raster forum a few days ago, but have not had any response as I think my question was poorly phrased and rather vague.  Also, I believe the Spatial Analyst forum may be a more appropriate place for my inquiry, so I apologize for the cross post!

I am looking at local map comparison metrics for maps with continuous attributes.  One of the metrics involves calculating the local covariance between two maps in a n x n moving window.  The formula for covariance is CV(x,y) = 1/N * SUM (x_i - x_bar)(y_i - y_bar). 

For example, in raster x, for the cell in row 1, column 1 of my moving window I need to subtract the local mean from this cell's value.  I do the same for the corresponding cell in raster y and multiply these values.  This calculation is performed for each cell in my moving window, the results for that neighborhood are summed, divided by n and assigned to the focal cell in the output raster.

I can envision a very slow, poor use of PixelBlocks where I define a PixelBlock of size n x n and access the results in the PixelBlock's safe array, perform the calculations and write to the output raster, but I feel this would be very inefficient.  I can also envision a method wherein I convert the source rasters to some form of text file and use clever indexing to get the pixel values I need, perform the calculation and write to the output.  If I use either of these approaches I will work in VBA.

I can sense that there is a way of doing this via map algebra, but I cannot quite come to a solution.  Perhaps someone (Bill Huber?!) has some ideas?

Edit - I have been reading through a lot of posts and my level of hope has increased somewhat.  In a post whuber stated, "For integer grids, though, focal stats with an appropriate custom neighborhood will pick out the neighboring value in any position you like."  Does anyone know how this is accomplished with integer grids?  Can the idea be extended to grids with real numbers?

Thanks for your help!

[DarrenBoss]

Thank you for your detailed response Bill!  I worked through another example with a 3x3 window and confirmed that this is exactly what I was hoping for.

Your assistance is greatly appreciated!!!

[DarrenBoss]

Hi Bill,

Thank you again for your assistance with the 'focal covariance' calculation. 

I would like to be able to use a variety of weights matrices when calculating the covariance, but I have not been able to determine a way of doing so.  The formula I now need to calculate is CV(x,y) = 1/N * SUM [w_i(xi - xbar)(yi - ybar)]  Where w_i is the weight for cell xi as specified in a kernel file/weights matrix.

I'll try to describe what I need to do with an example.

Example weights matrix/kernel file
[INDENT]3 3
0.5 1.0 0.5
1.0 1.0 1.0
0.5 1.0 0.5
[/INDENT]
Example grids x and y
[INDENT] is
1 2 0
1 0 2
3 0 0

is
1 3 5
1 3 5
1 3 5[/INDENT]

Determining weighted versions of xbar and ybar are easily accomplished with FocalMean and the custom weights matrix.  The problem is in calculating the properly weighted xybar.  One cannot simply calculate x*y and then determine the focal mean.  A properly weighted xybar would be generated as follows:

[INDENT]Apply weights to the window in
0.5 2.0 0.0
1.0 0.0 2.0
1.5 0.0 0.0[/INDENT]

[INDENT]Apply weights to the window in
0.5 3.0 2.5
1.0 3.0 5.0
0.5 3.0 2.5[/INDENT]

[INDENT]Multiply the weighted windows in and .
0.25 6.0 0.0
1.00 0.0 10.0
0.75 0.0 0.0[/INDENT]

My desired, weighted covariance is now simply the difference of FocalMean(weighted*weighted) - FocalMean(weighted)*FocalMean(weighted)

Can this be done with map algebra?

[DanPatterson_Retired]

you should be able to make a kernel file with your weights as inputs, see
http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//009z000000r7000000.htm
towards the bottom, they give an example of a kernel file which can be used with the focal mean.  This can be accessed through focal statistics in arctoolbox but I think you will have to preprocess before using map algebra

[DarrenBoss]

Thanks for your input Dan, but I am not sure I understand how it solves my problem.  I am aware of kernel files and provide an example of one in my previous post.  I am also familiar with focal statistics.  I think the preprocessing you mention is the 'black box' that I do not understand and the part I need a solution to.  The question is, how can I run a moving window over two rasters at the same time, apply the weights in the kernel file within each instance of the window, multiply the two windows and then obtain the focal mean of the resulting product?

I am hoping map algebra can be used in a clever way to achieve my goal!

[DarrenBoss]

So simple!

I would prefer that a user only inputs a path to a single weights file.  So I have to figure out how to read the file in Python, square all the values and write to a new file.  This should be fairly trivial.

Thanks for your help Bill!