gBay Wiki

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.

Download the script

Download the habitat suitability network and 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