nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Cushion isolated diagonal connections of same class cells in a categorial raster

Aline-Wbg

I have a categorical raster I use for a land use map. The categorical raster is further being polygonized and then adjusted in a graphical software for the target map. The raster datas' resolution is borderline for the target map scale, so we do not any generalization before hand and try to retain as much detail as possible.

But, in order to prevent polygons of the same class touching at an edge, we need to remove these cell constellations in the underlying raster.

Boundary Clean with sort type "ascending" removes some of these configurations, while retaining detail, but some diagonal connections remain.

Majority filtering with 8 neighbors is also retaining some of these vertical connections.

I thought about looking at a moving window of 4 cells where "if 2 cells have the same value and are diagonal to each other and the other two cells have another value" it will randomly assign one of the other two cells the value of the two cells being connected over a diagonal, creating a cushion that will not result in dangling connections after polygonization. If the other two cells also have a shared value then it will randomly choose which connection is being cushioned.

Has somebody already implemented something like this? I guess this cushioning could also introduce new isolated diagonal connections...

Find more posts tagged with

Accepted answers

Bsmith_DPW

I think I get it now - Focal Statistics probably won't work then, since the target cell is where we flag these cases, but it's actually an adjacent cell where we would edit the value. This numpy approach is getting closer to what you want. I flag and fix the isolated pixels, then cushion until there are no more diagonal conflicts

import arcpy
import numpy as np

# === CONFIGURATION ===
input_raster          = "Your_Raster_Layer_Name"  # Input layer in ArcGIS Pro
output_cleaned        = "Cushioned_Bridged_Raster" # Final repaired raster
output_isolated_flags = "Isolated_Flags_Raster"   # Output 1: Flag map for Phase 1 (1-pixel islands)
output_bridged_flags  = "Bridged_Flags_Raster"    # Output 2: Flag map for Phase 2 (diagonal bridges)
window_size           = 5                         # Neighborhood size for minority check (odd number)
max_passes            = 5                         # Max passes to clean cascading conflicts
# =====================

arcpy.env.overwriteOutput = True

# Read spatial metadata
desc = arcpy.Describe(input_raster)
cell_size = desc.meanCellWidth
lower_left = desc.extent.lowerLeft
spatial_ref = desc.spatialReference
nodata_val = desc.noDataValue

print("Loading raster into memory...")
r_obj = arcpy.Raster(input_raster)
arr = arcpy.RasterToNumPyArray(r_obj)

rows, cols = arr.shape
radius = window_size // 2

# ---------------------------------------------------------
# PHASE 1: REMOVE & FLAG ISOLATED 1-PIXEL ISLANDS
# ---------------------------------------------------------
print("Phase 1: Cleaning isolated single-pixel islands & recording flags...")
clean_arr = arr.copy()
isolated_flag_arr = np.zeros_like(arr, dtype=np.int8)

for r in range(1, rows - 1):
    for c in range(1, cols - 1):
        val = arr[r, c]
        if val == nodata_val:
            continue
        neighbors = [arr[r-1, c], arr[r+1, c], arr[r, c-1], arr[r, c+1]]
        
        # If cell matches NONE of its 4 orthogonal neighbors
        if all(n != val for n in neighbors if n != nodata_val):
            window = arr[max(0, r-1):min(rows, r+2), max(0, c-1):min(cols, c+2)]
            valid_vals = window[window != nodata_val]
            if len(valid_vals) > 0:
                vals, counts = np.unique(valid_vals, return_counts=True)
                clean_arr[r, c] = vals[np.argmax(counts)]
                isolated_flag_arr[r, c] = 1  # Flag Phase 1 edit

# Save baseline after isolated pixel cleaning
baseline_arr = clean_arr.copy()

# ---------------------------------------------------------
# HELPER FUNCTIONS FOR LOOK-AHEAD PREVENTATIVE CHECKS
# ---------------------------------------------------------
def is_2x2_conflict(p_tl, p_tr, p_bl, p_br):
    """Returns True if a 2x2 grid has an unbridged diagonal conflict."""
    if any(v == nodata_val for v in [p_tl, p_tr, p_bl, p_br]):
        return False
    main_diag = (p_tl == p_br) and (p_tr != p_tl) and (p_bl != p_tl)
    anti_diag = (p_tr == p_bl) and (p_tl != p_tr) and (p_br != p_tr)
    return main_diag or anti_diag

def test_cell_change_causes_conflict(grid, r, c, new_val):
    """Simulates changing grid[r, c] to new_val and checks if it creates a NEW conflict."""
    orig_val = grid[r, c]
    grid[r, c] = new_val
    new_conflict_created = False

    # Check all 4 overlapping 2x2 blocks containing (r, c)
    for dr in [-1, 0]:
        for dc in [-1, 0]:
            r_top, c_left = r + dr, c + dc
            if 0 <= r_top < grid.shape[0] - 1 and 0 <= c_left < grid.shape[1] - 1:
                if is_2x2_conflict(grid[r_top, c_left], grid[r_top, c_left+1],
                                  grid[r_top+1, c_left], grid[r_top+1, c_left+1]):
                    new_conflict_created = True
                    break
        if new_conflict_created:
            break

    grid[r, c] = orig_val  # Revert simulation
    return new_conflict_created

