Initial commit

README.md

+# cctrigrid
+
+A library for computing connected components on an unstructured triangular grid
+
+Authors: - Aiko Voigt (University of Vienna)
+         - Noam Rotberg (Otto von Guericke University of Magdeburg)
+         - Nicole Knopf (Karlsruhe Institute of Technology)
+         - Petra Schwer (Otto von Guericke University of Magdeburg)

cctrigrid/cctrigrid.py

+import numpy as np
+    
+import grid_functions as gridf
+import cubulation_functions as cubulf
+
+# make_cubulation
+#   Assign to each triangle a 3d-coordinate called cube coordinate.
+# INPUT: 
+#   start_triangle      int         the triangle index where the cubulatiion algorithm starts
+#   radius              int         the number of "outward" steps
+#   print_progress      Boolean     print progress of iteration (default False)
+# OUTPUT: 
+#   Returns a list of arrays of the following form:
+#   np.array([<triangle index>, <cube array>])
+#   where <cube array> is an array with three entries that encode the cube coordinate of the triangle with <triangle index>. 
+
+def make_cubulation(start_triangle, radius, print_progress=False): 
+    
+    # initialization
+    init = cubulf.EnhancedTri(start_triangle, np.array([0,0,0])) # assign coordinates (0,0,0) to start_triangle
+    cube_coordinates = [np.array([init.triangle, init.cube], dtype=object)]
+    visited_triangles = [init.triangle]
+    outmost = [init] # list for EnhancedTri that were added last round and thus lie "at the boundary of discovery"
+    
+    # colour the edges of start_triangle with 0,1,2
+    edge_colours = { gridf.get_edges_of_cell(start_triangle)[x] : x for x in (0,1,2) }
+    
+    # expand outwards 
+    for n in range(radius):
+        # update outmost triangles
+        outmost = cubulf.cubing_next_round(cube_coordinates, visited_triangles, outmost, edge_colours)
+        # optional: show progress
+        if print_progress: print('Round ', n+1, 'finished. Visited ', len(cube_coordinates), \
+                                 'triangles, thereof ', len(outmost), 'new.')
+    
+    # return coordinates after shifting them to positivity
+    return cubulf.shift_coordinates(cube_coordinates, radius)
+    
+    
+# make_connected_2d_components:
+#   From given cubuluation, compute connected components and output them as a list of lists
+# INPUT:
+#   cubulation          list of arrays      output of make_cubulation 
+#   field_cube          numpy array         input field in cubulated 3d-coordinates
+#   connectivity        string              use 'vertex' (default) or 'edge' connectivity
+# OUTPUT: 
+#   - Returns a list of lists. 
+#   - Each inner list corresponds to a connected component and contains 
+#     all triangle indices that belong to this connected component
+
+def make_connected_components_2d(cubulation, field_cube, connectivity = 'vertex'): 
+    
+    # compute connected components on field_cube
+    # components_cube is a 3-dimensional array that contains
+    # the connected component index of a triangle with cube coordinates (x,y,z) at position component_array[x][y][z]
+    components_cube = cubulf.make_connected_components(field_cube, connectivity)
+    
+    # a component list of triangle indices is generated for each connected component:
+    # returns a list containing all component lists
+    components = cubulf.make_list_of_components (cubulation, components_cube)
+   
+    return components
+
+
+# make_connected_3d_components:
+#   From given cubuluation, compute connected components and output them as a list of lists
+# INPUT:
+#   cubulation          list of arrays      output of make_cubulation 
+#   field_cube          numpy array         input field in dimensions of (lev x cubulated coordinates)
+#   connectivity        string              use 'vertex' (default) or 'edge' connectivity
+# OUTPUT: 
+#   - Returns a list of lists. 
+#   - Each inner list corresponds to a connected component and contains all tuples of (level, triangle index) that
+#     belong to this connected component
+
+def make_connected_components_3d(cubulation, field_cube, connectivity = 'vertex'):
+    
+    nlev = np.shape(field_cube)[0] # number of vertical levels
+    
+    # Step 1: compute connected components for each level separately, 
+    # in the following these are refered to as 2d connected components
+    components2d = []
+    for lev in range(nlev):
+        aux = make_connected_components_2d(cubulation, field_cube[lev], connectivity)
+        for items in aux:
+            components2d.append([lev]+sorted(items))
+    del aux
+    
+    # components is now a list of sublists, each sublist has the level as its first
+    # entry followed by the indices of connected triangles on that level 
+    # --> e.g., components[20] = [lev, 1058, 1059, 2099, 2100, ... ]
+    
+    # Step 2: we merge the 2d connected components in the vertical --> this will yield the 
+    # connected components in the 3-dimensional space ("the 3d connected components")
+    # to this end, we regard each 2d connected component as a node of an undirected graph and 
+    # use the connected component finding for graphs
+    
+    # Step 2a: we look for connections between nodes, i.e., cases when two 2d connected compoents
+    # are in adjacent levels and share at least on triangle with each other
+    # this creates a list of pairs that uses the indices components2d and are the edges of the graph
+    # ---> e.g., a pair might lool like pairs[0] = [0,1], in which case if would consist
+    #      of the first and second inner list of components2d
+    pairs =set() # we make use of sets as it automatically removes duplicate pairs
+    for item1 in components2d:
+        for item2 in components2d:
+            if np.abs(item1[0]-item2[0])==1: # the two 2d connected components need to be adjacent in the vertical
+            # is there overlap between the two 2d connected components?
+                if not set(item1[1:]).isdisjoint(item2[1:]): # --> enter if there is overlap
+                    pair = tuple(sorted([components2d.index(item1), components2d.index(item2)]))
+                    pairs.add(pair)
+    # turn pairs into a sorted list 
+    pairs = sorted(list(pairs))   
+
+    # we use the implmentation of the networkx library; this choice is 
+    # inspired by https://stackoverflow.com/a/62545221/10376696
+    import networkx as nx
+    # create ubdirected graph based on its edges, which are defined by pairs
+    graph = nx.Graph(pairs)
+   
+    # Step 2b: we also add all 2d connected components as nodes, this is needed
+    # to catch cases when a 2d connected component is not connected to another 2d connected
+    # components in a level above or below, i.e., it does not appear in "pairs"
+    graph.add_nodes_from(range(0, len(components2d)))
+                        
+    # Step 2c: we now use connected component search on the graph to identify 2d connected components 
+    # that belong together
+    # --> E.g., assume the graph consists of graph.nodes=[0, 1, 100] and 
+    #     graph.edges=[(0,1), (0,2), (1,3), (3,4), (7,8), (7,10), (10,11), (15,16), (20,24)],
+    #     then this should result in [{100}, {0,1,2,3,4}, {7,8,10,11}, {15,16}, {20,24}].
+    #     That is, there are five 3d connected components.
+    # "result" is a list of sets, which each set defining a 3d connected component
+    # for each set, the entries refer to the indices of components2d
+    result = list(nx.connected_components(graph))
+    del graph
+    
+    # Step 3 (final step): we create the list of 3d connected components, with each
+    # 3d connected component being given as a sublist of (lev, triangle) tuples
+    components3d=list()
+    # loop over the sets of result
+    for cc3d in result:
+        cc3d_tuples = [] # a list of (lev, triangle) tuples that belong to this 3d connected component
+        # loop over the 2d connected components, or more precisely their indices in the list components
+        for cc2d in cc3d: 
+            # loop over the triangles that belong to a 2d connected component,
+            # remember that the first entry is the level (see step 1)
+            lev = (components2d[cc2d])[0]
+            tri = (components2d[cc2d])[1:]
+            for ind_tri in tri: 
+                # generate the tuples and append them to the list cc3d_tuples that define the 3d connected component
+                aux = tuple((lev,ind_tri))
+                cc3d_tuples.append(aux)
+        components3d.append(cc3d_tuples)
+    
+    return components3d
\ No newline at end of file

