This script calculates the distance from each cell to the nearest cell of a specific type (e.g. a road).

The script can be used used in the Carbon storage and wood production example, to calculate the distance to roads (node Distance_road) directly in gBay.

Download the script

Download the carbon sequestration network and spatial data to run the network with gBay.

Input:

• discrete raster (1 band for each state)
• number of the state of interest (e.g. roads = state number 0)

Output:

• hard evidence on a continuous node

  ### Calculates the distance to the nearest cell with a specific state (e.g. road)
  ### Input: discrete raster of categories (e.g. 1 = road, 0 = no road)
  ### Sets evidence on node continuous node (Distance_road)
  
  # import packages
  import os
  import sys
  import numpy
  from osgeo.gdalconst import *
  from osgeo import gdal
  import math
  from node_utils import *
  from scipy.spatial import distance
  
  # set names for labelling the script while running
  SCRIPT_NAME = "CONT_DISTANCE: "  # to add as a prefix for all messages
  CELLS_PRINTING = 10
  
  # SET FUNCTION PARAMETERS
  # name of input node
  input_name = "road"
  # number of the state in the input raster that defines the cells (e.g. roads), the distance to which we are interested in
  state_number = 0;
  # name of output node 
  output_name = "Distance_road"
  
  
  # function that finds cells of interest
  def isRoad(node, cell):
      return (getStateHighestLikelihood(node, cell) == state_number)
  
  # function that finds distance to the nearest cell of interest
  def findCloserRoadCell(cell, road_list, width):
  
      if (len(road_list) == 0):
          print("findCloserRoadCell: There are no roads in the map.")
          return None
  
      return min(distance.cdist([cell], road_list, 'euclidean')[0])
  
  # function that writes distance to the output node
  def process(dataset,nodes_data,iteration):
  
       # only use it at the beginning
      if (iteration != -1):
          return []
  
      node_road = getNodeByName(nodes_data, input_name)
      if (not node_road):
          print(SCRIPT_NAME, "ERROR Node", input_name, "is not in nodes_data")
          return []
  
      total_cells = dataset.RasterXSize * dataset.RasterYSize ;
  
      # get pixel size
      pixelsize = dataset.GetGeoTransform()[1]
  
      print(SCRIPT_NAME, "pixelsize:", pixelsize)
  
      road_cells = [];
  
      for cell in range(total_cells ):
      
          if (isRoad(node_road,cell)):
              road_cells.append(cell);
           
      print(SCRIPT_NAME, "There are", len(road_cells), "cells that are roads")
  
      distances=[]
  
      # Transform to 2-dimensional array
      road_coords = list(map(lambda x: [x / dataset.RasterXSize, x % dataset.RasterXSize], road_cells))
      
      # append distances to new array
      for cell in range(total_cells ):
          
          if (isNODATA(node_road, cell)):
              distances.append(PY_NODATA_VALUE)
          elif isRoad(node_road, cell):
              distances.append(0)
          else:
              distances.append(findCloserRoadCell([cell / dataset.RasterXSize, cell % dataset.RasterXSize], road_coords, dataset.RasterXSize) * pixelsize)
      
  
      new_nodes_data = [{'name' : output_name, 'type': PY_CONTINUOUS, 'data': distances}]
  
  
      for i in range(CELLS_PRINTING):
          print(SCRIPT_NAME, "cell:", i, ":" ,  new_nodes_data[0]['data'][i])
      
  
      if (validResultData(new_nodes_data, total_cells)):
          return new_nodes_data
  
      else:
          print ("Some error here: ")