Differences

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

--- calculate_distance [2019/02/26 11:53] – created stritiha
+++ calculate_distance [2023/04/21 15:48] (current) – irladmin
@@ Line 1: / Line 1: @@
 ==== Calculate distance to specific cell type ====
+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.
+{{ :Distance.zip | Download the script}}
+{{ :carbon_network.zip |Download the carbon sequestration network}} and {{ :carbon_data.zip |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: ")