python array processing: how to generate new values on pixel by pixel basis
Question: My problem is about array manipulation in python. Want I want to do is to
analyze a multiband raster with python and output another raster depending on
the results. Real case example: I have a raster with 3 bands and I want to see
if the different band values for the same pixel are in increasing or
decreasing order (from band 1 to band 3). If the values are in increasing
order, a new value, lets say "1" will be assigned to that particular pixel in
a new output raster (or can be also a new band in the same raster). If the
values are in decreasing order, the value "2" will be assigned. I have defined
the functions and read the array values but I m stuck at recursively build
lists for every pixels, and assign new values for a new raster.

I will post the code that I have until now:

    
    
    import arcpy
    import numpy as np
    
    arcpy.env.workspace = r"\\workspace"
    
    def decrease(A):
        return all(A[i] >= A[i+1] for i in range(len(A)-1)) # function to check if a list has decreasing values
    def increase(A):
        return all(A[i] <= A[i+1] for i in range(len(A)-1))  # function to check if a list has increasing values
    def no0(A,val):
        return[value for value in A if value !=val]   # function to check if a list has 0 values, and if so, strip them
    
    myArray = arcpy.RasterToNumpyArray(r"raster_3bands.img")  # my 3 bands raster
    
    print myArray.shape  # to see how many dimensions the array has
    for z in range(3): # where z is the number of bands, in my case 3
        for y in range(6): #the number of columns, i.e. the vertical axis, in my case 6
            for x in range(6): #the number of rows,i.e. the horizontal axis, in my case 6.
                a = myArray[z,y,x] # get all the values for the 3 bands
                print "The value for axes:  ", z,y,x, " is: ",a
    

What I am missing: 1\. the way to build for every pixel a list which will
store the three correspondent band values, so that I can later run the
Decrease and Increase functions on those lists 2\. a way to assign new values,
pixel by pixel and to store that array as a new raster.

Many thanks for your patience to read this,

Bogdan

So, here is the code for a 3 bands raster:

    
    
    import arcpy, os
    import numpy as np
    myArray = arcpy.RasterToNumPyArray(r"\\nas2\home\bpalade\Downloads\test.img")
    nbands = 3
    
    print "Array dimensions are: ", myArray.shape
    print "\n"
    print "The array is: "
    print "\n"
    print myArray
    
    
    increasing_pixels = np.product([(myArray[i] <= myArray[i+1]) for i in range(nbands-1)],axis=0)
    decreasing_pixels = np.product([(myArray[i] >= myArray[i+1]) for i in range(nbands-1)],axis=0)
    
    new_band = np.zeros_like(myArray[0])
    
    new_band[increasing_pixels] = 1
    new_band[decreasing_pixels] = 2
    
    print "\n"
    print "The new array is: "
    print "\n"
    print new_band
    

The result is attached as jpeg

I do not see the jpeg, so I copy/paste from my results window:

    
    
    Array dimensions are:  (3L, 4L, 4L)
    
    
    The array is: 
    
    
    [[[60 62 62 60]
      [64 64 63 60]
      [62 62 58 55]
      [59 57 54 50]]
    
     [[53 55 55 55]
      [57 57 56 55]
      [55 55 51 50]
      [52 50 47 45]]
    
     [[35 37 37 36]
      [39 39 38 36]
      [37 37 33 31]
      [34 32 29 26]]]
    
    
    The new array is: 
    
    
    [[1 1 1 1]
     [2 2 2 2]
     [0 0 0 0]
     [0 0 0 0]]
    >>> 
    


Answer: The first think to think about here is the structure of the code. When using
numpy (and especially in a case like this) it is rare that you need to apply
loops. You can simply apply vectorized operations along the axes.

In this case you can use array broadcasting to find which sets of pixels are
increasing and decreasing.

    
    
    increasing_pixels = (myArray[0] <= myArray[1]) * (myArray[1] <= myArray[2])
    

This will give you a boolean array where `True` pixels are those which
increase across the three bands. The same can be done with decreasing pixels
(this time `True` implies pixels that decrease across the band):

    
    
    decreasing_pixels = (myArray[0] >= myArray[1]) * (myArray[1] >= myArray[2])
    

These arrays can act as a mask for your raster, allowing you to identify the
pixels that you wish to change value. Below is an example of creating a new
band based on the increasing or decreasing values of pixels across the
original three bands.

    
    
    new_band = np.zeros_like(myArray[0])
    
    new_band[increasing_pixels==True] = 1
    new_band[decreasing_pixels==True] = 2
    

**EDIT:**

For an arbitrary number of bands, the above comparisons can be replaced by:

    
    
    increasing_pixels = np.product([(myArray[i] <= myArray[i+1]) for i in range(nbands-1)],axis=0)
    decreasing_pixels = np.product([(myArray[i] >= myArray[i+1]) for i in range(nbands-1)],axis=0)
    