def get_class_count(array, center_r, center_c, target_val):
    """Counts occurrences of target_val in a local neighborhood window."""
    r_min, r_max = max(0, center_r - radius), min(array.shape[0], center_r + radius + 1)
    c_min, c_max = max(0, center_c - radius), min(array.shape[1], center_c + radius + 1)
    sub_grid = array[r_min:r_max, c_min:c_max]
    return np.sum(sub_grid == target_val)

# ---------------------------------------------------------
# PHASE 2: LOOK-AHEAD MINORITY BRIDGING (MULTI-PASS LOOP)
# ---------------------------------------------------------
current_grid = clean_arr.copy()

for pass_num in range(1, max_passes + 1):
    read_arr = current_grid.copy()
    write_arr = current_grid.copy()
    changes_in_pass = 0

    print(f"Phase 2 (Pass {pass_num}/{max_passes}): Scanning and evaluating conflicts...")

    for r in range(rows - 1):
        for c in range(cols - 1):
            p_tl = read_arr[r, c]
            p_tr = read_arr[r, c+1]
            p_bl = read_arr[r+1, c]
            p_br = read_arr[r+1, c+1]

            if any(v == nodata_val for v in [p_tl, p_tr, p_bl, p_br]):
                continue

            candidate_r, candidate_c = None, None

            # --- Case A: Main Diagonal Conflict (Top-Left == Bottom-Right) ---
            if (p_tl == p_br) and (p_tr != p_tl) and (p_bl != p_tl):
                target_color = p_tl
                
                tr_creates_conflict = test_cell_change_causes_conflict(read_arr, r, c+1, target_color)
                bl_creates_conflict = test_cell_change_causes_conflict(read_arr, r+1, c, target_color)

                if not tr_creates_conflict and bl_creates_conflict:
                    candidate_r, candidate_c = r, c+1
                elif not bl_creates_conflict and tr_creates_conflict:
                    candidate_r, candidate_c = r+1, c
                else:
                    cnt_tr = get_class_count(read_arr, r, c+1, p_tr)
                    cnt_bl = get_class_count(read_arr, r+1, c, p_bl)
                    if cnt_tr <= cnt_bl:
                        candidate_r, candidate_c = r, c+1
                    else:
                        candidate_r, candidate_c = r+1, c

                write_arr[candidate_r, candidate_c] = target_color
                changes_in_pass += 1

            # --- Case B: Anti-Diagonal Conflict (Top-Right == Bottom-Left) ---
            elif (p_tr == p_bl) and (p_tl != p_tr) and (p_br != p_tr):
                target_color = p_tr
                
                tl_creates_conflict = test_cell_change_causes_conflict(read_arr, r, c, target_color)
                br_creates_conflict = test_cell_change_causes_conflict(read_arr, r+1, c+1, target_color)

                if not tl_creates_conflict and br_creates_conflict:
                    candidate_r, candidate_c = r, c
                elif not br_creates_conflict and tl_creates_conflict:
                    candidate_r, candidate_c = r+1, c+1
                else:
                    cnt_tl = get_class_count(read_arr, r, c, p_tl)
                    cnt_br = get_class_count(read_arr, r+1, c+1, p_br)
                    if cnt_tl <= cnt_br:
                        candidate_r, candidate_c = r, c
                    else:
                        candidate_r, candidate_c = r+1, c+1

                write_arr[candidate_r, candidate_c] = target_color
                changes_in_pass += 1

    current_grid = write_arr.copy()
    print(f"   Pass {pass_num} complete. Conflicts resolved/bridged: {changes_in_pass}")

    if changes_in_pass == 0:
        print("   Topology converged! No further diagonal conflicts remain.")
        break

# ---------------------------------------------------------
# PHASE 3: DERIVE FLAGS & EXPORT TO ARCGIS PRO
# ---------------------------------------------------------
print("Writing rasters back to disk...")

# Phase 2 Bridged Flag (differs from baseline)
bridged_flag_arr = np.where(current_grid != baseline_arr, 1, 0).astype(np.int8)

# Output 1: Phase 1 Isolated Flags
out_iso_flags = arcpy.NumPyArrayToRaster(isolated_flag_arr, lower_left, cell_size, cell_size, 0)
arcpy.DefineProjection_management(out_iso_flags, spatial_ref)
out_iso_flags.save(output_isolated_flags)

# Output 2: Phase 2 Bridged Flags
out_bridged_flags = arcpy.NumPyArrayToRaster(bridged_flag_arr, lower_left, cell_size, cell_size, 0)
arcpy.DefineProjection_management(out_bridged_flags, spatial_ref)
out_bridged_flags.save(output_bridged_flags)

