Differences

This shows you the differences between two versions of the page.

--- calculate_mean_value_in_neighbouring_cells [2019/02/25 17:33] – stritiha
+++ calculate_mean_value_in_neighbouring_cells [2023/04/21 16:04] (current) – irladmin
@@ Line 1: / Line 1: @@
-==== Mean of neighboring cells ====
+==== 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).
@@ Line 9: / Line 9: @@
 Input:
-  * continuous raster (1 band) [Forest_cover]
+  * continuous raster (1 band)
-  * window size (distance) [500 m]
+  * window size (distance)
 Output:
-  * hard evidence on a continuous node [Forest_cover_neighb]
+  * 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