RiverJoin_OpenInArcGIS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extract the FIM HF reaches which have corresponding SWORD reaches, and join them with the attributes of corresponding SWORD and IRIS reaches "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting the output paths and local variables [Modify based on your case. The (most possible) modifying lines have been indicated with a comment [Replace/Redefine...]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import arcpy\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "# Enable overwriting of output files\n",
    "arcpy.env.overwriteOutput = True\n",
    "\n",
    "# Set geodatabase path to store the output geo-data\n",
    "# The 1-4 paths for saving the intermediate reach and point feature classes, you can use separated paths or an overall one\n",
    "# The 5-6 paths for saving the all extracted FIM HF flowlines and those with joined IRIS attributes, respectively\n",
    "# You can rename the geodatabases (.gdb)\n",
    "# [Replace with your own paths]\n",
    "output_gdb_sword = r\"...\\Intermediate_SWORD.gdb\"\n",
    "output_gdb_point = r\"...\\Intermediate_Points.gdb\"\n",
    "output_gdb_reach = r\"...\\Intermediate_Reaches.gdb\"\n",
    "output_gdb_reach_point_perp = r\"...\\Intermediate_Reaches_Points_Perp.gdb\" # Path to save the intermediate outputs when using the perpendicular lines to extract the target flowlines\n",
    "output_gdb_final_reach = r\"...\\Extracted_Final_Reaches.gdb\"\n",
    "output_gdb_final_reach_iris = r\"...\\Extracted_Final_Reaches_Join_IRIS_SWORD.gdb\"\n",
    "\n",
    "# Set local variables\n",
    "arcpy.env.workspace = output_gdb_final_reach_iris"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: First extract FIM HF flowlines within SWORD flowline buffer [Normally no need to adjust]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_stream_network_within_sword_buffer(huc_id, original_fc, sword_fc, buffer_distance, unique_id):\n",
    "    # Step 1: Create Midpoints for Each Reach\n",
    "    midpoints_fc = f\"{output_gdb_point}/midpoints_{huc_id}\"\n",
    "    arcpy.management.FeatureVerticesToPoints(original_fc, midpoints_fc, \"MID\")\n",
    "    \n",
    "    # Step 2: \n",
    "    # Option1: Create Buffer Around SWORD Stream Network\n",
    "    # Extract numeric value and remove the ' Meters' part\n",
    "    buffer_distance_value = buffer_distance.split()[0]  # Gets 'buffer_distance'\n",
    "    sword_buffer_fc = f\"{output_gdb_sword}/sword_buffer_{buffer_distance_value}m\"\n",
    "    # print(sword_fc)\n",
    "    # print(sword_buffer_fc)\n",
    "    # print(buffer_distance)\n",
    "    arcpy.analysis.Buffer(sword_fc, sword_buffer_fc, buffer_distance)\n",
    "\n",
    "    # # Option2: Use the Existing Buffer Shapefile of SWORD Stream Network, then please comment the codes for Option 1\n",
    "    # sword_buffer_fc = \"Z:\\Documents\\RiverSlope\\Paper_NewRiverSlopeDataset\\Codes_Uploaded\\Demonstration_SpatialJoin\\Intermediate_SWORD.gdb/sword_buffer_120m\"\n",
    "   \n",
    "    # Step 3: Spatially Select Midpoints within SWORD Buffer\n",
    "    selected_midpoints_fc = f\"{output_gdb_point}/selected_midpoints_{huc_id}\"\n",
    "    arcpy.analysis.SpatialJoin(midpoints_fc, sword_buffer_fc, selected_midpoints_fc, join_type=\"KEEP_COMMON\", match_option=\"WITHIN\")\n",
    "    \n",
    "    # Step 4: Select and Extract Matching Reaches Using ORIG_FID\n",
    "    # Create a set of unique ORIG_FID values from the selected midpoints\n",
    "    orig_fid_list = [row[0] for row in arcpy.da.SearchCursor(selected_midpoints_fc, unique_id)]\n",
    "    # print(orig_fid_list, \"\\n\")\n",
    "\n",
    "    for field in arcpy.ListFields(original_fc):\n",
    "        if field.name == unique_id:\n",
    "            if field.type == \"String\":\n",
    "                # Build a where clause to select the matching reaches\n",
    "                where_clause = \"{} IN ({})\".format(unique_id, ','.join([f\"'{fid}'\" for fid in orig_fid_list]))\n",
    "                # print((f\"WHERE CLAUSE (string): {where_clause}\"))\n",
    "            else:\n",
    "                where_clause = f\"{unique_id} IN ({','.join(map(str, orig_fid_list))})\"\n",
    "                # print((f\"WHERE CLAUSE (other): {where_clause}\"))\n",
    "    \n",
    "    # Select the original stream reaches that correspond to the selected midpoints\n",
    "    arcpy.management.MakeFeatureLayer(original_fc, f\"original_lyr_{huc_id}\", where_clause)\n",
    "    # print(\"Select the original stream reaches successed!\", \"\\n\")\n",
    "    \n",
    "    # Copy the selected reaches to the output feature class\n",
    "    output_fc = f\"{output_gdb_reach}/extracted_reaches_{huc_id}\"\n",
    "    arcpy.management.CopyFeatures(f\"original_lyr_{huc_id}\", output_fc)\n",
    "    \n",
    "    print(f\"Step 1: First extracted flowlines saved to {output_fc}\")\n",
    "    \n",
    "    return output_fc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Trace downstream for each first extracted flowline to get the not initially extracted flowlines [Normally no need to adjust]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reconstruct_disconnected_segments(extract_stream_fc, unique_ids_current_nextDown, unique_id, unique_id_nextDown, original_fc, huc_id):\n",
    "    disconnected_segments = []\n",
    "    \n",
    "    # Initialize lists to collect the HydroID and NextDownID values\n",
    "    hydro_ids = []\n",
    "    nextdown_ids = []\n",
    "    \n",
    "    # Step 1: Identify ends of disconnected segments\n",
    "    with arcpy.da.SearchCursor(extract_stream_fc, unique_ids_current_nextDown) as cursor:\n",
    "        for row in cursor:\n",
    "            hydro_ids.append(row[0])     # Collect HydroID values\n",
    "            nextdown_ids.append(row[1])  # Collect NextDownID values\n",
    "\n",
    "    # Convert the lists to NumPy arrays\n",
    "    hydro_ids_np = np.array(hydro_ids)\n",
    "    nextdown_ids_np = np.array(nextdown_ids)\n",
    "    \n",
    "    disconnected_ends = np.unique(nextdown_ids_np[~np.isin(nextdown_ids_np, hydro_ids_np)])  # IDs that should be upstream but are missing\n",
    "    # print(f\"disconnected_ends: {disconnected_ends}\")\n",
    "\n",
    "    # Step 2: Trace upstream and find missing reaches\n",
    "    for end_id in disconnected_ends:\n",
    "        current_id = end_id\n",
    "        missing_reaches = []\n",
    "\n",
    "        while current_id is not None:\n",
    "            # Check if this reach exists in the original dataset\n",
    "\n",
    "            for field in arcpy.ListFields(original_fc):\n",
    "                if field.name == unique_id:\n",
    "                    if field.type == \"String\":\n",
    "                        where_clause = f\"{unique_id} = '{current_id}'\"\n",
    "                        # print(f\"where_clause (string): {where_clause}\")\n",
    "                    else:\n",
    "                        where_clause = f\"{unique_id} = {current_id}\"\n",
    "                        # print(f\"where_clause (other): {where_clause}\")\n",
    "                    unique_id_field_type = field.type\n",
    "            \n",
    "            found_row = False  # Flag to detect if any rows were found\n",
    "            \n",
    "            with arcpy.da.SearchCursor(original_fc, [unique_id, unique_id_nextDown], where_clause) as cursor:\n",
    "                for row in cursor:\n",
    "                    # print(row, \"\\n\")\n",
    "                    found_row = True                    \n",
    "                    if current_id not in hydro_ids: \n",
    "                        missing_reaches.append(current_id)  # Add to missing reaches if it's not in the extracted network\n",
    "#                         print(current_id)\n",
    "                   \n",
    "                    current_id = row[1]  # Move to the next downstream reach\n",
    "#                     print(current_id, \"\\n\")\n",
    "                    \n",
    "            if not found_row:\n",
    "#                 print(f\"No matching record found for HydroID {current_id}. Exiting loop.\")\n",
    "                break  # Exit the while loop if no matching rows were found\n",
    "\n",
    "            if current_id in hydro_ids or current_id is None:\n",
    "                break  # Stop if we've reached a connected reach or the end of the network\n",
    "\n",
    "        if missing_reaches:\n",
    "            disconnected_segments.append(missing_reaches)\n",
    "            \n",
    "    # Step 3: Add missing reaches to the extracted network\n",
    "    missing_reaches_fc = f\"{output_gdb_reach}/missing_reaches_{huc_id}\"\n",
    "    arcpy.management.MakeFeatureLayer(original_fc, f\"original_lyr_{huc_id}\")\n",
    "    \n",
    "    for segment in disconnected_segments:\n",
    "        # print(f\"segment: {segment}\")\n",
    "        \n",
    "        if unique_id_field_type == \"String\":\n",
    "            where_clause = \"{} IN ({})\".format(unique_id, ','.join([f\"'{fid}'\" for fid in segment]))\n",
    "        else:\n",
    "            where_clause = f\"{unique_id} IN ({','.join(map(str, segment))})\"\n",
    "                        \n",
    "        arcpy.management.SelectLayerByAttribute(f\"original_lyr_{huc_id}\", \"ADD_TO_SELECTION\", where_clause)\n",
    "    \n",
    "    arcpy.management.CopyFeatures(f\"original_lyr_{huc_id}\", missing_reaches_fc)\n",
    "    \n",
    "    # Merge the missing reaches back into the extracted stream network\n",
    "    final_fc = f\"{output_gdb_reach}/traced_extracted_reaches_{huc_id}\" # Define the path to save the final reconstructed network\n",
    "    arcpy.management.Merge([extract_stream_fc, missing_reaches_fc], final_fc)\n",
    "    \n",
    "    print(f\"Step 2: Traced and first extracted flowlines saved to {final_fc}\")\n",
    "\n",
    "    return final_fc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Check if there are unnormal situations [Modify based on your case]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Situation 1: If count of extracted_reaches == 0 --> no SWORD flowline intersects with FIM HF flowlines in this HUC --> do not take any FIM HF flowlines\n",
    "### Situation 2: If count of traced_extracted_reaches = count of extracted_reaches + count of all original FIM HF flowlines in this HUC --> the extracted_reaches already have all the FIM HF flowlines corresponding to SWORD flowlines in this HUC --> just use the extracted_reaches (this is due to the deficiency of our downstream tracing algorithm - Step 2)\n",
    "### Situation 3: If not the previous two cases, suggesting there still are possibly upstream reaches and others to be extracted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_flowlines(extract_stream_fc, traced_extracted_reaches, original_fc, huc_id):\n",
    "    # Get count of features in each shapefile\n",
    "    count_extract = int(arcpy.management.GetCount(extract_stream_fc)[0])\n",
    "    count_trace_extract = int(arcpy.management.GetCount(traced_extracted_reaches)[0])\n",
    "    count_original = int(arcpy.management.GetCount(original_fc)[0])\n",
    "\n",
    "    # print(f\"   Extracted reaches count: {count_extract}\")\n",
    "    # print(f\"   Traced and extracted reaches count: {count_trace_extract}\")\n",
    "    # print(f\"   Original reaches count: {count_original}\")\n",
    "\n",
    "    # Situation 1: If extract_stream_fc count is 0, skip further processing\n",
    "    if count_extract == 0:\n",
    "        situation_identify = \"no_intersection\"\n",
    "        print(\"Step 3: There may be no intersection of FIM HF and SWORD flowlines in this HUC.\")\n",
    "\n",
    "    # Situation 2: If traced_extracted_reaches count = extract_stream_fc + original_fc, copy the extracted streams\n",
    "    elif count_trace_extract == count_extract + count_original:\n",
    "        output_fc = f\"{output_gdb_final_reach}\\\\all_extracted_reaches_{huc_id}\"\n",
    "        arcpy.management.CopyFeatures(extract_stream_fc, output_fc)\n",
    "\n",
    "        # Enable automatic addition of outputs to the map\n",
    "        arcpy.env.addOutputsToMap = True\n",
    "\n",
    "        # Define fields for duplicate detection\n",
    "        # [Replace based on your case]\n",
    "        fields = [\"HydroID\", \"From_Node\", \"To_Node\"]\n",
    "        arcpy.management.DeleteIdentical(output_fc, fields)\n",
    "\n",
    "        situation_identify = \"all_extracted\"\n",
    "        \n",
    "        print(f\"Step 3: All corresponding flowlines extracted have been removed duplicate records and saved to {output_fc}\")\n",
    "        \n",
    "    else:\n",
    "        # Situation 3: Otherwise, print that more flowlines need to be extracted\n",
    "        situation_identify = \"extract_upstream_others\"\n",
    "\n",
    "        print(\"Step 3: Maybe the most upstream flowlines (and other) still need to be extracted using Step 4.\")\n",
    "\n",
    "    return situation_identify\n",
    "    # print(f\"Identified situation: {situation_identify}\")\n",
    "\n",
    "    # Disable automatic addition of outputs to the map\n",
    "    arcpy.env.addOutputsToMap = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Extract the (possibly) upstream (and other) flowlines that have not yet been extracted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4.1: Remove the perpendicular lines intersected with traced_extracted_reaches [Normally no need to adjust]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def remove_touching_perpendiculars(traced_extracted_reaches, perpendicular_lines, huc_id):\n",
    "    # Make a copy of perpendicular_lines for modification\n",
    "    output_fc = f\"{output_gdb_reach_point_perp}\\\\perp_lines_notcross_{huc_id}\"\n",
    "    arcpy.CopyFeatures_management(perpendicular_lines, output_fc)\n",
    "\n",
    "    # Create a feature layer from the copied perpendicular lines\n",
    "    arcpy.MakeFeatureLayer_management(output_fc, \"perpendicular_lines_layer\")\n",
    "\n",
    "    # Select perpendicular lines that TOUCH the flowlines in traced_extracted_reaches\n",
    "    arcpy.management.SelectLayerByLocation(\"perpendicular_lines_layer\", \n",
    "                                           overlap_type=\"INTERSECT\", \n",
    "                                           select_features=traced_extracted_reaches, \n",
    "                                           selection_type=\"NEW_SELECTION\")\n",
    "\n",
    "    # Count the selected features\n",
    "    count = int(arcpy.management.GetCount(\"perpendicular_lines_layer\")[0])\n",
    "    \n",
    "    if count > 0:\n",
    "        # Delete the selected features that touch flowlines\n",
    "        arcpy.management.DeleteFeatures(\"perpendicular_lines_layer\")\n",
    "        # print(f\"Removed {count} perpendicular lines that touched traced flowlines.\")\n",
    "\n",
    "    # print(f\"     Step 4.1: Saved perpendicular lines not touching traced flowlines to {output_fc}\")\n",
    "\n",
    "    return output_fc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4.2: Do the intersect join of all original FIM HF flowlins of this HUC and remaining perpendicular lines [Normally no need to adjust]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to perform the spatial join for a given HUC ID\n",
    "def perform_spatial_join_with_perp(huc_id, original_fc, perp_lines_notcross):\n",
    "    \n",
    "    output_fc = os.path.join(output_gdb_reach_point_perp, f\"original_reaches_{huc_id}_join_perp\")\n",
    "    \n",
    "    # Perform the spatial join\n",
    "    arcpy.analysis.SpatialJoin(\n",
    "        target_features = original_fc,\n",
    "        join_features = perp_lines_notcross,\n",
    "        out_feature_class = output_fc,\n",
    "        join_type = \"KEEP_ALL\",  # Keeps all features\n",
    "        match_option = \"INTERSECT\",  # Match on intersect\n",
    "        join_operation = \"JOIN_ONE_TO_ONE\"\n",
    "    )\n",
    "    \n",
    "    # print(f\"     Step 4.2: Spatial join output saved at {output_fc}\")\n",
    "    \n",
    "    return output_fc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4.3: Set filter to find only the really needed reaches (e.g. most upstream, the main stream reaches corresponding to SWORD reaches) [Modify based on your case]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def filter_intersected_reaches(huc_id, ori_reaches_join_perp):\n",
    "\n",
    "    # Set paths\n",
    "    output_fc = f\"{output_gdb_reach_point_perp}/original_reaches_{huc_id}_join_perp_Seled_count3\"  # Replace with the desired output path for selected features\n",
    "\n",
    "    # Create a feature layer from the spatial join output\n",
    "    arcpy.management.MakeFeatureLayer(ori_reaches_join_perp, f\"joined_layer_{huc_id}\")\n",
    "\n",
    "    # Apply the selection based on Join_Count and order_ conditions\n",
    "    # Select Join_Count >= 3 OR (Join_Count >= 1 AND order_ > 2)\n",
    "    # [Redefine the filter rules based on your case, after having some local inspections]\n",
    "    where_clause = \"\"\"(\"Join_Count\" >= 3) OR (\"Join_Count\" >= 1 AND \"order_\" > 2)\"\"\"\n",
    "#     where_clause = \"\"\"(\"Join_Count\" >= 3)\"\"\" # can reduce the wrongly got reaches\n",
    "    arcpy.management.SelectLayerByAttribute(f\"joined_layer_{huc_id}\", \"NEW_SELECTION\", where_clause)\n",
    "\n",
    "    # Save the selected features to a new feature class\n",
    "    arcpy.management.CopyFeatures(f\"joined_layer_{huc_id}\", output_fc)\n",
    "    \n",
    "    # print(f\"     Step 4.3: Selected features saved to {output_fc}\")\n",
    "    \n",
    "    return output_fc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4.4: Combine all the extracted reaches including traced_extracted_reaches (Steps 1-2), join_perp_Seled (Steps 4.1-4.3). [Modify based on your case]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def merge_and_remove_duplicates(traced_extracted_reaches, huc_id):\n",
    "    # Define output feature class\n",
    "    all_extracted_reaches = f\"{output_gdb_final_reach}\\\\all_extracted_reaches_{huc_id}\"\n",
    "\n",
    "    # Merge the two feature classes\n",
    "    arcpy.management.Merge([join_perp_Seled_reaches, traced_extracted_reaches], all_extracted_reaches)\n",
    "    print(f\"Step 4: Merged all features saved to {all_extracted_reaches}\")\n",
    "\n",
    "    # Disable automatic addition of outputs to the map\n",
    "    arcpy.env.addOutputsToMap = True\n",
    "\n",
    "    # Define fields for duplicate detection\n",
    "    # [Replace based on your case]\n",
    "    fields = [\"HydroID\", \"From_Node\", \"To_Node\"]\n",
    "\n",
    "    # Use Delete Identical tool to remove duplicates based on the specified fields\n",
    "    arcpy.management.DeleteIdentical(all_extracted_reaches, fields)\n",
    "    # print(f\"Removed duplicate records based on {fields}\")\n",
    "\n",
    "    # Get the field names for both feature classes\n",
    "    all_extracted_fields = {field.name for field in arcpy.ListFields(all_extracted_reaches)}\n",
    "    original_fields = {field.name for field in arcpy.ListFields(original_fc)}\n",
    "\n",
    "    # Identify fields to remove (present in all_extracted_fields but not in original_fields)\n",
    "    # [Replace the specific not-in fields based on your case, if needed]\n",
    "    fields_to_remove = [field for field in all_extracted_fields if field not in original_fields and field not in ['Shape_Length', 'Shape', 'geom_Length', 'OBJECTID']]\n",
    "    \n",
    "    # # Print fields that will be removed\n",
    "    # print(\"Fields to be removed:\", fields_to_remove)\n",
    "\n",
    "    if fields_to_remove:\n",
    "        # Delete extra fields\n",
    "        arcpy.management.DeleteField(all_extracted_reaches, fields_to_remove)\n",
    "        # print(f\"Removed fields: {fields_to_remove}\")\n",
    "\n",
    "    return all_extracted_reaches"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5: Join reach attributes of IRIS, SWORD to all extracted reaches [Modify based on your case]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [You may still need to replace some specific field names based on your case in the following lines]\n",
    "\n",
    "def join_sword_iris_to_reaches(huc_id, all_extract_reaches, sword_fc, iris_fc, iris_reach_id_field):\n",
    "\n",
    "    # Disable automatic addition of outputs to the map\n",
    "    arcpy.env.addOutputsToMap = False\n",
    "\n",
    "    # Step 5.1: Generate Midpoints for All Extracted Reaches\n",
    "    all_reach_midpoints_fc = f\"{output_gdb_point}/all_extracted_reaches_midp_{huc_id}\"\n",
    "    arcpy.management.FeatureVerticesToPoints(all_extract_reaches, all_reach_midpoints_fc, \"MID\")\n",
    "\n",
    "    # Step 5.2: Find the Closest SWORD Reach for Each Midpoint\n",
    "    arcpy.analysis.Near(all_reach_midpoints_fc, sword_fc, location=\"NO_LOCATION\", angle=\"NO_ANGLE\", method=\"PLANAR\")\n",
    "    \n",
    "    # Step 5.2.1: Add a new field to the target feature class to store NEAR_FID first\n",
    "    arcpy.management.AddField(all_extract_reaches, \"NEAR_FID\", \"BIGINTEGER\")\n",
    "\n",
    "    # Step 5.2.2: Join Midpoints with Reaches Based on OBJECTID\n",
    "    # Join midpoints to reaches based on HydroID\n",
    "    joined_layer = arcpy.management.AddJoin(all_extract_reaches, \"OBJECTID\", all_reach_midpoints_fc, \"ORIG_FID\")\n",
    "    \n",
    "#     fields = [f.name for f in arcpy.ListFields(joined_layer)]\n",
    "#     print(fields)\n",
    "\n",
    "    # Copy the NEAR_FID from midpoints to the original reaches\n",
    "    arcpy.management.CalculateField(joined_layer, \"NEAR_FID\", f\"!all_extracted_reaches_midp_{huc_id}.NEAR_FID!\", \"PYTHON3\")\n",
    "\n",
    "    # Remove the join after copying NEAR_FID\n",
    "    arcpy.management.RemoveJoin(joined_layer, f\"all_extracted_reaches_midp_{huc_id}\")\n",
    "\n",
    "    # Step 5.3: Join SWORD Attributes to the Extracted Reaches Using NEAR_FID\n",
    "    sword_oid_field = \"OBJECTID\"  # Assuming SWORD's unique ID field is 'OBJECTID'\n",
    "    arcpy.management.AddJoin(joined_layer, \"NEAR_FID\", sword_fc, sword_oid_field)\n",
    "    \n",
    "#     fields = [f.name for f in arcpy.ListFields(joined_layer)]\n",
    "#     print(fields)\n",
    "    \n",
    "#     fields = [f.name for f in arcpy.ListFields(iris_fc)]\n",
    "#     print(fields)\n",
    "\n",
    "    # Step 5.4: Join IRIS Layer to Extracted Reaches\n",
    "    arcpy.management.AddJoin(joined_layer, \"main_reaches.reach_id\", iris_fc, iris_reach_id_field)   # [Replace with that in your data]\n",
    "    \n",
    "#     fields = [f.name for f in arcpy.ListFields(joined_layer)]\n",
    "#     print(fields)\n",
    "\n",
    "    # Enable automatic addition of outputs to the map\n",
    "    arcpy.env.addOutputsToMap = True\n",
    "\n",
    "    # Step 5.5: Export the joined layer to GDB\n",
    "    # Define the output shapefile path\n",
    "    output_filepath = os.path.join(output_gdb_final_reach_iris, f\"all_extracted_reaches_iris_{huc_id}\")\n",
    "\n",
    "    # Export the joined layer to a shapefile\n",
    "    arcpy.conversion.ExportFeatures(joined_layer, output_filepath)\n",
    "    print(f\"Step 5: All extracted reaches joined with IRIS saved to {output_filepath}\")\n",
    "\n",
    "    # Disable automatic addition of outputs to the map\n",
    "    arcpy.env.addOutputsToMap = False\n",
    "\n",
    "    # Enable overwriting of output files\n",
    "    arcpy.env.overwriteOutput = True\n",
    "    \n",
    "    # Define the output CSV file path\n",
    "    output_csv = os.path.join(output_gdb_final_reach_iris, f\"csv_all_extracted_reaches_iris_{huc_id}\") \n",
    "\n",
    "    # Export the attribute table to a CSV file\n",
    "    arcpy.conversion.ExportTable(output_filepath, output_csv)\n",
    "    # print(f\"Step 5: csv of all extracted reaches joined with IRIS saved to {output_csv}\")\n",
    "\n",
    "    # return(output_filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running all the functions [Modify based on your case]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Step 1: First extracted flowlines saved to Z:\\Documents\\RiverSlope\\Paper_NewRiverSlopeDataset\\Codes_Uploaded\\Demonstration_SpatialJoin\\Intermediate_Reaches.gdb/extracted_reaches_01020004\n",
      "Step 2: Traced and first extracted flowlines saved to Z:\\Documents\\RiverSlope\\Paper_NewRiverSlopeDataset\\Codes_Uploaded\\Demonstration_SpatialJoin\\Intermediate_Reaches.gdb/traced_extracted_reaches_01020004\n",
      "Step 3: Maybe the most upstream flowlines (and other) still need to be extracted using Step 4.\n",
      "Step 4: Merged all features saved to Z:\\Documents\\RiverSlope\\Paper_NewRiverSlopeDataset\\Codes_Uploaded\\Demonstration_SpatialJoin\\Extracted_Final_Reaches.gdb\\all_extracted_reaches_01020004\n",
      "Step 5: All extracted reaches joined with IRIS saved to Z:\\Documents\\RiverSlope\\Paper_NewRiverSlopeDataset\\Codes_Uploaded\\Demonstration_SpatialJoin\\Extracted_Final_Reaches_Join_IRIS_SWORD.gdb\\all_extracted_reaches_iris_01020004\n"
     ]
    }
   ],
   "source": [
    "# Disable automatic addition of outputs to the map\n",
    "arcpy.env.addOutputsToMap = False\n",
    "\n",
    "## Step 1: First extracted reaches using SWORD buffer\n",
    "huc_id = \"03020201\"   # Unique HUC8 id [Global unique identifier, especially useful when you iteratively process multiple datasets; replace with your own (data) id] \n",
    "# HUC 03020201 FIM hydrofabric flowline (reach) path [Replace with your path]\n",
    "original_fc = rf\"...\\InputGeoData.gdb/main_demDerived_reaches_split_filtered_addedAttributes_crosswalked_0\"\n",
    "unique_id = \"HydroID\"   # Unique identifier of FIM flowline [Replace with that in your data]\n",
    "\n",
    "# SWORD flowlines path, i.e. the flowlines your want to join to your target flowlines (FIM HF flowlines here) [Replace with yours]\n",
    "sword_fc = r\"...\\InputGeoData.gdb/main_reaches\"\n",
    "buffer_distance = \"100 Meters\"   # [Replace with the one you want specify]\n",
    "\n",
    "extract_stream_fc = extract_stream_network_within_sword_buffer(huc_id, original_fc, sword_fc, buffer_distance, unique_id)\n",
    "\n",
    "## Step 2: Trace downstream to find the initially not extracted\n",
    "unique_id_nextDown = \"NextDownID\"  # Next downstream FIM flowline unique identifier [Replace with that in your data]\n",
    "unique_ids_current_nextDown = [unique_id, unique_id_nextDown] # Unique identifier of FIM flowline and its downstream\n",
    "traced_extracted_reaches = reconstruct_disconnected_segments(extract_stream_fc, unique_ids_current_nextDown, unique_id, unique_id_nextDown, original_fc, huc_id)\n",
    "\n",
    "## Step 3: Deal with the unnormal situations\n",
    "situation = process_flowlines(extract_stream_fc, traced_extracted_reaches, original_fc, huc_id)\n",
    "\n",
    "## Step 4: Find upstream and other reaches using perpendicular lines of SWORD reaches\n",
    "if situation == \"no_intersection\":\n",
    "    \n",
    "    pass\n",
    "\n",
    "elif situation == \"all_extracted\" or situation == \"extract_upstream_others\":\n",
    "    \n",
    "    if situation == \"extract_upstream_others\":\n",
    "\n",
    "        # Step 4.1: Remove the perp lines intersected with already extracted FIM flowlines\n",
    "        # Due to it's difficult for us to generate perpendicular lines of SWORD reaches in ArcGIS, so we generated them separately in Python, then used the resulting feature class here\n",
    "        # Python codes to generate the perp lines has been provided as 'Step4.0_PerpendicularLines_Generation.py' in this GitHub\n",
    "        perpendicular_lines = r'...\\InputGeoData.gdb\\perpendicular_lines_len0006_space_sword_nodes_cleaned'  # [Replace with your actual path]\n",
    "        perp_lines_notcross = remove_touching_perpendiculars(traced_extracted_reaches, perpendicular_lines, huc_id)\n",
    "    \n",
    "        # Step 4.2: Join the remaining perp lines with original FIM flowlines\n",
    "        ori_reaches_join_perp = perform_spatial_join_with_perp(huc_id, original_fc, perp_lines_notcross)\n",
    "    \n",
    "        # Step 4.3: Filter the joined FIM flowlines to get the wanted ones\n",
    "        join_perp_Seled_reaches = filter_intersected_reaches(huc_id, ori_reaches_join_perp)\n",
    "    \n",
    "        # Step 4.4: Merge all the flowlines extracted from different steps\n",
    "        merge_and_remove_duplicates(traced_extracted_reaches, huc_id)\n",
    "\n",
    "    ## Step 5: Join IRIS and SWORD to all extracted reaches\n",
    "    all_extract_reaches = f\"{output_gdb_final_reach}\\\\all_extracted_reaches_{huc_id}\"   # [Modify this if necessary]\n",
    "    iris_fc = r\"...\\InputGeoData.gdb/IRIS_v26\"   # The flwolines you want to link via the common flowlines, i.e. IRIS here [Replace with yours if necessary]\n",
    "    iris_reach_id_field = \"reach_id\"  # The common field between SWORD and IRIS [Replace with that in your data, if necessary]\n",
    "    join_sword_iris_to_reaches(huc_id, all_extract_reaches, sword_fc, iris_fc, iris_reach_id_field)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ArcGISPro",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "version": "3.11.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}