# Output 3: Cleaned & Cushioned Raster
out_raster = arcpy.NumPyArrayToRaster(current_grid, lower_left, cell_size, cell_size, nodata_val)
arcpy.DefineProjection_management(out_raster, spatial_ref)
out_raster.save(output_cleaned)

print(f"\nSuccess! Generated 3 output layers:\n"
      f"1. Phase 1 Isolated Flags:  {output_isolated_flags}\n"
      f"2. Phase 2 Bridged Flags:   {output_bridged_flags}\n"
      f"3. Final Cleaned Raster:    {output_cleaned}")

All comments

DanPatterson

The strategies for trying to clean up rasters are outlined, as you have employed some.

Generalization of classified raster imagery | ArcGIS Pro documentation

You might want to consider a regiongroup, extract by attributes (based on some threshold size), then nibble to fill in the areas. This may affect other areas that don't meet at a corner, but you will have to assess.

If it is one particular class you could even attempt a shrink/expand.

Even in your first image, there are some conditions that were missed when doing a manual inspection.

Nothing will be perfect, but employing the closest will leave less vector cleanup.

good luck

Aline-Wbg

Hi Dan thank you a lot for your response. We actually do employ the workflow you suggested for small-scale maps.

For our large-scale maps we leaned away from applying it as too much detail (that we need) was being removed. So in this case we directly polygonize and then smooth it in a special software that retains topology and area. Here we really would like to just "pad" corners where polygons of the same class touch but the two other cells in a 2x2 configuration are in another class before we proceed to generalization.....(We also remove single-pixel polygons, these little crumbs).

The full raster generalization workflow (which I agree with you on is generally better than doing the clean-up in vector form) actually does remove a lot of those diagonal situations (not all of them) but also detail (which should be avoided in our case for these specific maps).

Do you have an idea how to find those pixel configurations and change the neighbouring ones according to my sketch?

And

DanPatterson

You could explore the concept by using a 2x2 kernel with a focal "variety"

Focal Statistics (Spatial Analyst Tools) | ArcGIS Pro documentation

if the variety is 2 or 3 would identify those locations but you would still be tasked with how examine those diagonals with a weight kernel, of which you have 2

0 1

1 0

0 1

I see nothing that would handle this situation in one pass

Aline-Wbg

Thank you Dan! Yes it is tricky....

Bsmith_DPW

An interesting problem. I had the same initial thought as Dan, with two different 2x2 kernels. But what if instead we evaluate a 3x3 window, and define a diagonal conflict as "matches a diagonal neighbor but does not match either of the orthogonal neighbors". E.g., matches top left, but neither top nor left. This script might do the trick (thanks Gemini!). It assigns values based on the majority in a 5x5 region but you can set that to whatever you'd like. I don't know if you'll need to run it more than once to resolve new issues introduced in the correction itself.

import os
import arcpy
from arcpy.sa import *

# === CONFIGURATION ===
input_raster = "Your_Raster_Layer_Name"  # Change to your layer name in ArcGIS Pro
output_flags = "Conflict_Flags_Raster"   # Output 1: 1 = Conflict found, 0 = Clean
output_fixed = "Fixed_Topology_Raster"   # Output 2: The final repaired raster
# =====================

# Check out the Spatial Analyst extension
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True
scratch_dir = arcpy.env.scratchFolder

# Helper function to generate 3x3 single-neighbor kernel text files
def make_kernel(row, col, name):
    path = os.path.join(scratch_dir, f"{name}.txt")
    matrix = [["0","0","0"], ["0","0","0"], ["0","0","0"]]
    matrix[row][col] = "1"
    with open(path, "w") as f:
        f.write("3 3\n")
        for r in matrix:
            f.write(" ".join(r) + "\n")
    return path

print("Analyzing neighborhood structures...")
# Isolate orthogonal neighbors
r_T = FocalStatistics(input_raster, NbrIrregular(make_kernel(0, 1, "k_T")), "MAXIMUM")
r_B = FocalStatistics(input_raster, NbrIrregular(make_kernel(2, 1, "k_B")), "MAXIMUM")
r_L = FocalStatistics(input_raster, NbrIrregular(make_kernel(1, 0, "k_L")), "MAXIMUM")
r_R = FocalStatistics(input_raster, NbrIrregular(make_kernel(1, 2, "k_R")), "MAXIMUM")

# Isolate diagonal neighbors
r_TL = FocalStatistics(input_raster, NbrIrregular(make_kernel(0, 0, "k_TL")), "MAXIMUM")
r_TR = FocalStatistics(input_raster, NbrIrregular(make_kernel(0, 2, "k_TR")), "MAXIMUM")
r_BL = FocalStatistics(input_raster, NbrIrregular(make_kernel(2, 0, "k_BL")), "MAXIMUM")
r_BR = FocalStatistics(input_raster, NbrIrregular(make_kernel(2, 2, "k_BR")), "MAXIMUM")

