Calculating beta-diversity w raster data [renamed]

Question

Using ArcGIS 10, it is not possible to generate a raster attribute table for 32-bit float values. Instead, I have tried the RasterToNumPyArray method, which crashes since the file input raster is so large (I'm talking >7GB).

I've also tried converting the raster into an ASCII text file, and from there reading the file, line by line, into numpy arrays:

e.g.

import numpy as np
file = open(r"c:/yourfile.txt", 'r')
array_list = []
for line in file:
    array = np.loadfromstring(line.strip(), dtype = np.float32, sep = " ")
    array_list.append(array)

My main issue is that this is still pretty slow since the file has 71058 columns and 27820 rows.

Any advice? The troubling thing is that I can generate a list of unique values in ArcMap from symbolization, but there's no easy real way to move them into a table (besides copy and paste of each value).

Many thanks!

Answer 1

The comments indicate that this question seeks a method to compute "beta diversity" indexes for pairs of neighboring 5 x 5 blocks. This can be done by creating binary indicator grids (1=presence, 0=absence) for each of the 11 categories. The block maximum of such a grid will be a binary indicator of the presence of a category within each block. What we need is to XOR these block indicators within each pair of adjacent blocks: this is the contribution to the beta diversity index. The local sum of such "focal XORs" is the beta diversity.

Most raster GISes do not directly support a focal XOR, but one can be created from a focal sum: the focal sum within a 1 x 2 (or 2 x 1) neighborhood equals 0, 1, or 2; reclassifying the 2 as a 0 gives the focal XOR.

Here is a graphical representation of the workflow for a pair of vertically adjacent blocks:

The two blocks of numbers depict a 10 x 5 portion of the original grid. The parentheses show the two 5 x 5 blocks that comprise it. The numbers, which range from 1 through 11, are codes for the 11 classes.

Beneath these two blocks are color renderings of the values.

The next panel shows the 11 indicator grids. Colored squares have 1's and white squares have 0's. (The coloring scheme is the same as before.) To create an indicator grid for a class, equate the original grid to the class code, as in

[MyGrid] == 3

for the indicator grid of class 3. This is a fast operation, in part because the resulting grid is in integer format and will compress greatly with the native run length encoding: expect each one to occupy much less than 200 MB as a temporary dataset on disk. (The sizes might be only a megabyte or so if there is strong spatial correlation among the classes.)

The final panel shows the block maxima of these indicator grids. Note that the dimensions of this grid are each one-fifth those of the original grid, so there are just 1/25th as many cells. These block maxima should be quick to compute. (When using non-ESRI software, instead compute block sums using a fast Fourier transform and reclassify them by comparing them to 0. This will be extremely fast.)

The focal sums of these block maxima, using a 2 x 1 neighborhood, would be

2 2 2 2 2 2 1 2 2 1 0

Reclassifying 2->0 in each one gives

0 0 0 0 0 0 1 0 0 1 0

The local sum of these values equals 2: that's the beta diversity index for this pair of blocks. This approach produces an entire grid of such indexes, one value for every pair of vertically adjacent blocks (using a 2 x 1 neighborhood) or for every pair of horizontally adjacent blocks (using a 1 x 2 neighborhood).