==== Mean of neighbouring cells ====
This script calculates the mean of a continuous variable (e.g. % forest cover) in a specified window around each cell. It runs before the first iteration (i.e. to prepare the input data).

The script is used in the [[Habitat suitability]] example, to calculate the mean forest cover in a 500m-window around each cell. 

{{ :Mean_neighbour.zip | Download the script}}

{{ :habitat_netowrk.zip | Download the habitat suitability network}} and {{ :habitat_data.zip |spatial data}} to run this script in gBay. 

Input: 
  * continuous raster (1 band) 
  * window size (distance)

Output: 
  * hard evidence on a continuous node 

  ### Calculates the mean of a continuous variable (e.g. % forest cover)
  ### in a specified window around each cell. 
  ### Input raster: continuous raster (1 band) [Forest_cover]
  ### Sets hard evidence on continuous node [Forest_cover_neighb]
  # import packages
  import numpy
  import math
  from node_utils import *
  # set names for labelling the script while running
  SCRIPT_NAME = "CONT_NEIGHBOUR_MEAN_TEST: "  # to add as a prefix for all messages
  CELLS_PRINTING = 3
  # SET FUNCTION PARAMETERS
  # neighbourhood distance - size of the window [same units as raster]
  d = 500
  # name of input node
  input_name = "Forest_cover"
  # name of output node 
  output_name = "Forest_cover_neighb"
  # DEFINE THE FUNCTION
  def process(dataset,nodes_data,iteration):
    # only use it at the beginning
    if (iteration != -1):
        return []
    # get pixel size
    pixelsize = dataset.GetGeoTransform()[1]
    # number of pixels to take into account for the neighborhood
    npix = int(d/pixelsize)
    print(SCRIPT_NAME, "Taking into account:", npix, "pixels")
    input_data = getNodeByName(nodes_data, input_name)
    if (not input_data):
        print(SCRIPT_NAME, "ERROR Node", input_name, "is not in nodes_data")
        return []
    # create a list with only the prob of the state of interest (road)
    # create empty list
    data_list = []
    
    for cell in range(dataset.RasterXSize * dataset.RasterYSize ):
    
        # get value of state
        value = getValue(input_data,cell) 
        # append the value for this cell to the list 
        data_list.append(value)
    
    # check if the values are correct     
    print("Values:",(data_list[0:10]))
        
    # create 2D array from the vector (easier for edge problems)
    data_array = numpy.array(data_list)
    # replace no data values with nan
    data_array[data_array==PY_NODATA_VALUE]=numpy.nan
    wide = data_array.reshape(dataset.RasterYSize,dataset.RasterXSize)
    
    # create empty neighborhood array
    neighb = numpy.empty((dataset.RasterYSize,dataset.RasterXSize))
    
    # loop through cells and calculate mean of neighbouring cells
    for i in range(wide.shape[0]): 
        for j in range (wide.shape[1]):
            if i >= npix and j >= npix: 
                neighblist = wide[i-npix : i+npix+1, j-npix : j+npix+1] 
            elif  i >= npix:
                neighblist = wide[i-npix : i+npix+1,  : j+npix+1] 
            elif j >= npix:
                neighblist = wide[: i+npix+1, j-npix : j+npix+1] 
            else:
                neighblist = wide[: i+npix+1, : j+npix+1]         
            
            average =  neighblist[~numpy.isnan(neighblist)].mean()
    
            neighb[i,j] = average    
    # check maximum and minimum
    print("Maximum value of output is",neighb[~numpy.isnan(neighb)].max())
    print("Minimum value of output is",neighb[~numpy.isnan(neighb)].min())
    # write the values into the node_data format 
    # into a continuous node format
    
     # translate to 1D list
    neighb_list = neighb.flatten().tolist()
    print("Length of outputs is:", len(neighb_list))
    
    new_nodes_data = [{'name' : output_name, 'type': PY_CONTINUOUS, 'data': neighb_list}]
        
    for i in range(CELLS_PRINTING):
        print("cell:", i, ":" ,  new_nodes_data[0]['data'][i])
    return new_nodes_data