cctrigrid/config.py

+ds_grid = None
\ No newline at end of file

cctrigrid/cubulation_functions.py

+import numpy as np
+import cc3d # external library for connected component labeling on cubic grids
+
+import grid_functions as gridf
+
+
+# class to enhance a triangle with its cube coordinate
+#   for a given triangle create an object consisting of the triangle index and the triangle's cube coordinate
+class EnhancedTri:
+    def __init__(self, triangle, cube):
+        self.triangle = triangle
+        self.cube = cube
+    def print(self):
+        print("Triangle index:", self.triangle, " Cube coordinate:", self.cube)
+        
+
+# colouring edges of a triangle
+#   situation:  the colours of two edges are already known, 
+#               the third edge must be coloured
+#   called from colour_new_edges in this special situation
+# INPUT:
+#   first_known_edge      int     index of a already coloured edge
+#   second_known_edge     int     index of a already coloured edge
+#   new_edge              int     index of the not yet coloured edge
+#   edge_colours          dict    dictionary storing to each edge index the colour
+# OUTPUT:
+#   None, but updates edge_colours
+def colour_exactly_one_new_edge(first_known_edge, second_known_edge, new_edge, edge_colours):
+    
+    # get the colours of the already coloured edges
+    first_colour = edge_colours[first_known_edge]
+    second_colour = edge_colours[second_known_edge]
+    
+    # determine new colour, i.e. the one not taken from (0,1,2) by first_colour and second_colour
+    for i in (0,1,2): 
+        if not (first_colour == i or second_colour == i):
+            edge_colours[new_edge] = i  # colour found -> update edge_colours
+            return
+        
+        
+# colouring edges of a triangle
+#   from old_triangle with already coloured edges determine the edge_colours of the adjacent new_triangle
+# INPUT:
+#   old_triangle_edges      set     contains the edge indices of old_triangle's edges
+#   new_triangle_edges      set     contains the edge indices of new_triangle's edges
+#   joint_edge              int     index of the shared edge from old_triangle and new_triangle
+#   edge_colours            dict    dictionary storing to each edge index the colour
+# OUTPUT:
+#   None, but updates edge_colours
+def colour_new_edges(old_triangle_edges, new_triangle_edges, joint_edge, edge_colours):
+    
+    # determine the edges of new_triangle that are not shared with old_triangle
+    new_edges   = new_triangle_edges.difference( old_triangle_edges ) 
+    first_edge  = new_edges.pop()
+    second_edge = new_edges.pop()
+    
+    # are the found edges uncoloured?
+    first_is_coloured = first_edge in edge_colours
+    second_is_coloured = second_edge in edge_colours
+    
+    # case: first_edge is coloured
+    if first_is_coloured: 
+        if second_is_coloured: 
+            return # both known -> nothing to colour
+        else: 
+            # colour second_edge using the colours of first_edge and joint_edge
+            colour_exactly_one_new_edge( first_edge, joint_edge, second_edge, edge_colours )
+    
+    # case: second_edge is coloured
+    if second_is_coloured: 
+        # colour first_edge using the colours of second_edge and joint_edge
+        colour_exactly_one_new_edge( second_edge, joint_edge, first_edge, edge_colours )
+    
+    # case: both edges are uncoloured
+    else: 
+        # get old_triangle's edges that are not joined with new_triangle
+        comparison_edges = old_triangle_edges.difference( new_triangle_edges )
+        first_comparison_edge = comparison_edges.pop()
+        second_comparison_edge = comparison_edges.pop()
+        
+        # get vertices of first_edge and first_comparison_edge
+        first_vertices = gridf.get_vertices_of_edge( first_edge )
+        first_comparison_vertices = gridf.get_vertices_of_edge( first_comparison_edge )
+        
+        # if first_edge and first_comparison_edge have no shared vertices:
+        #   the two edges must be parallel and thus, have the same colour
+        #   else, first_edge must be parallel to second_comparison_edge and thus, have the same colour
+        parallel = set(first_vertices).isdisjoint( set(first_comparison_vertices) )
+        
+        if parallel: # parallel: same colour
+            edge_colours[ first_edge ] = edge_colours[ first_comparison_edge ]
+            edge_colours[ second_edge ] = edge_colours[ second_comparison_edge ] 
+        
+        else: # not parallel: first_edge and second_comparison_edge are parallel and same-coloured
+            edge_colours[ first_edge ] = edge_colours[ second_comparison_edge ]
+            edge_colours[ second_edge ] = edge_colours[ first_comparison_edge ]
+            
+        
+# compute cube coordinate of a new triangle
+#   given the old_triangle's cube_coordinate derive the cube coordinate of the adjacent new_triangle 
+#   by using the colour of their joint_edge 
+# INPUT:
+#   old                     EnhancedTri     the old EnhancedTri with known cube_coordinates
+#   joint_edge              int             index of the shared edge from old_triangle and new_triangle
+#   edge_colours            dict            dictionary storing to each edge index the colour
+# OUTPUT:
+#   new_cube_coordinate, the cube_coordinate of new_triangle
+def determine_cube_coordinates (old, old_triangle_edges, new_triangle_edges, joint_edge, edge_colours):
+    
+    # get colour of the joint_edge and call it direction
+    direction = edge_colours[joint_edge]
+    # the direction encodes the parallel class of the joint edge
+    # this encodes also the direction in which one has to walk from old_triangle to new_triangle
+    # this gives the coordinate in which the cube_coordinates of old_- and new_triangle differ
+    
+    # determine direction-vector (1,0,0), (0,1,0) or (0,0,1) by which the cube_coordinates will differ
+    direction_vector = np.array([0,0,0])
+    direction_vector[direction] = 1 
+    
+    # use invariant: sum of coordinates musst be 0 or 1
+    # determine wether the direction vector has to be added or subtracted
+    # by using the invariant that the sum of all valid cube_coordinates has to be 0 or 1
+    old_coordinate_sum = sum(old.cube[i] for i in (0,1,2))
+    if old_coordinate_sum == 0:
+        # old triangle's coordinate sum is 0 thus new ones has to be 1 -> add
+        new_cube_coordinate = old.cube + direction_vector
+    else: 
+        # old triangle's coordinate sum is 1 thus new ones has to be 0 -> subtract
+        new_cube_coordinate = old.cube - direction_vector
+    
+    return new_cube_coordinate
+
+
+# iterate over the outmost triangles and compute the cube_coordinates of their adjacent triangles
+# INPUT:
+#   cube_coordinates        array       EnhancedTri's of the triangles whose cube_coordinates are already computed
+#   visited_triangles       list        indices of the triangles whose cube_coordinates are already computed
+#   outmost                 list        EnhancedTri's of the triangles considered in the round before   
+#   edge_colours            dict        dictionary storing to each edge index the colour
+# OUTPUT:
+#   updated outmost, EnhancedTri's of the new triangles considered this round
+def cubing_next_round(cube_coordinates, visited_triangles, outmost, edge_colours):
+    
+    # list for triangles that will be outmost in next iteration
+    new_outmost = []
+    
+    for old in outmost: # consider all EnhancedTri's at the border of visited
+        for neigh in gridf.get_neighbors_of_cell(old.triangle ): # consider all neighbours
+            new = EnhancedTri(neigh, 0)
+            if new.triangle == -9999: # use feature: triangles at the grid border claim to be adjacent to -9999
+                break
+            else:
+                if new.triangle not in visited_triangles: # EnhancedTri new has no cube_coordinate yet
+
+                    # add new to new_outmost for next round hereafter and update visited_triangles
+                    new_outmost.append( new )
+                    visited_triangles.append ( new.triangle )
+
+                    # preparation: obtain edges of new_triangle and old_triangle and derive their joint_edge
+                    old_triangle_edges = set(gridf.get_edges_of_cell(old.triangle)) # this is a set for python reasons...
+                    new_triangle_edges = set(gridf.get_edges_of_cell(new.triangle))
+                    joint_edge = ( new_triangle_edges & old_triangle_edges ).pop()
+
+                    # colour the edges of new.triangle
+                    colour_new_edges (old_triangle_edges, new_triangle_edges, joint_edge, edge_colours)
+
+                    # get cube coordinates for EnhancedTri new and update the array of cube_coordinates
+                    new.cube = determine_cube_coordinates (old, old_triangle_edges, new_triangle_edges, joint_edge, edge_colours)
+                    cube_coordinates.append(np.array([new.triangle, new.cube], dtype=object))
+    
+    return new_outmost
+
+
+# shift cube coordinates such that they are all positive
+# INPUT:
+#   cube_coordinates    array   stores triangle - cube_coordinate pairs
+#   radius              int     same as in make_cube_coordiantes
+# OUTPUT:
+#   shifted cube_coordinates
+def shift_coordinates (cube_coordinates, radius):
+    shift = int(radius/2)
+    
+    # update all coordinates
+    for entry in cube_coordinates:
+        cube = entry[1] # get cube_coordinate 
+        
+        # consider all three directions
+        for direction in (0,1,2):
+            cube[direction] += shift # shift coordinate
+    
+    return cube_coordinates
+
+
+# compute connected components
+# INPUT:
+#   field_cube     array    obtained as part of field prepaation
+#   connectivity   string   use 'vertex' (default) or 'edge' connectivity
+# OUTPUT:
+#   component_cube, 3d array containing the component indices
+def make_connected_components(field_cube, connectivity = 'vertex'):
+    
+    import sys
+    
+    # translate connectivity - string into input value for 3d connected component labeling
+    if connectivity == 'edge':
+        connectivity_value = 6
+    else:
+        connectivity_value = 26     
+        if connectivity != 'vertex':    # use default 'vertex' connectivity for invalid input
+            sys.exit("Invalid input: only 'vertex' or 'edge' connectivity are allowed.")
+            
+    # calling external 3d connected component labeling
+    component_cube = cc3d.connected_components(field_cube, connectivity = connectivity_value)
+    
+    return component_cube
+
+
+# which triangle belongs to which component?
+# INPUT:
+#   cubulation         array       output of make_cube_coordinates   
+#   component_cube     3d arry     output of make_connected_components
+# OUTPUT:
+#   Returns a list of lists. 
+#   Each inner list corresponds to one connected component and contains the indices
+#   of all triangles that belong to it
+def make_list_of_components (cubulation, component_cube):
+    
+    # number of found components on cubical grid
+    ncomp = component_cube.max()
+    
+    # initialize list of ncomp lists, with each list holding the triangles that
+    # belong to a specific connected component
+    component_list = []
+    for n in range(ncomp):
+        component_list.append([])
+    
+    # determine component
+    for entry in cubulation:
+        triangle  = entry[0]        # index of triangle
+        cube      = entry[1]        # cube coordinate
+        component = component_cube[cube[0]][cube[1]][cube[2]]
+        if component != 0:  # triangle belongs to a connected component
+            # add triangle to its component list
+            component_list[component-1].append(triangle)
+    
+    return component_list
\ No newline at end of file

