some updates

Author: rvalavi

Cargo.lock

@@ -1,6 +1,6 @@
 [package]
 name = "rust_conn"
-version = "0.2.2"
+version = "0.2.3"
 edition = "2021"
 
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

Cargo.toml

@@ -1,3 +1,8 @@
+# Version 0.2.3
+* Fixed the isolated pixel calculations by duplicating its values to a fake neighbour.
+* Increased the capacity of graph to take more neighbours for high resoltion rasters.
+* Changed the default of `sigma` paramter to 1.
+
 # Version 0.2.2
 * Removed the dependency on GDAL, replacing it with `rasterio`.
 * Updated logic to ignore overview level 1 (both when reading and generating) since it corresponds to the original resolution and requires no additional computation.
@@ -7,3 +12,11 @@
 * Added Euclidean and Haversine distances for projected and unprojected grids, respectively, instead of using cell units.
 * Some changes in the arguments names including `max_cost`, `window_size` and `outer_window`.
 * Some small fixes for graph calculation.
+
+# Version 0.2.0
+Added geographic distance (in km) instead of cell units
+Replaced the scale argument with max_cost
+Various minor improvements
+
+# Version 0.1.0
+* Initial version.

NEWS.md

@@ -13,7 +13,7 @@ def connectedness(
         window_size=3, 
         outer_window=9,
         levels=[2, 4, 8, 16],
-        sigma=None,
+        sigma=1,
         option=3,
         n_threads=None,
         filename=""
@@ -62,7 +62,7 @@ def connectedness(
         Default is [2, 4, 8, 16].
     sigma : float, optional
         Standard deviation of the Gaussian kernel used for smoothing. 
-        Default is None (smoothing disabled).
+        Default is 1. Zero or None for disabling smoothing.
     option : int, optional
         Option flag to generate the connected condition from connectedness and input condition.
         Default is 3:
@@ -92,7 +92,7 @@ def connectedness(
         print(f"Notice: 'outer_window' was smaller than 'window_size' and has been adjusted to {window_size}.")
         outer_window = window_size
 
-    # Read raster overviews as a dictionary
+    # Read raster overviews as a dictionary; this checks levels as well
     data_dict, tran_dict, is_geo = read_raster(file=condition_file, gdf=polygon_mask, levels=levels, expand_px=outer_window)
 
     conn_array = connectivity(
@@ -108,7 +108,7 @@ def connectedness(
     )
 
     # Smooth the output array with Gaussian filtering
-    if sigma is not None:
+    if sigma is not None and sigma != 0:
         sigma = max(sigma, 1)
         conn_array = smoothing_filter(conn_array, sigma=sigma)
 
@@ -132,7 +132,7 @@ def beri(
         window_size=3, 
         outer_window=9,
         levels=[2, 4, 8, 16],
-        sigma=None,
+        sigma=1,
         n_threads=None,
         filename=""
     ):
@@ -187,7 +187,7 @@ def beri(
         Default is [2, 4, 8, 16].
     sigma : float, optional
         Standard deviation for the Gaussian kernel if smoothing is applied to input layers.
-        Default is None (no smoothing).
+        Default is 1. Zero or None for disabling smoothing.
     n_threads : int, optional
         The number of CPU cores for parallel processing. 
         Default is None (all available cores).
@@ -212,7 +212,7 @@ def beri(
         print(f"Notice: 'outer_window' was smaller than 'window_size' and has been adjusted to {window_size}.")
         outer_window = window_size
 
-    # Read raster overview as a dictionary
+    # Read raster overview as a dictionary; this checks levels as well.
     data_dict, tran_dict, is_geo  = read_raster(file=condition_file, gdf=polygon_mask, levels=levels, expand_px=outer_window)
 
     # Insert current climate as the first element in the list (this is important) before reading
@@ -240,7 +240,7 @@ def beri(
     )
 
     # Smooth the output array with Gaussian filtering
-    if sigma is not None:
+    if sigma is not None and sigma != 0:
         out_array = smoothing_filter(out_array, sigma=sigma)
 
     if len(filename) > 3:

connectivity/core.py

@@ -63,7 +63,6 @@ def create_overviews(input_raster, output_raster=None, overview_levels=[2, 4, 8,
         print(f"Error creating overviews: {e}")
 
 
-
 def overview_info(file_path):
     """Display information about all available overviews and their actual dimensions.
      
@@ -94,7 +93,6 @@ def overview_info(file_path):
                 print(f"  Level {level}: {w} x {h}")
 
 
-
 def read_raster(file, gdf=None, levels=None, expand_px=3):
     """Reads specified overview levels from a multi-band Cloud-Optimized GeoTIFF (COG) file
     and stores them in a dictionary.
@@ -223,7 +221,6 @@ def read_raster(file, gdf=None, levels=None, expand_px=3):
     return data_dict, tran_dict, is_geo
 
 
-
 def write_raster(in_array, outfile="output.tif", template="somefile.tif"):
     """Write a numpy array to a GeoTIFF file using the geographic transformation
     and projection information from a template raster file.

connectivity/rastio.py

@@ -71,14 +71,14 @@ pub fn multi_level_graph(
             }
         }
 
-        // Process points - this could be parallelized for large datasets
+        // Process points
         for point_idx in 0..num_points {
             let i = i_array[point_idx];
             let j = j_array[point_idx];
             let u = point_idx as u16 + level as u16 * 100;
 
             // Process neighbors at current level
-            process_current_level_neighbors(
+            current_level_neighbors(
                 i, j, u,
                 &COLS, &ROWS,
                 &node_mapping,
@@ -90,7 +90,7 @@ pub fn multi_level_graph(
 
             // Process connections to higher level if needed
             if level < max_level && edge_indices.contains(&(i, j)) {
-                process_higher_level_connections(
+                higher_level_connections(
                     i, j,
                     &node_mapping,
                     &node_mapping_higher,
@@ -108,7 +108,7 @@ pub fn multi_level_graph(
 }
 
 
-/// Create node mapping efficiently
+/// Create node mapping (the unique ID of each node)
 fn create_node_mapping(
     i_array: &[i32],
     j_array: &[i32],
@@ -136,11 +136,10 @@ fn create_node_mapping(
 /// output graph: (u, v) (adj_cond, cond, dist, similarities)
 /// u: source, v: destination
 #[inline]
-fn process_current_level_neighbors(
+fn current_level_neighbors(
     i: i32, 
     j: i32, 
-    u: u16, 
-    // level: i32,
+    u: u16, // source node
     i_ngb: &[i32], 
     j_ngb: &[i32],
     node_mapping: &HashMap<(i32, i32), (u16, f32, Vec<f32>)>,
@@ -170,11 +169,11 @@ fn process_current_level_neighbors(
 }
 
 
-/// Process connections to higher level
+/// Process connections to higher level (e.g. level 2 to level 4)
 /// output graph: (u, v) (adj_cond, cond, dist, similarities)
 /// u: source, v: destination
 #[inline]
-fn process_higher_level_connections(
+fn higher_level_connections(
     i: i32, j: i32,
     node_mapping: &HashMap<(i32, i32), (u16, f32, Vec<f32>)>,
     node_mapping_higher: &HashMap<(i32, i32), (u16, f32, Vec<f32>)>,
@@ -188,7 +187,7 @@ fn process_higher_level_connections(
         let higher_level: i32 = level * 2;
 
         // Get all higher neighbors at once
-        let higher_neighbors = get_higher_neighbors(i, j);
+        let higher_neighbors = get_edge_neighbors(i, j);
 
         // Get the Affines for distance calc
         let transform: &Affine = transforms.get(&level).unwrap();
@@ -214,9 +213,9 @@ fn process_higher_level_connections(
 }
 
 
-/// Get the 3 possible neighbours in the higher level 
+/// Get the 3 possible neighbours of edge cells in the higher level 
 /// i.e. the link between a level edge to its higher level cells
-fn get_higher_neighbors(i: i32, j: i32) -> [(i32, i32); 3] {
+fn get_edge_neighbors(i: i32, j: i32) -> [(i32, i32); 3] {
     // Higher level cell containing the target cell
     let target_higher = (i >> 1, j >> 1);    
     // 8 neighbor offsets: N, S, W, E, NW, NE, SW, SE
@@ -259,7 +258,7 @@ fn get_higher_neighbors(i: i32, j: i32) -> [(i32, i32); 3] {
 }
 
 
-/// Get edge cells efficiently using HashSet
+/// Efficiently get edge cells in a level neighbourhood using HashSet
 fn get_edge_indices(i_arr: &[i32], j_arr: &[i32]) -> HashSet<(i32, i32)> {
     let i_min = *i_arr.iter().min().unwrap_or(&0);
     let i_max = *i_arr.iter().max().unwrap_or(&0);

src/graph.rs

@@ -87,7 +87,7 @@ pub fn multi_level_window(
         (current_height / agg_factor, current_width / agg_factor)
     };
 
-    // Iterate over the NxN neighborhood
+    // Iterate over the NxN neighborhood (this is in potential higher level)
     for di in -radius..=radius {
         for dj in -radius..=radius {
             let higher_i = higher_center_i + di;