# Define the 4 conflict conditions relative to the center cell
print("Identifying diagonal conflicts...")
r = Raster(input_raster)
conflict_condition = (
    ((r == r_BR) & (r != r_R) & (r != r_B)) |  # Bottom-Right conflict
    ((r == r_BL) & (r != r_L) & (r != r_B)) |  # Bottom-Left conflict
    ((r == r_TR) & (r != r_R) & (r != r_T)) |  # Top-Right conflict
    ((r == r_TL) & (r != r_L) & (r != r_T))    # Top-Left conflict
)

# OUTPUT 1: Generate and save the Flag Raster
print("Saving conflict flag map...")
flag_raster = Con(conflict_condition, 1, 0)
flag_raster.save(output_flags)

# Calculate the replacement raster
print("Calculating majority replacement map...")
maj_ras = FocalStatistics(input_raster, NbrRectangle(5, 5, "CELL"), "MAJORITY")

# OUTPUT 2: Generate and save the Fixed Raster
print("Saving fixed topology raster...")
fixed_raster = Con(conflict_condition, maj_ras, r)
fixed_raster.save(output_fixed)

print(f"\nSuccess!")

Aline-Wbg

Hi @BrennanSmith1 thank you for your reply and interesting idea.

In this picture you can see two cells that would be highlighted as a cell being in conflict. If we replace them with a neighborhood majority (here 3x3) the two red zones will be disconnected and there could be the possibility that the lower red zone of which we only see two pixels here is actually really big and it would make sense to keep the connection between those two zones and cushion their connection (replacing either the yellow pixel or the green pixel with red in their 2x2 neighborhood).

With your / an Geminis 🙂 approach we get rid of these connections..

It would be really nice to only cushion these connections and not getting rid of them.

PS: What I forgot to mention: I get rid of isolated pixels (pixels that potentially only have diagonal connections to same-class pixels or are fully isolated) beforehand as I am not interested to cushion connections where one zone is just one pixel....

Bsmith_DPW

import arcpy
import numpy as np

# === CONFIGURATION ===
input_raster          = "Your_Raster_Layer_Name"  # Input layer in ArcGIS Pro
output_cleaned        = "Cushioned_Bridged_Raster" # Final repaired raster
output_isolated_flags = "Isolated_Flags_Raster"   # Output 1: Flag map for Phase 1 (1-pixel islands)
output_bridged_flags  = "Bridged_Flags_Raster"    # Output 2: Flag map for Phase 2 (diagonal bridges)
window_size           = 5                         # Neighborhood size for minority check (odd number)
max_passes            = 5                         # Max passes to clean cascading conflicts
# =====================

arcpy.env.overwriteOutput = True

# Read spatial metadata
desc = arcpy.Describe(input_raster)
cell_size = desc.meanCellWidth
lower_left = desc.extent.lowerLeft
spatial_ref = desc.spatialReference
nodata_val = desc.noDataValue

print("Loading raster into memory...")
r_obj = arcpy.Raster(input_raster)
arr = arcpy.RasterToNumPyArray(r_obj)

rows, cols = arr.shape
radius = window_size // 2

# ---------------------------------------------------------
# PHASE 1: REMOVE & FLAG ISOLATED 1-PIXEL ISLANDS
# ---------------------------------------------------------
print("Phase 1: Cleaning isolated single-pixel islands & recording flags...")
clean_arr = arr.copy()
isolated_flag_arr = np.zeros_like(arr, dtype=np.int8)

for r in range(1, rows - 1):
    for c in range(1, cols - 1):
        val = arr[r, c]
        if val == nodata_val:
            continue
        neighbors = [arr[r-1, c], arr[r+1, c], arr[r, c-1], arr[r, c+1]]
        
        # If cell matches NONE of its 4 orthogonal neighbors
        if all(n != val for n in neighbors if n != nodata_val):
            window = arr[max(0, r-1):min(rows, r+2), max(0, c-1):min(cols, c+2)]
            valid_vals = window[window != nodata_val]
            if len(valid_vals) > 0:
                vals, counts = np.unique(valid_vals, return_counts=True)
                clean_arr[r, c] = vals[np.argmax(counts)]
                isolated_flag_arr[r, c] = 1  # Flag Phase 1 edit

# Save baseline after isolated pixel cleaning
baseline_arr = clean_arr.copy()

# ---------------------------------------------------------
# HELPER FUNCTIONS FOR LOOK-AHEAD PREVENTATIVE CHECKS
# ---------------------------------------------------------
def is_2x2_conflict(p_tl, p_tr, p_bl, p_br):
    """Returns True if a 2x2 grid has an unbridged diagonal conflict."""
    if any(v == nodata_val for v in [p_tl, p_tr, p_bl, p_br]):
        return False
    main_diag = (p_tl == p_br) and (p_tr != p_tl) and (p_bl != p_tl)
    anti_diag = (p_tr == p_bl) and (p_tl != p_tr) and (p_br != p_tr)
    return main_diag or anti_diag

