sparse_window_demo_2016_05_03.py

sparse_window_demo_2016_05_03.py

"""
Script:   sparse_window_demo.py
Author:   Dan.Patterson@carleton.ca
Modified: 2016_05_03
Purpose:
- Conversion of point data to raster
- To produce a moving window equivalent to find the cells that are 
  surrounded by a particular value
References:
- http://desktop.arcgis.com/en/arcmap/latest/tools/conversion-toolbox/
       point-to-raster.htm
- http://stackoverflow.com/questions/32108409/identifying-pairs-of-
       python-array-cells-separated-by-maximum-distance-in-a-large
Example:  slicing options
>>> a = np.arange(25).reshape((5,5))
>>> a
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19],
       [20, 21, 22, 23, 24]])
(1)....             (2)....             (3)....             (4)....
a                   a[0:,:3]            a[:2,:3]            a[:2,:2]
array([[0, 1, 2],   array([[0, 1, 2],   array([[0, 1, 2],   array([[0, 1],
       [3, 4, 5],          [3, 4, 5],          [3, 4, 5]])         [3, 4]])
       [6, 7, 8]])         [6, 7, 8]])  
Notes:
  distance - given in cell widths, numeric value not relevant
  A 3*3 window gives 1 unit of distance (d) in the cardinal directions and sqrt(2)*d 
  on the diagonals based on cell to cell centers.  A 3*3 window represents the extents
  of the cells.
"""
import numpy as np

np.set_printoptions(edgeitems=10,linewidth=80,precision=2,suppress=True,threshold=100)

def plot_sparse(arr):
    """plot a sparse matrix"""
    import matplotlib.pyplot as plt
    # binary,gray,spectral,prism,flag,rainbow,Paired,Blues,
    plt.imshow(arr, cmap="binary", interpolation='nearest')
    plt.show()    

def sparse_sub(a,val=1,rows=3,cols=3):
    """subsample a sparse grid array with a rows*cols window for cells
       that contain val
    """
    w0 = rows-2; w1 = cols-1   # window slice parameters
    r,c = np.where(a==val)     # determine the cells with val
    width, hght = a.shape               
    #pnts = list(zip(r,c))      # rows-columns to pnts...use below
    pnts = np.zeros((len(r),2),dtype=np.int32)  # (N,2) zeros array
    pnts[:,0] = r
    pnts[:,1] = c
    results=[]
    for p in pnts:
        r,c = p
        r0 = max(0,r-w0); c0 = max(0,c-w0)         # window bounds
        r1 = min(r+w1,width); c1 = min(c+w1,hght)
        #print("r0,r1,c0,c1 \n{},{},{},{}".format(r0,r1,c0,c1))
        result = np.where(a[r0:r1,c0:c1]==1)       # val in window
        if result:
            pairs = zip(result[0]+r0,result[1]+c0)     # zip row/column values
            results.append([p,len(pairs),list(pairs)]) # store results
            #print("p: {} ...pairs: {}".format(p,list(pairs)))
    return results

def grid_XYZ(a, sparse=True, sorted=True):
    """make a mesh grid from a 2D raster/grid data structure"""
    rows, cols = a.shape     
    XX, YY = np.meshgrid(np.arange(a.shape[1]), np.arange(a.shape[0]))
    dt = [("X",'int64'),("Y",'int64'),("Z",'int64')]
    #XYZ = np.vstack((XX.ravel(), YY.ravel(), a.ravel()), dtype=dt).T
    XYZ = np.array(list(zip(XX.ravel(),YY.ravel(),a.ravel())),dtype=dt)
    if sparse:
        #idx =  np.where(XYZ[:,2] != 0)  # pick all non-zero Z vlues
        idx =  np.where(XYZ['Z'] != 0) 
        XYZ = XYZ[idx]
    if sorted:
       XYZ = XYZ[np.argsort(XYZ ,order=('X','Y'))]
    return XYZ
    
def XYZ_grid(XYZ):
    """convert XYZ values to a sparse array"""
    x = XYZ['X']
    y = XYZ['Y']
    z = XYZ['Z']
    xyz = XYZ.tolist()
    x_min = x.min(); x_max = x.max()
    y_min = y.min(); y_max = y.max()
    a = np.zeros((y_max+1,x_max+1), dtype='int')
    for row in xyz:
        i,j,k = row
        a[j][i] = k
    return a
    
def data1():
    # Example 2-D habitat array (1 = data) from link
    a = np.array([[0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
                  [0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1],
                  [0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0],
                  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
                  [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
                  [1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1],
                  [1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1],
                  [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
                  [1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0],
                  [1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                  [1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]])
    return a
    
def data2():
    # Example 2-D  mine
    a = np.array([[0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
                  [0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1],
                  [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
                  [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                  [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
                  [0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
                  [0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0],
                  [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
                  [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
                  [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1]])
    return a

def data3():
    # Example 2-D  mine
    a = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
                  [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                  [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
                  [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]])
    return a
    
def data_xy():
    """make an array out of xy coordinates"""
    N = 50
    id = np.arange(N)
    x = np.random.randint(2,N) # generate N 0's or 1's
    y = np.random.randint(N)   # ditto
    xy = np.zeros((N,N), dtype='int32')
    return x, y, xy
    
def main(arr=None,to_graph=True):
    """run the demo with sample data and optional graphing"""
    if arr == None:
        arr = data1()
    results = sparse_sub(arr,val=1,rows=3,cols=3)  
    print("Input array...\n{}".format(arr))
    print("{:<12}{:<5}{:}".format("row/col","N","pnts"))
    for i in results:
        print("{:<12}{:<5}{}".format(*i))
    if to_graph:
        plot_sparse(arr)

         
#----------------------------------------------------------------------
# 
if __name__=="__main__":
    """See documentation above.
    """
    #arr = main()  # run on the test ndarray
    ## run a sample on a slice of the test data ie...
    #arr = main(data()[:3],True)  # ie, first 3 rows
    ## run on subsample on your own data
    #main(my_ndarray,False)
    ## another test
    #
    a = data3()
    a_f = grid_XYZ(a, sparse=False, sorted=False)  #  X,Y,Z,XYZ
    #print("\nFull array of x,y,z...\n{}\n".format(a_f))
    a_s = grid_XYZ(a, sparse=True, sorted=False)
    print("\nSparse array...\n{}".format(a_s))
    b = XYZ_grid(a_s)
    main(a)