cctrigrid/grid_functions.py

+# need to get access to the grid dataset
+import config
+ds_grid = config.ds_grid
+
+# get_neighbors_of_cell:
+#   returns edge neighbors of a given triangle
+def get_neighbors_of_cell(triangle):
+    return ds_grid.neighbor_cell_index[:,triangle].values
+   
+# get_edges_of_cell:
+#   returns the three edges of a given triangle
+def get_edges_of_cell(triangle):
+    return ds_grid.edge_of_cell[:, triangle].values
+    return edges.values
+
+
+# get_vertices_of_edge:
+#   returns the two vertices of a given edge
+def get_vertices_of_edge(edge):
+    return ds_grid.edge_vertices[:, edge].values
\ No newline at end of file

cctrigrid/prepare.py

+import numpy as np
+import xarray as xr
+
+# path_input:
+#   create globale variables to open the xArray dataset from gridfile and datafile
+# INPUT:
+#   path        string      location where gridfile and datafile are stored
+#   gridfile    string      name of the file containing the grid
+#   datafile    string      name of the file containing the field
+
+# we only need the following information about the horizontal triangular grid:
+# - the edge neighbors of a given triangle
+# - the three edges of a given triangle
+# - the two vertices of a given edge 
+
+
+# implementation for ICON model
+def grid_icon(path, file):
+    
+    ds_grid = xr.open_dataset(path + file)
+    ds_grid = ds_grid[['neighbor_cell_index', 'edge_of_cell', 'edge_vertices']]
+    
+    # in the ICON model grid, the indexing of the triangle cells and vertices 
+    # starts with 1 and not with 0 as assumed by cctrigrid --> we need to subtract 1 here
+    ds_grid['neighbor_cell_index'] = ds_grid['neighbor_cell_index'] - 1
+    ds_grid['edge_of_cell'] = ds_grid['edge_of_cell'] - 1
+    ds_grid['edge_vertices'] = ds_grid['edge_vertices'] - 1
+    
+    # after the substraction of 1, triangles at the grid border claim to be adjacent to -1
+    # --> we set these values to -9999 to clearly flag that there is no neighboing triangle
+    #     note: the -9999 value as a flag is used in the function "cubing_next_round"
+    ds_grid['neighbor_cell_index'] = xr.where(ds_grid['neighbor_cell_index']!=-1, ds_grid['neighbor_cell_index'], -9999)
+    
+    return ds_grid
+
+# prepare_field
+#   returns the data field as a numpy array, also applies threshold-based mask
+# INPUT:
+#   file: file containing the data field
+#   var:  variabe name
+#   threshold: set field to one if field>=threshold, and to zero otherwise
+# OUTPUT:
+#   field both on the triangular grid as well as on the cubulated grid
+
+# implementation for ICON model
+
+# single-level data
+def field_icon(path, file, var, threshold, cubulation):
+    
+    # field on triangular grid
+    field = xr.open_dataset(path + file)[var].squeeze().values
+    field[field>=threshold] = 1
+    field[field< threshold] = 0
+    
+    # determine size of array that holds cubulation
+    array_size=0
+    for ind in [0,1,2]:
+        xlist=list()
+        for entry in cubulation:
+            xlist.append(entry[1][ind])
+        array_size=max(array_size, max(xlist))
+    array_size=array_size+1 # need to add one as cubulation indices start with 1 instead of 0
+    
+    # map field from triangular grid to field on cubic grid
+    field_cube = np.zeros((array_size, array_size, array_size), dtype = 'int')
+    for entry in cubulation:
+        triangle = entry[0]     # index of triangle
+        cube = entry[1]         # cube coordinate
+        field_cube[cube[0], cube[1], cube[2]] = field[triangle]
+    
+    return field, field_cube
+
+# data on multiple levels
+def field_icon_lev(path, file, var, threshold, cubulation):
+    
+    # field on triangular grid
+    field = xr.open_dataset(path + file)[var].squeeze().values
+    field[field>=threshold] = 1
+    field[field< threshold] = 0
+    
+    nlev = np.shape(field)[0]
+    
+    # determine size of array that holds cubulation
+    array_size=0
+    for ind in [0,1,2]:
+        xlist=list()
+        for entry in cubulation:
+            xlist.append(entry[1][ind])
+        array_size=max(array_size, max(xlist))
+    array_size=array_size+1 # need to add one as cubulation indices start with 1 instead of 0
+    
+    # map field from triangular grid to field on nlev x cubic grid
+    field_cube = np.zeros((nlev, array_size, array_size, array_size), dtype = 'int')
+    for lev in range(nlev):
+        for entry in cubulation:
+            triangle = entry[0]     # index of triangle
+            cube = entry[1]         # cube coordinate
+            field_cube[lev, cube[0], cube[1], cube[2]] = field[lev, triangle]
+    
+    return field, field_cube
\ No newline at end of file