def test_cell_change_causes_conflict(grid, r, c, new_val):
    """Simulates changing grid[r, c] to new_val and checks if it creates a NEW conflict."""
    orig_val = grid[r, c]
    grid[r, c] = new_val
    new_conflict_created = False

    # Check all 4 overlapping 2x2 blocks containing (r, c)
    for dr in [-1, 0]:
        for dc in [-1, 0]:
            r_top, c_left = r + dr, c + dc
            if 0 <= r_top < grid.shape[0] - 1 and 0 <= c_left < grid.shape[1] - 1:
                if is_2x2_conflict(grid[r_top, c_left], grid[r_top, c_left+1],
                                  grid[r_top+1, c_left], grid[r_top+1, c_left+1]):
                    new_conflict_created = True
                    break
        if new_conflict_created:
            break

    grid[r, c] = orig_val  # Revert simulation
    return new_conflict_created

def get_class_count(array, center_r, center_c, target_val):
    """Counts occurrences of target_val in a local neighborhood window."""
    r_min, r_max = max(0, center_r - radius), min(array.shape[0], center_r + radius + 1)
    c_min, c_max = max(0, center_c - radius), min(array.shape[1], center_c + radius + 1)
    sub_grid = array[r_min:r_max, c_min:c_max]
    return np.sum(sub_grid == target_val)

# ---------------------------------------------------------
# PHASE 2: LOOK-AHEAD MINORITY BRIDGING (MULTI-PASS LOOP)
# ---------------------------------------------------------
current_grid = clean_arr.copy()

for pass_num in range(1, max_passes + 1):
    read_arr = current_grid.copy()
    write_arr = current_grid.copy()
    changes_in_pass = 0

    print(f"Phase 2 (Pass {pass_num}/{max_passes}): Scanning and evaluating conflicts...")

    for r in range(rows - 1):
        for c in range(cols - 1):
            p_tl = read_arr[r, c]
            p_tr = read_arr[r, c+1]
            p_bl = read_arr[r+1, c]
            p_br = read_arr[r+1, c+1]

            if any(v == nodata_val for v in [p_tl, p_tr, p_bl, p_br]):
                continue

            candidate_r, candidate_c = None, None

            # --- Case A: Main Diagonal Conflict (Top-Left == Bottom-Right) ---
            if (p_tl == p_br) and (p_tr != p_tl) and (p_bl != p_tl):
                target_color = p_tl
                
                tr_creates_conflict = test_cell_change_causes_conflict(read_arr, r, c+1, target_color)
                bl_creates_conflict = test_cell_change_causes_conflict(read_arr, r+1, c, target_color)

                if not tr_creates_conflict and bl_creates_conflict:
                    candidate_r, candidate_c = r, c+1
                elif not bl_creates_conflict and tr_creates_conflict:
                    candidate_r, candidate_c = r+1, c
                else:
                    cnt_tr = get_class_count(read_arr, r, c+1, p_tr)
                    cnt_bl = get_class_count(read_arr, r+1, c, p_bl)
                    if cnt_tr <= cnt_bl:
                        candidate_r, candidate_c = r, c+1
                    else:
                        candidate_r, candidate_c = r+1, c

                write_arr[candidate_r, candidate_c] = target_color
                changes_in_pass += 1

            # --- Case B: Anti-Diagonal Conflict (Top-Right == Bottom-Left) ---
            elif (p_tr == p_bl) and (p_tl != p_tr) and (p_br != p_tr):
                target_color = p_tr
                
                tl_creates_conflict = test_cell_change_causes_conflict(read_arr, r, c, target_color)
                br_creates_conflict = test_cell_change_causes_conflict(read_arr, r+1, c+1, target_color)

                if not tl_creates_conflict and br_creates_conflict:
                    candidate_r, candidate_c = r, c
                elif not br_creates_conflict and tl_creates_conflict:
                    candidate_r, candidate_c = r+1, c+1
                else:
                    cnt_tl = get_class_count(read_arr, r, c, p_tl)
                    cnt_br = get_class_count(read_arr, r+1, c+1, p_br)
                    if cnt_tl <= cnt_br:
                        candidate_r, candidate_c = r, c
                    else:
                        candidate_r, candidate_c = r+1, c+1

                write_arr[candidate_r, candidate_c] = target_color
                changes_in_pass += 1

    current_grid = write_arr.copy()
    print(f"   Pass {pass_num} complete. Conflicts resolved/bridged: {changes_in_pass}")

    if changes_in_pass == 0:
        print("   Topology converged! No further diagonal conflicts remain.")
        break

# ---------------------------------------------------------
# PHASE 3: DERIVE FLAGS & EXPORT TO ARCGIS PRO
# ---------------------------------------------------------
print("Writing rasters back to disk...")

