from geosoft import gxpy import geosoft.gxpy.gdb as gxdb import numpy as np import pandas as pd from scipy.spatial import cKDTree import math import time import os import warnings from itertools import combinations # Ignore warnings that may arise from NumPy operations on NaN/dummy values warnings.filterwarnings("ignore", category=UserWarning) warnings.filterwarnings("ignore", category=RuntimeWarning) # Initialize Geosoft context once try: gxc = gxpy.gx.GXpy() except Exception as e: print(f"Warning: Could not initialize Geosoft context: {e}") # ----------------------------- # Configuration (USER MUST EDIT) # ----------------------------- # !!! UPDATE THESE PATHS !!! # This path must be correct for the script to run Stage 1 and 3 GDB_PATH = r"C:\Users\NivedithaPoolla\Documents\P110-2025\Test Oasis Project\Database\P285-2025_MAG_Channel_Turning_Processed_Data_Split_A1.gdb" WORKING_DIR = r"C:\Users\NivedithaPoolla\Documents\P110-2025\Test Oasis Project\Database" TEMP_CSV_INPUT = os.path.join(WORKING_DIR, "temp_gdb_export_input.csv") # Channels required for the processing logic CHECK_COLUMN = "Signal" ALT_COLUMN = "Alt" ESSENTIAL_CHANNELS = ["X", "Y", CHECK_COLUMN, ALT_COLUMN] # KD-Tree parameters RADIUS_M = 0.5 # Reverting to 0.5m, as it worked in the collapsed CSV script PARALLEL_TOLERANCE_DEG = 10.0 FLAG_CHANNELS = [ "is_intersection", "is_minimum_signal_line", "is_maximum_alt_line" ] # CRITICAL: Tolerance for robust floating-point comparison FLOAT_TOLERANCE = 1e-6 # ----------------------------- # Helper Functions (Endpoint Angle Check) # ----------------------------- def get_endpoint_angle(coords): """Calculates the angle based only on the start and end points of the segment.""" if len(coords) < 2: return np.nan dx = coords[-1][0] - coords[0][0] dy = coords[-1][1] - coords[0][1] if (dx == 0 and dy == 0): return np.nan return math.atan2(dy, dx) def are_parallel(coords1, coords2, tolerance_degrees=30.0): """Checks if two segments are approximately parallel using their endpoint angles.""" angle1 = get_endpoint_angle(coords1) angle2 = get_endpoint_angle(coords2) if np.isnan(angle1) or np.isnan(angle2): return False diff_rad = abs(angle1 - angle2) % math.pi if diff_rad > math.pi / 2: diff_rad = math.pi - diff_rad diff_deg = math.degrees(diff_rad) return diff_deg <= tolerance_degrees def flatten(arr): """Flatten nested lists/arrays from Geosoft VV channels into 1d numpy float array.""" out = [] for item in arr: if item is None: continue if isinstance(item, (list, np.ndarray)): out.extend(flatten(item)) else: try: out.append(float(item)) except: pass return np.array(out, dtype=float) # -------------------------------------------------------- # STAGE 1: EXPORT GDB TO CSV # -------------------------------------------------------- def export_gdb_to_csv(gdb_path, csv_path, essential_channels): """Exports all lines and channels from a GDB to a single CSV file, including NaNs.""" t_start = time.time() print(f"\n--- STAGE 1: EXPORTING GDB to CSV ---") # ... (GDB Export logic as before, ensuring NaNs are included in the export) try: with gxdb.Geosoft_gdb.open(gdb_path) as gdb_file: all_channels = gdb_file.list_channels() lines = gdb_file.list_lines() print(f"Found {len(lines)} lines in GDB.") export_channels = [ch for ch in all_channels if ch not in ['Line', 'row_id']] records = [] global_row_id = 0 for idx, line in enumerate(lines, start=1): gdb_file.current_line = line if not all(c in all_channels for c in essential_channels): continue try: line_data = {ch: flatten(gdb_file.read_channel(line, ch)) for ch in export_channels} n = min(len(arr) for arr in line_data.values()) if line_data else 0 if n == 0: continue for i in range(n): record = {"row_id": int(global_row_id), "Line": line} for ch in export_channels: record[ch] = line_data[ch][i] records.append(record) global_row_id += 1 except Exception as e: print(f"[{idx}/{len(lines)}] Error reading line {line}: {e}") continue except Exception as e: print(f"FATAL GDB Export Error: {e}") return None if not records: print("No valid data exported.") return None df_export = pd.DataFrame.from_records(records) if 'X' not in df_export.columns or 'Y' not in df_export.columns: raise ValueError("Export failed: X and Y channels missing.") os.makedirs(os.path.dirname(csv_path) or '.', exist_ok=True) df_export.to_csv(csv_path, index=False) t_end = time.time() print(f"✔ Export successful. {len(df_export)} points written to: {csv_path}") print(f"Stage 1 Time: {t_end - t_start:.2f} seconds.") return csv_path # -------------------------------------------------------- # STAGE 2: CSV PROCESSING (Simple Geometry, Robust Flags) # -------------------------------------------------------- def mark_intersections_kdtree(input_csv_path, radius, parallel_tolerance_deg, output_path): """ Identifies intersection points using cKDTree by collapsing lines (simple geometry) and applying robust, tie-breaking flag attribution. """ t_start = time.time() print(f"\n--- STAGE 2: PROCESSING CSV with KD-Tree (Simple Geometry) ---") SIGNAL_COL = CHECK_COLUMN ALT_COL = ALT_COLUMN FLAG_COLUMN_INT = FLAG_CHANNELS[0] FLAG_SIG = FLAG_CHANNELS[1] FLAG_ALT = FLAG_CHANNELS[2] # 1. Load the full data (for final flagging) df_with_dummies = pd.read_csv(input_csv_path) df_with_dummies['row_id'] = df_with_dummies['row_id'].astype(int) # 2. Prepare the collapsed data (matching the successful CSV script) # This prepares the data for the KD-Tree search. df_valid = df_with_dummies.dropna(subset=["X", "Y"]).copy().reset_index(drop=True) # We must use the original row_id from df_with_dummies for final flagging, # but the temporary index of df_valid for the KD-Tree coordinates. df_valid['original_row_id'] = df_valid['row_id'].astype(int) line_ids = df_valid["Line"].unique().tolist() print(f"{len(line_ids)} unique lines detected in cleaned data.") line_coords = {} line_info = {} # Stores original row_id # Pre-process coordinates and row info by line for line in line_ids: group = df_valid[df_valid["Line"] == line].sort_values(by='original_row_id') line_coords[line] = group[["X", "Y"]].to_numpy() line_info[line] = group["original_row_id"].to_numpy() # 3. KD-Tree Search for Intersections and Mapping intersection_map = [] total_match_pairs = 0 print(f"Computing line intersections (Tolerance: {parallel_tolerance_deg} deg, Radius: {radius} m)...") for line1, line2 in combinations(line_ids, 2): A_coords = line_coords[line1] B_coords = line_coords[line2] A_ids = line_info[line1] # Original row_ids B_ids = line_info[line2] # Original row_ids if len(A_coords) < 2 or len(B_coords) < 2: continue # Parallel check runs on the entire collapsed line using endpoint angle if are_parallel(A_coords, B_coords, tolerance_degrees=parallel_tolerance_deg): continue tree = cKDTree(B_coords) matches = tree.query_ball_point(A_coords, r=radius) n_matches = sum(len(js) for js in matches) if n_matches == 0: continue total_match_pairs += n_matches intersection_key = '_'.join(sorted([line1, line2])) for i, js in enumerate(matches): if js: # Use original row_ids here intersection_map.append((A_ids[i], intersection_key, line1)) for j in js: intersection_map.append((B_ids[j], intersection_key, line2)) if not intersection_map: print("\nNo intersecting points found. Initializing flags to zero.") df_with_dummies[FLAG_COLUMN_INT] = 0 df_with_dummies[FLAG_SIG] = 0 df_with_dummies[FLAG_ALT] = 0 else: # 4. Create Intersection Map DataFrame and Calculate Flags df_map = pd.DataFrame(intersection_map, columns=['row_id', 'key', 'Line']).drop_duplicates() print(f"Marking {df_map['row_id'].nunique()} unique survey points as intersecting.") # Merge signal/alt values from the original df_with_dummies df_map = df_map.merge(df_with_dummies[['row_id', SIGNAL_COL, ALT_COL]], on='row_id', how='left') # Calculate means per line per intersection key line_means = df_map.groupby(['key', 'Line']).agg( mean_signal=(SIGNAL_COL, 'mean'), mean_alt=(ALT_COL, 'mean') ).reset_index() # --- MIN SIGNAL CALCULATION (ROBUST TIE-BREAKING) --- line_means["min_signal"] = line_means.groupby("key")["mean_signal"].transform("min") # FIX: Use np.isclose to handle floating point equality for ties (robust flagging) line_means[FLAG_SIG] = np.isclose( line_means["mean_signal"], line_means["min_signal"], atol=FLOAT_TOLERANCE ).astype(int) # --- MAX ALTITUDE CALCULATION (ROBUST TIE-BREAKING) --- line_means["max_alt"] = line_means.groupby("key")["mean_alt"].transform("max") # FIX: Use np.isclose to handle floating point equality for ties (robust flagging) line_means[FLAG_ALT] = np.isclose( line_means["mean_alt"], line_means["max_alt"], atol=FLOAT_TOLERANCE ).astype(int) # 5. Determine the final point-level flags based on line winners winning_sig_lines = line_means[line_means[FLAG_SIG] == 1][['key', 'Line']] winning_alt_lines = line_means[line_means[FLAG_ALT] == 1][['key', 'Line']] # Merge the winning line IDs back into df_map to retrieve ALL points df_sig_winners = df_map.merge(winning_sig_lines, on=['key', 'Line'], how='inner') sig_row_ids = df_sig_winners['row_id'].unique() df_alt_winners = df_map.merge(winning_alt_lines, on=['key', 'Line'], how='inner') alt_row_ids = df_alt_winners['row_id'].unique() # 6. Mark the original DataFrame (df_with_dummies) df_with_dummies[FLAG_COLUMN_INT] = 0 df_with_dummies[FLAG_SIG] = 0 df_with_dummies[FLAG_ALT] = 0 intersection_row_ids = df_map['row_id'].unique() df_with_dummies.loc[df_with_dummies['row_id'].isin(intersection_row_ids), FLAG_COLUMN_INT] = 1 df_with_dummies.loc[df_with_dummies['row_id'].isin(sig_row_ids), FLAG_SIG] = 1 df_with_dummies.loc[df_with_dummies['row_id'].isin(alt_row_ids), FLAG_ALT] = 1 # 7. Save Output base = os.path.splitext(os.path.basename(input_csv_path))[0] output_filename = f"{base}_MARKED_{int(radius * 100)}cm_GDB_FINAL.csv" outpath = os.path.join(output_path, output_filename) df_with_dummies.to_csv(outpath, index=False) t_end = time.time() print(f"✔ Processing complete. Saved marked CSV: {outpath}") print(f"Stage 2 Time: {t_end - t_start:.2f} seconds.") return outpath # -------------------------------------------------------- # STAGE 3: IMPORT CSV INTO GDB # -------------------------------------------------------- def import_csv_to_gdb(gdb_path, marked_csv_path, flag_channels): """Reads marked CSV and writes new flag channels back to the GDB.""" t_start = time.time() print(f"\n--- STAGE 3: IMPORTING CSV to GDB ---") # ... (Import logic as before) df_marked = pd.read_csv(marked_csv_path) line_flag_data = {line: group.set_index('row_id') for line, group in df_marked.groupby('Line')} try: with gxdb.Geosoft_gdb.open(gdb_path) as gdb_file: existing_chs = gdb_file.list_channels() for ch in flag_channels: if ch in existing_chs: gdb_file.delete_channel(ch) gdb_file.new_channel(ch, dtype=np.float32) gdb_lines = gdb_file.list_lines() x_ch = "X" for line in gdb_lines: gdb_file.current_line = line if line in line_flag_data: df_line = line_flag_data[line] try: gdb_point_count = len(flatten(gdb_file.read_channel(line, x_ch))) except Exception: gdb_point_count = 0 if gdb_point_count == 0: continue if len(df_line) != gdb_point_count: print( f"WARNING: Skipping line {line}. CSV length ({len(df_line)}) does not match GDB length ({gdb_point_count}).") continue for ch in flag_channels: flag_list = df_line[ch].astype(float).tolist() gdb_file.write_channel(line, ch, flag_list) print(f"✔ Flag channels successfully written to GDB: {gdb_path}") except Exception as e: print(f"FATAL GDB Import Error: {e}") print("Please check file permissions (GDB not open elsewhere) or channel types.") t_end = time.time() print(f"Stage 3 Time: {t_end - t_start:.2f} seconds.") # -------------------------------------------------------- # MAIN SCRIPT EXECUTION # -------------------------------------------------------- if __name__ == "__main__": t_total_start = time.time() os.makedirs(WORKING_DIR, exist_ok=True) # 1. EXPORT GDB marked_csv_input_path = export_gdb_to_csv(GDB_PATH, TEMP_CSV_INPUT, ESSENTIAL_CHANNELS) if marked_csv_input_path is None: print("\nProcess failed during GDB export. Cannot continue.") exit() # 2. PROCESS CSV try: marked_csv_output_path = mark_intersections_kdtree( input_csv_path=marked_csv_input_path, radius=RADIUS_M, parallel_tolerance_deg=PARALLEL_TOLERANCE_DEG, output_path=WORKING_DIR ) except Exception as e: print(f"\n--- FATAL ERROR IN STAGE 2 (Processing) ---: {e}") exit() # 3. IMPORT GDB try: import_csv_to_gdb(GDB_PATH, marked_csv_output_path, FLAG_CHANNELS) except Exception as e: print(f"\n--- FATAL ERROR IN STAGE 3 (GDB Import) ---: {e}") exit() # 4. Cleanup try: os.remove(TEMP_CSV_INPUT) print(f"\nCleanup: Removed temporary file {TEMP_CSV_INPUT}") except Exception as e: print(f"\nCleanup warning: Could not remove temporary file: {e}") t_total_end = time.time() print("\n--- ALL STAGES COMPLETE ---") print(f"Total Elapsed Time: {(t_total_end - t_total_start) / 60.0:.2f} minutes.")