# Phase 2 Bridged Flag (differs from baseline)
bridged_flag_arr = np.where(current_grid != baseline_arr, 1, 0).astype(np.int8)

# Output 1: Phase 1 Isolated Flags
out_iso_flags = arcpy.NumPyArrayToRaster(isolated_flag_arr, lower_left, cell_size, cell_size, 0)
arcpy.DefineProjection_management(out_iso_flags, spatial_ref)
out_iso_flags.save(output_isolated_flags)

# Output 2: Phase 2 Bridged Flags
out_bridged_flags = arcpy.NumPyArrayToRaster(bridged_flag_arr, lower_left, cell_size, cell_size, 0)
arcpy.DefineProjection_management(out_bridged_flags, spatial_ref)
out_bridged_flags.save(output_bridged_flags)

# Output 3: Cleaned & Cushioned Raster
out_raster = arcpy.NumPyArrayToRaster(current_grid, lower_left, cell_size, cell_size, nodata_val)
arcpy.DefineProjection_management(out_raster, spatial_ref)
out_raster.save(output_cleaned)

print(f"\nSuccess! Generated 3 output layers:\n"
      f"1. Phase 1 Isolated Flags:  {output_isolated_flags}\n"
      f"2. Phase 2 Bridged Flags:   {output_bridged_flags}\n"
      f"3. Final Cleaned Raster:    {output_cleaned}")

Aline-Wbg

This worked! Well done. The only thing that can happen is that new solo pixels are introduced so I made a loop where it runs as many times up until no new flags are added anymore and also if there's two isolated pixels next to each other instead of flipping back and forth (exchanging values) it will prevent that. The code also removes solo pixels that are at the raster border.

import arcpy
import numpy as np
import os

# === CONFIGURATION ===
input_raster          =  

# Use ArcGIS Pro's physical scratch folder so it saves as proper .tif files without GDB errors
scratch_dir           = arcpy.env.scratchFolder

output_cleaned        = os.path.join(scratch_dir, "Cushioned_Bridged_Raster.tif")
output_isolated_cleaned = os.path.join(scratch_dir, "Phase1_Isolated_Cleaned_Raster.tif")
output_isolated_flags = os.path.join(scratch_dir, "Isolated_Flags_Raster.tif")
output_bridged_flags  = os.path.join(scratch_dir, "Bridged_Flags_Raster.tif")

window_size           = 5                         
max_macro_iterations  = 10                        # Safety limit for full cleaning cycles
max_passes            = 5                         # Max passes for diagonal bridging per cycle
# =====================

arcpy.env.overwriteOutput = True

# Read spatial metadata
desc = arcpy.Describe(input_raster)
cell_size = desc.meanCellWidth
lower_left = desc.extent.lowerLeft
spatial_ref = desc.spatialReference
nodata_val = desc.noDataValue

print("Loading raster into memory...")
r_obj = arcpy.Raster(input_raster)
arr = arcpy.RasterToNumPyArray(r_obj)
rows, cols = arr.shape
radius = window_size // 2

current_grid = arr.copy()
total_isolated_flags = np.zeros_like(arr, dtype=np.int8)
total_bridged_flags = np.zeros_like(arr, dtype=np.int8)

# ---------------------------------------------------------
# HELPER FUNCTIONS 
# ---------------------------------------------------------
def is_2x2_conflict(p_tl, p_tr, p_bl, p_br):
    if all(v == nodata_val for v in [p_tl, p_tr, p_bl, p_br]):
        return False
    if p_tl == p_tr == p_bl == p_br:
        return False

    main_diag = (p_tl != nodata_val) and (p_tl == p_br) and (p_tr != p_tl) and (p_bl != p_tl)
    anti_diag = (p_tr != nodata_val) and (p_tr == p_bl) and (p_tl != p_tr) and (p_br != p_tr)
    
    return main_diag or anti_diag

def test_cell_change_causes_conflict(grid, r_idx, c_idx, new_val):
    orig_val = grid[r_idx, c_idx]
    grid[r_idx, c_idx] = new_val
    new_conflict_created = False

    for dr in [-1, 0]:
        for dc in [-1, 0]:
            r_top, c_left = r_idx + dr, c_idx + dc
            if 0 <= r_top < grid.shape[0] - 1 and 0 <= c_left < grid.shape[1] - 1:
                if is_2x2_conflict(grid[r_top, c_left], grid[r_top, c_left+1],
                                  grid[r_top+1, c_left], grid[r_top+1, c_left+1]):
                    new_conflict_created = True
                    break
        if new_conflict_created:
            break

    grid[r_idx, c_idx] = orig_val
    return new_conflict_created

def get_class_count(array, center_r, center_c, target_val):
    r_min, r_max = max(0, center_r - radius), min(array.shape[0], center_r + radius + 1)
    c_min, c_max = max(0, center_c - radius), min(array.shape[1], center_c + radius + 1)
    sub_grid = array[r_min:r_max, c_min:c_max]
    return np.sum((sub_grid == target_val) & (sub_grid != nodata_val))

# ---------------------------------------------------------
# MASTER ITERATIVE CONVERGENCE LOOP
# ---------------------------------------------------------
print("Starting iterative convergence loop...")

for macro_iter in range(1, max_macro_iterations + 1):
    print(f"\n--- Macro Iteration {macro_iter} ---")
    grid_before_cycle = current_grid.copy()
    
    # -----------------------------------------------------
    # PHASE 1: REMOVE & FLAG ISOLATED PIXELS (Edge-Inclusive)
    # -----------------------------------------------------
    mask_raster_arr = np.zeros_like(arr, dtype=np.int8)
    
    for r_idx in range(rows):
        for c_idx in range(cols):
            center = current_grid[r_idx, c_idx]
            if center == nodata_val:
                continue
                
            top    = current_grid[r_idx-1, c_idx] if r_idx > 0 else nodata_val
            bottom = current_grid[r_idx+1, c_idx] if r_idx < rows - 1 else nodata_val
            left   = current_grid[r_idx, c_idx-1] if c_idx > 0 else nodata_val
            right  = current_grid[r_idx, c_idx+1] if c_idx < cols - 1 else nodata_val
            
            neighbors = [top, bottom, left, right]
            valid_neighbors = [n for n in neighbors if n != nodata_val]
            
            if len(valid_neighbors) == 0:
                continue
                
            if all(n != center for n in valid_neighbors):
                has_isolated_partner = False
                for dr, dc in [(-1,0), (1,0), (0,-1), (0,1)]:
                    nr, nc = r_idx + dr, c_idx + dc
                    if 0 <= nr < rows and 0 <= nc < cols and current_grid[nr, nc] == center:
                        n_top    = current_grid[nr-1, nc] if nr > 0 else nodata_val
                        n_bottom = current_grid[nr+1, nc] if nr < rows - 1 else nodata_val
                        n_left   = current_grid[nr, nc-1] if nc > 0 else nodata_val
                        n_right  = current_grid[nr, nc+1] if nc < cols - 1 else nodata_val
                        
                        n_valid = [nn for nn in [n_top, n_bottom, n_left, n_right] if nn != nodata_val]
                        if len(n_valid) > 0 and all(nn != center for nn in n_valid):
                            has_isolated_partner = True
                            break
                
                if not has_isolated_partner:
                    mask_raster_arr[r_idx, c_idx] = 1

    stable_arr = np.where(mask_raster_arr == 1, nodata_val, current_grid)

    filled_arr = stable_arr.copy()
    iter_iso_flags = np.zeros_like(arr, dtype=np.int8)
    iso_changes = 0

    for r_idx in range(rows):
        for c_idx in range(cols):
            if mask_raster_arr[r_idx, c_idx] == 1:
                top    = stable_arr[r_idx-1, c_idx] if r_idx > 0 else nodata_val
                bottom = stable_arr[r_idx+1, c_idx] if r_idx < rows - 1 else nodata_val
                left   = stable_arr[r_idx, c_idx-1] if c_idx > 0 else nodata_val
                right  = stable_arr[r_idx, c_idx+1] if c_idx < cols - 1 else nodata_val
                
                clean_ortho = [o for o in [top, bottom, left, right] if o != nodata_val]
                
                if len(clean_ortho) > 0:
                    vals, counts = np.unique(clean_ortho, return_counts=True)
                    filled_arr[r_idx, c_idx] = vals[np.argmax(counts)]
                    iter_iso_flags[r_idx, c_idx] = 1
                    iso_changes += 1

    clean_arr = filled_arr.copy()
    baseline_arr = clean_arr.copy()

    # -----------------------------------------------------
    # PHASE 2: LOOK-AHEAD MINORITY BRIDGING (NoData-Inclusive)
    # -----------------------------------------------------
    pass_grid = clean_arr.copy()
    bridged_changes_in_cycle = 0

    for pass_num in range(1, max_passes + 1):
        read_arr = pass_grid.copy()
        write_arr = pass_grid.copy()
        changes_in_pass = 0

        for r_idx in range(rows - 1):
            for c_idx in range(cols - 1):
                p_tl = read_arr[r_idx, c_idx]
                p_tr = read_arr[r_idx, c_idx+1]
                p_bl = read_arr[r_idx+1, c_idx]
                p_br = read_arr[r_idx+1, c_idx+1]

                if all(v == nodata_val for v in [p_tl, p_tr, p_bl, p_br]):
                    continue

                candidate_r, candidate_c = None, None

                # Case A: Main Diagonal Conflict (TL == BR)
                if (p_tl != nodata_val) and (p_tl == p_br) and (p_tr != p_tl) and (p_bl != p_tl):
                    target_color = p_tl
                    tr_cc = test_cell_change_causes_conflict(read_arr, r_idx, c_idx+1, target_color)
                    bl_cc = test_cell_change_causes_conflict(read_arr, r_idx+1, c_idx, target_color)

                    if not tr_cc and bl_cc:
                        candidate_r, candidate_c = r_idx, c_idx+1
                    elif not bl_cc and tr_cc:
                        candidate_r, candidate_c = r_idx+1, c_idx
                    else:
                        if p_tr == nodata_val:
                            candidate_r, candidate_c = r_idx, c_idx+1
                        elif p_bl == nodata_val:
                            candidate_r, candidate_c = r_idx+1, c_idx
                        else:
                            cnt_tr = get_class_count(read_arr, r_idx, c_idx+1, p_tr)
                            cnt_bl = get_class_count(read_arr, r_idx+1, c_idx, p_bl)
                            candidate_r, candidate_c = (r_idx, c_idx+1) if cnt_tr <= cnt_bl else (r_idx+1, c_idx)

                    write_arr[candidate_r, candidate_c] = target_color
                    changes_in_pass += 1

                # Case B: Anti-Diagonal Conflict (TR == BL)
                elif (p_tr != nodata_val) and (p_tr == p_bl) and (p_tl != p_tr) and (p_br != p_tr):
                    target_color = p_tr
                    tl_cc = test_cell_change_causes_conflict(read_arr, r_idx, c_idx, target_color)
                    br_cc = test_cell_change_causes_conflict(read_arr, r_idx+1, c_idx+1, target_color)

                    if not tl_cc and br_cc:
                        candidate_r, candidate_c = r_idx, c_idx
                    elif not br_cc and tl_cc:
                        candidate_r, candidate_c = r_idx+1, c_idx+1
                    else:
                        if p_tl == nodata_val:
                            candidate_r, candidate_c = r_idx, c_idx
                        elif p_br == nodata_val:
                            candidate_r, candidate_c = r_idx+1, c_idx+1
                        else:
                            cnt_tl = get_class_count(read_arr, r_idx, c_idx, p_tl)
                            cnt_br = get_class_count(read_arr, r_idx+1, c_idx+1, p_br)
                            candidate_r, candidate_c = (r_idx, c_idx) if cnt_tl <= cnt_br else (r_idx+1, c_idx+1)

                    write_arr[candidate_r, candidate_c] = target_color
                    changes_in_pass += 1

        pass_grid = write_arr.copy()
        bridged_changes_in_cycle += changes_in_pass
        if changes_in_pass == 0:
            break

    current_grid = pass_grid.copy()
    
    total_isolated_flags = np.where((iter_iso_flags == 1) | (total_isolated_flags == 1), 1, 0).astype(np.int8)
    cycle_bridged_flags = np.where((current_grid != baseline_arr), 1, 0).astype(np.int8)
    total_bridged_flags = np.where((cycle_bridged_flags == 1) | (total_bridged_flags == 1), 1, 0).astype(np.int8)

    print(f"   Isolated pixels fixed: {iso_changes} | Diagonal bridges fixed: {bridged_changes_in_cycle}")

    if iso_changes == 0 and bridged_changes_in_cycle == 0:
        print("\nFull structural convergence reached! Both flag categories are empty.")
        break

# ---------------------------------------------------------
# EXPORT TO SCRATCH FOLDER & ADD ALL LAYERS TO MAP
# ---------------------------------------------------------
print("Exporting rasters to scratch folder and pushing to active map...")

out_iso_cleaned_raster = arcpy.NumPyArrayToRaster(clean_arr, lower_left, cell_size, cell_size, nodata_val)
arcpy.DefineProjection_management(out_iso_cleaned_raster, spatial_ref)
out_iso_cleaned_raster.save(output_isolated_cleaned)

out_iso_flags = arcpy.NumPyArrayToRaster(total_isolated_flags, lower_left, cell_size, cell_size, 0)
arcpy.DefineProjection_management(out_iso_flags, spatial_ref)
out_iso_flags.save(output_isolated_flags)

out_bridged_flags = arcpy.NumPyArrayToRaster(total_bridged_flags, lower_left, cell_size, cell_size, 0)
arcpy.DefineProjection_management(out_bridged_flags, spatial_ref)
out_bridged_flags.save(output_bridged_flags)

out_raster = arcpy.NumPyArrayToRaster(current_grid, lower_left, cell_size, cell_size, nodata_val)
arcpy.DefineProjection_management(out_raster, spatial_ref)
out_raster.save(output_cleaned)

p = arcpy.mp.ArcGISProject("CURRENT")
m = p.activeMap
if m:
    m.addDataFromPath(output_cleaned)
    m.addDataFromPath(output_isolated_cleaned)
    m.addDataFromPath(output_isolated_flags)
    m.addDataFromPath(output_bridged_flags)

print(f"\nSuccess! All scratch rasters generated and added to active map view.")

Join the Conversation

Signed in members can post, follow updates, and more. New here? Register a free account.

Quick Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.