Update Features and Analysis in Arcpy

Arcpy generated Point Feature Data from Ground water contamination, PFAS chemical levels and MNPCA clean up activity sites

Multiple Shapefile data sources access programmatically online and converted to points in Arcpy, then mapped with ArcGIS pro

What is Arcpy?:¶

Arcpy is a Python library for automating GIS tasks in Arc GIS Pro

Collection of modules, functions, and classes from the ArcGIS toolbox
It allows you to perform geoprocessing, analysis, data management and mapping automation using python

In [1]:

import arcpy # allows for access to ArcGIS pro geoprocessing tools and workflow automation
import os # enables interaction with local system resources (paths to folders and files)
import requests # access data from the web
import zipfile ## process and extract downloaded zipfiles 
import pandas as pd ## use to check spreadsheet formatting


initial_dir= os.getcwd() # capture the initial directory before setting environment for arcpy
print("Initial working directory: ", initial_dir)

Initial working directory:  c:\Projects\my_git_pages_website\Py-and-Sky-Labs\content\ArcPY

What we did last time¶

Created a custom function to pull a countytownship shapefile and load dataset into our local geodatabase
Perform several QC checks on attributes and geometry
Created a unique key by joining attributes
Performed Spatial Join within Arcpy

What we will do this time:¶

Update Features using Arcpy
Perform Spatial Analysis within Arcpy taking each water well's latest reading per polygon

Create and re-use Custom Functions to Save time¶

When you create a python script, the functions within those scripts can be loaded into other scripts
- To import your script, do as you would any library in python
  - python import my_script as ms
- Now the functions within this script are available to use
  - python ms.my_function()
Not used here but when multiple script hold different tools you'd be better served using a package and __init__.py setup

In [2]:

## Let's load in our custom tool
import custom_arcpy_tools as cat

Initial working directory:  c:\Projects\my_git_pages_website\Py-and-Sky-Labs\content\ArcPY

In [ ]:

# note: importing the custom tool sets a variable called initial_dir to the path of the script
        # this is intentionally done to support its functionality i.e., creating folders to store data

In [3]:

cat.listFC_dataset_Attributes()

## We purposefully caused an error!
# Why? To demonstrate that the error is due to workspace not being set in this new notebook or script

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_18168\2419765190.py in <cell line: 0>()
----> 1 cat.listFC_dataset_Attributes()
      2 
      3 ## We purposefully caused an error!
      4 # Why? To demonstrate that the error is due to workspace not being set in this new notebook or script

c:\Projects\my_git_pages_website\Py-and-Sky-Labs\content\ArcPY\custom_arcpy_tools.py in listFC_dataset_Attributes(workspace)
     92     if not arcpy.env.workspace:
     93 
---> 94         raise ValueError("No workspace is set! Please specify a workspace") #stops the program, calls out this error
     95 
     96 

ValueError: No workspace is set! Please specify a workspace

In [ ]:

# See the ValueError? 
### Don't remember the arguments for this function? 
### Use this:
#  help(function_name)

In [ ]:

help(cat.listFC_dataset_Attributes) ## use help(function_name) to get the documentation for the function

Help on function listFC_dataset_Attributes in module custom_arcpy_tools:

listFC_dataset_Attributes(workspace=None)
    List feature classes in the workspace and checks their key attributes, such as SpatialReference , 
    geometry type, and coordinate systems
    Parameters:
    -----------
    workspace: str, optional
        The file path to the workspace (e.g., geodatabase path) containing the feature classes
        If none, the function will use the current ArcPy workspace
        If no workspace is set will raise an error telling you to set a a workspace
    
    Returns:
    --------
    None , if no workspace exists
    Otherwise,    Prints the details of all feature classes in the wksp , including
                - Feature Class Name
                - spatial reference name
                - spatial reference type
                - Geometry
                - Well-known ID of the spatial reference
                Also, warns if different spatial references or coordinate systems exists in the dataset
    Raises: 
    -------
    Value Error: 
        If no wksp is set or if not feature classes are found in the workspace
    Exception:
        For general errors that may occur during processing
    
    Example:
    # To list fc in a geodatabase
    >>> listFC_dataset_Attributes(r"c:/Path/To/Your/Workspace.gdb")
    
    # To use the current wksp already set in Arcpy
    >>> listFC_dataset_Attributes()

In [ ]:

# Let's change the directory to where we have a geodatabase with some feature classes to work with
# Note: % is a magic command to access local system file and folders
%cd ./02 Result

c:\Projects\my_git_pages_website\Py-and-Sky-Labs\content\ArcPY\02 Result

In [4]:

# list the files
%ls

 Volume in drive C is OS
 Volume Serial Number is DA19-2472

 Directory of c:\Projects\my_git_pages_website\Py-and-Sky-Labs\content\ArcPY

03/28/2025  03:15 PM    <DIR>          .
04/02/2025  09:38 AM    <DIR>          ..
04/02/2025  02:34 PM    <DIR>          __pycache__
03/17/2025  03:46 PM    <DIR>          01 Data
03/21/2025  10:29 AM    <DIR>          02 Result
03/11/2025  12:38 AM    <DIR>          03 Map
03/06/2025  11:44 PM            39,845 ArcGIS_Documentation.png
04/02/2025  10:20 AM           486,810 Arcpy_Spatial_Join_part4.ipynb
04/02/2025  03:27 PM           693,787 Arcpy_UpdateFeatures_JoinAnalysis_part5.ipynb
03/11/2025  09:50 AM                84 category.md
04/02/2025  09:48 AM            23,301 custom_arcpy_tools.py
03/17/2025  10:12 AM             1,318 MN--Ground Water Contamination data.txt
03/28/2025  03:15 PM           488,824 SpatialJoinResult_JoinCount0.png
03/17/2025  10:08 PM            50,674 what_is_arcpy.ipynb
03/21/2025  10:12 PM            62,557 what_is_arcpy_part2.ipynb
03/26/2025  11:38 PM            55,590 what_is_arcpy_part3.ipynb
              10 File(s)      1,902,790 bytes
               6 Dir(s)  434,748,108,800 bytes free

In [ ]:

## should see a file ending with .gdb. Copy the folder path and add its name below

In [5]:

## Use the gdb name in the folder, and set the path to the gdb

###======================    Set up file paths


gdb_path = r"c:\Projects\my_git_pages_website\Py-and-Sky-Labs\content\ArcPY\02 Result\MN_WaterContamination.gdb"

##===================== Define the path to script and project folder

# set the folder path as the working environment, and store it for later use
try:
    arcpy.env.workspace = gdb_path
    wksp = arcpy.env.workspace
except Exception as e:
    print(f"Error setting up the path, {e}, Exception type: {type(e).__name__}") # python's built-in error messaging, {e} prints description and type(e).__name__ category

print("Working environment is here",wksp)

## Modify the env attribute to allow overwriting output files
arcpy.env.overwriteOutput= True

Working environment is here c:\Projects\my_git_pages_website\Py-and-Sky-Labs\content\ArcPY\02 Result\MN_WaterContamination.gdb

In [6]:

arcpy.Exists(wksp) ## similar to using print statement following if not os.path.exists(gdb_path)

Out[6]:

True

In [7]:

# Now. Retry the function that failed earlier

cat.listFC_dataset_Attributes()

Workspace is set here c:\Projects\my_git_pages_website\Py-and-Sky-Labs\content\ArcPY\02 Result\MN_WaterContamination.gdb

Feature class: , atlas_contamination_areas, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPolygon, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , atlas_contamination_area_lines, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPolyline, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , atlas_contamination_sites, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPoint, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , atlas_contamination_sites_noca, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPoint, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , atlas_contamination_source_areas, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPolygon, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , atlas_contamination_wells_result_sum, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPoint, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , MN_PFAS_Levels, Spatial Reference name :, GCS_North_American_1983, Spatial Ref Type : Geographic, GeometryPoint, WKID: 4269
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , MN_PFAS_Levels_proj, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPoint, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , MN_PFAS_Levels_by_City, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPolygon, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , city_township_unorg, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPolygon, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , MN_WaterContamination_SpatialJoin, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPoint, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , MN_latest_WaterSamples_Joined_lyr, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPoint, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , overlap_check1, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPolygon, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , overlap_check2, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPolygon, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 


14 Feature classes found in the workspace path

Warning Different wkid detected among feature Classes:  {26915, 4269}
Warning Different Coordinate Systems types found among Feature Classes:  {'Projected', 'Geographic'}

Data Sources:¶

Geography: MN

Input Name	file types	Coordinate System	Source
MN PFAS levels	1 excel file	GCS: NAD 83	Report
MN Groundwater contamination atlas	6 shape files	PCS: NAD 83 UTM Zone 15N	MetaData
MN Cities and townships	1 shapefile	PCS: NAD 83 UTM Zone 15N	MetaData

We have already pulled and done some initial processing of all these datasets in the previous walkthroughs
CTU poygon dataset was spatially joined to the MN wells dataset, bucketing wells by CTU
Now will go into analysis of joined data

In [8]:

# quick reminder of where our current workspace is set
print("Data Folder is here" , wksp)
# let's call it gdb for clarity
gdb=wksp

Data Folder is here c:\Projects\my_git_pages_website\Py-and-Sky-Labs\content\ArcPY\02 Result\MN_WaterContamination.gdb

In [9]:

## Count initial records
# first check the join fc: County/Township/Unorganized territory (CTU)

fc= "city_township_unorg"

print(f"Total records in {fc} feature class is:  {int(arcpy.management.GetCount(fc)[0])}")

Total records in city_township_unorg feature class is:  2744

In [10]:

## check the target input fc

for fc in arcpy.ListFeatureClasses("*wells_result_sum*"):
    print(f"Feature class: {fc}\n Total Record Count: {int(arcpy.management.GetCount(fc)[0])}")

Feature class: atlas_contamination_wells_result_sum
 Total Record Count: 17443

In [11]:

fc ="city_township_unorg"
cat.showFieldinfo(fc)

*Summary*
 Fields in the Polygon Feature Class: city_township_unorg

Name                 Type         Length   Unique_Values
 OBJECTID             | OID        | 4     | 2744      
 Shape                | Geometry   | 0     | 2743      
 GNIS_FEATU           | Integer    | 4     | 2693      
 FEATURE_NA           | String     | 254   | 2249      
 CTU_CLASS            | String     | 25    | 3         
 COUNTY_GNI           | Integer    | 4     | 87        
 COUNTY_COD           | String     | 2     | 87        
 COUNTY_NAM           | String     | 100   | 87        
 POPULATION           | Integer    | 4     | 1312      
 SHAPE_Leng           | Double     | 8     | 2743      
 Shape_Length         | Double     | 8     | 2743      
 Shape_Area           | Double     | 8     | 2743      
 Combined_GNIS_IDs    | String     | 50    | 2743      

Total Record Count: 2744, Total Geometry Count(by WKT): 2743
Warning 1 Potential Duplicate Features Detected in the Feature Class

In [ ]:

# Number of unique geometries or shapes in the layer (polygons in this case) matches total record count
    # Each row represents a unique CTU location, as evident by the number of unique geometries and lack of overlap between them
    # Combined_GNIS_ID is the corresponding unique attribute field marking each of those unique polygons

In [12]:

cat.open_arcgis_documentation()

Out[12]:

Approach¶

Add Field --> Summarize Data to get the Latest Sampling Date by well and polygon

In [ ]:

## We have already spatially joined a point layer of ground water well sites and a polygon layer of CTU locations
# Next identify unique well sites within each CTU 
# Define a unique well:
#   A spatially unique well is a unique combination of CTU Feature ID (GNIS_FEATURE_ID) + County ID (COUNY_GNIS_ID) + Well ID (sys_loc_co)

## We need to 
# 1) Create a new field that will mark each each spatially unique well within each polygon of the CTU layer
# 2) Join modified CTU layer with the Wells Data
# 3) Create a table or view of the joined dataset that selects only the most recent sampling date for each unique well

In [15]:

import custom_arcpy_tools as cat

cat.listFC_dataset_Attributes(wksp)

Workspace is set here c:\Projects\my_git_pages_website\Py-and-Sky-Labs\content\ArcPY\02 Result\MN_WaterContamination.gdb

Feature class: , atlas_contamination_areas, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPolygon, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , atlas_contamination_area_lines, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPolyline, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , atlas_contamination_sites, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPoint, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , atlas_contamination_sites_noca, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPoint, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , atlas_contamination_source_areas, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPolygon, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , atlas_contamination_wells_result_sum, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPoint, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , MN_PFAS_Levels, Spatial Reference name :, GCS_North_American_1983, Spatial Ref Type : Geographic, GeometryPoint, WKID: 4269
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , MN_PFAS_Levels_proj, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPoint, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , MN_PFAS_Levels_by_City, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPolygon, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , city_township_unorg, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPolygon, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , MN_WaterContamination_SpatialJoin, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPoint, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , MN_latest_WaterSamples_Joined_lyr, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPoint, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , overlap_check1, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPolygon, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 

Feature class: , overlap_check2, Spatial Reference name :, NAD_1983_UTM_Zone_15N, Spatial Ref Type : Projected, GeometryPolygon, WKID: 26915
------------------------------------------------------------------------------------------------------------------------------------------------------ 


14 Feature classes found in the workspace path

Warning Different wkid detected among feature Classes:  {26915, 4269}
Warning Different Coordinate Systems types found among Feature Classes:  {'Projected', 'Geographic'}

Add a Join Key Field¶

In [13]:

fc = "MN_WaterContamination_SpatialJoin"

for field in arcpy.ListFields(fc):
    print(field.name)

OBJECTID
Shape
Join_Count
TARGET_FID
facility_c
facility_i
facility_t
facility_n
sys_loc_co
loc_name
loc_type
loc_type_2
well_statu
geologic_u
depth_of_w
bottom_of_
surf_elev
latest_sam
number_of_
result_cou
non_detect
detect_bel
exceedance
exceedan_1
exceedan_2
exceedan_3
exceedan_4
exceedan_5
exceedan_6
detect_b_1
detect_b_2
detect_b_3
detect_b_4
detect_b_5
latitude
longitude
eda_url
data_downl
GNIS_FEATU
FEATURE_NA
CTU_CLASS
COUNTY_GNI
COUNTY_COD
COUNTY_NAM
POPULATION
SHAPE_Leng
Combined_GNIS_IDs
JoinKey

In [14]:

fc = "MN_WaterContamination_SpatialJoin"
newfield= "JoinKey"

try:
    if newfield not in [field.name for field in arcpy.ListFields(fc)]:

        ## a new field
        arcpy.management.AddField(fc,newfield,"TEXT",field_length=50)
        print(f"New {newfield} field added")
        with arcpy.da.UpdateCursor(fc,["sys_loc_co","Combined_GNIS_IDs", newfield]) as cursor:
            for row in cursor:
                well_id= str(row[0])
                combined_id= str(row[1])
                row[2] =  f"{well_id}_{combined_id}"
                cursor.updateRow(row)
        print("Successfully updated the table for the new field")
    else:
        print(f"{newfield} aleady exists in the target dataset {fc}")
        
except Exception as e:
    print("Error occured", e , type(e).__name__)
except arcpy.ExecuteError:
    print("Arcpy Error", arcpy.GetMessages(2))

JoinKey aleady exists in the target dataset MN_WaterContamination_SpatialJoin

Summary Statistics¶

`arcpy.analysis.Statistics()`¶

Calculate summary statistics for fields in a table
arcpy.analysis.Statistics(in_table, out_table, statistics_fields, {case_field}, {concatenation_separator})
- Key Arguments:
  - in_table: input table with fields to be used to compute statistics
  - out_table: output table to store the calculated results
  - statistics_field: Field or fields containing attribute values to be used and statistic to be used
    - Null values are excluded
    - Can do multiple combinations of statistics and fields
    - Format: [[field, {statistic_type}],...]
  - case_field: Field or combination of fields that will be used to group the data when calculating statistics
    - if not specified the statistics are calculated for all records
    - the output will be a single record

In [ ]:

## Calculate summary stats

fc = "MN_WaterContamination_SpatialJoin"
newfield= "JoinKey" # Field to be used to summarize or group the data
statisticfield= "latest_sam" # date field

try:
    ### TOOL: arcpy.analysis.Statistics(in_table, out_table, statistics_fields, {case_field}, {concatenation_separator})
    arcpy.analysis.Statistics(fc, "MN_WaterContamination_Stats",[[statisticfield,"MAX"]], case_field= newfield)
    print("Statistics Successfully calculated for", fc)
except Exception as e:
    print("Error occured", e, type(e).__name__)
except arcpy.ExecuteError:
    print("Arcy Error", arcpy.GetMessages(2))

Statistics Successfully calculated for MN_WaterContamination_SpatialJoin

In [17]:

arcpy.Exists("MN_WaterContamination_Stats")

Out[17]:

True

In [18]:

fc = "MN_WaterContamination_Stats"
record_count = arcpy.management.GetCount(fc)[0]

print("Record Count", record_count)

Record Count 11275

In [ ]:

### 11275 is the count of records for unique wells within each CTU by their most recent sampling date
    # This is the level or grain our datasets needs to be in to map well water site contamination by CTU 
    # one row = a unique well within a unique CTU for the most recent sampling date

Join Statistics Table to Spatial Data Layer¶

`arcpy.management.AddJoin()`¶

Join a layer to another layer or to a table based on a common field
- supports Feature layers, table views, and raster layers with attribute tables
- Result is temporary and must be saved (e.g., CopyFeatures, ExportFeatures, ExportTable)
arcpy.management.AddJoin(in_layer_or_view, in_field, join_table, join_field, {join_type}, {index_join_fields}, {rebuild_index}, {join_operation})

In [ ]:

## We have already spatially joined a point layer of ground water well sites and a polygon layer of CTU locations
# Next identify unique well sites within each CTU 
# Define a unique well:
#   A spatially unique well is a unique combination of CTU Feature ID (GNIS_FEATURE_ID) + County ID (COUNY_GNIS_ID) + Well ID (sys_loc_co)

## We need to 
# 1) Create a new field that will mark each each spatially unique well within each polygon of the CTU layer
# 2) Join modified CTU layer with the Wells Data
# 3) Create a table or view of the joined dataset that selects only the most recent sampling date for each unique well


fc = "MN_WaterContamination_SpatialJoin"
newfield= "JoinKey"

try:
    if newfield not in [field.name for field in arcpy.ListFields(fc)]:

        ## a new field
        arcpy.management.AddField(fc,newfield,"TEXT",field_length=50)
        print(f"New {newfield} field added")
        with arcpy.da.UpdateCursor(fc,["sys_loc_co","Combined_GNIS_IDs", newfield]) as cursor:
            for row in cursor:
                well_id= str(row[0])
                combined_id= str(row[1])
                row[2] =  f"{well_id}_{combined_id}"
                cursor.updateRow(row)
        print("Successfully updated the table for the new field")
    else:
        print(f"{newfield} aleady exists in the target dataset {fc}")
        
except Exception as e:
    print("Error occured", e , type(e).__name__)
except arcpy.ExecuteError:
    print("Arcpy Error", arcpy.GetMessages(2))


## Calculate summary stats

fc = "MN_WaterContamination_SpatialJoin"
newfield= "JoinKey" # Field to be used to summarize or group the data
statisticfield= "latest_sam" # date field

try:
    ### TOOL: arcpy.analysis.Statistics(in_table, out_table, statistics_fields, {case_field}, {concatenation_separator})
    arcpy.analysis.Statistics(fc, "MN_WaterContamination_Stats",[[statisticfield,"MAX"]], case_field= newfield)
    print("Statistics Successfully calculated for", fc)
except Exception as e:
    print("Error occured", e, type(e).__name__)
except arcpy.ExecuteError:
    print("Arcy Error", arcpy.GetMessages(2))

### 11275 is the count of records for unique wells within each CTU by their most recent sampling date
    # This is the level or grain our datasets needs to be in to map well water site contamination by CTU 
    # one row = a unique well within a unique CTU for the most recent sampling date
### Join two tables based on a common field and its values

in_fc= "MN_WaterContamination_SpatialJoin" # input table
in_join_key= "JoinKey"

join_table = "MN_WaterContamination_Stats" # Join table
join_table_key= "JoinKey"

try:
    ### TOOL: arcpy.management.AddJoin(in_layer_or_view, in_field, join_table, join_field, {join_type}, {index_join_fields}, {rebuild_index}, {join_operation})
    arcpy.management.AddJoin(in_fc,in_join_key, join_table, join_table_key)

    print(f"Successfully Joined {in_fc} and {join_table}")
    print(arcpy.GetMessages(0))
    
except arcpy.ExecuteError:
    print("Arcy Error: ", arcpy.GetMessages(2))
except Exception as e:
    print("Error occured", e, type(e).__name__)
### We added the Max Date from the stats table to our joined data, now we will use this to refine our joined dataset
    # if the Joined Layer has different dates for each well, then we should be able to filter it using this max date field
    # would give us the grain of data we need: one row is one unique well within one unique CTU for the latest sampling date

Successfully Joined MN_WaterContamination_SpatialJoin and MN_WaterContamination_Stats
Start Time: Friday, March 28, 2025 4:56:50 PM
Succeeded at Friday, March 28, 2025 4:56:50 PM (Elapsed Time: 0.06 seconds)

In [ ]:

### We added the Max Date from the stats table to our joined data, now we will use this to refine our joined dataset
    # if the Joined Layer has different dates for each well, then we should be able to filter it using this max date field
    # would give us the grain of data we need: one row is one unique well within one unique CTU for the latest sampling date

Select the latest water sampling dates using Select By Attributes¶

`arcpy.management.SelectLayerByAttribute()`¶

Adds updates, or removes a selection based on an attribute query
Only works on feature layers or table views!
If the query is present only those features will be used in downstream processes
arcpy.management.SelectLayerByAttribute(in_layer_or_view, {selection_type}, {where_clause}, {invert_where_clause})

In [58]:

in_fc= "MN_WaterContamination_SpatialJoin" # input table
out_lyr= "MN_WaterContamination_joined_layer"
# 
### Join two tables based on a common field and its values
in_join_key= "JoinKey"

join_table = "MN_WaterContamination_Stats" # Join table
join_table_key= "JoinKey"

try:
    ### TOOL: arcpy.management.AddJoin(in_layer_or_view, in_field, join_table, join_field, {join_type}, {index_join_fields}, {rebuild_index}, {join_operation})
    arcpy.management.MakeFeatureLayer(in_fc, out_lyr)
    
    arcpy.management.AddJoin(out_lyr,in_join_key, join_table, join_table_key)
    
    print(f"Successfully Joined {out_lyr} and {join_table}")
    print(arcpy.GetMessages(0))
    where_c = "latest_sam = MAX_latest_sam"
    arcpy.management.SelectLayerByAttribute(out_lyr, "NEW_SELECTION", where_c)
    print(f"Successfully selected features in {out_lyr}")
    print(arcpy.GetMessages(0))

    arcpy.management.CopyFeatures(out_lyr, "MN_latest_WaterSamples_Joined_lyr")
    print(f"Successfully copied features in {out_lyr}")
    print(arcpy.GetMessages(0))
    
except arcpy.ExecuteError:
    print("Arcy Error: ", arcpy.GetMessages(2))
except Exception as e:
    print("Error occured", e, type(e).__name__)

Successfully Joined MN_WaterContamination_joined_layer and MN_WaterContamination_Stats
Start Time: Friday, March 28, 2025 7:13:36 PM
Succeeded at Friday, March 28, 2025 7:13:36 PM (Elapsed Time: 0.07 seconds)
Successfully selected features in MN_WaterContamination_joined_layer
Start Time: Friday, March 28, 2025 7:13:37 PM
Succeeded at Friday, March 28, 2025 7:13:37 PM (Elapsed Time: 0.52 seconds)
Successfully copied features in MN_WaterContamination_joined_layer
Start Time: Friday, March 28, 2025 7:13:37 PM
Succeeded at Friday, March 28, 2025 7:13:38 PM (Elapsed Time: 1.24 seconds)

Check Results¶

Attribute table from a dataset in ArcGIS Pro generated by spatial join and enriched for additional attributes

The attribute table of a joined dataset comprised of water well points enriched with latest Date for each well ID
Our goal is to have one row for each unique well site with its latest sampling date

In [ ]:

### When we added the Max Sampling date, did we change the grain of the well water point layer from unique well IDs to unique wells IDs within a CTU polygon?
## No.

# What about those well IDs with multiple records of the same dates for the same CTU?

In [ ]:

## We could have gone a different way and deduplicated the dataset before adding the table, giving us the desired result of one row per unique well
# But for learning purposes our process is more informative
# Now let's understand a bit more about the data and clean up this result

In [ ]:

### check duplicate dates for each unique well id.
from collections import Counter


fc= "MN_WaterContamination_SpatialJoin" # we could have used the newly created layer 


dup_date_counter = Counter() # Counter dictionary enables us to collect and count the number of times unique keys appear (key:key_counts)

field_name= "sys_loc_co"


with arcpy.da.SearchCursor(fc, [field_name,"JoinKey","latest_sam"]) as cursor:
    for row in cursor:
        key = (row[0],row[1],row[2])
        dup_date_counter[key]+= 1 # count each time we see the same well id + ctu_id + date

count_printedrows = 0 # counter to limit printed duplicates

for id,count in dup_date_counter.items():
     
     ##  print the well_id_ctu_id_latest_sample date combos that appear more than once
     if count >1 and count_printedrows <20:
        count_printedrows +=1
        print(f"Unique well CTU date combo: {id}, appears {count} times")

Unique well CTU date combo: ('257221', '257221_2393644_659526', datetime.datetime(2007, 3, 22, 0, 0)), appears 2 times
Unique well CTU date combo: ('MNPCAP17884', 'MNPCAP17884_2397369_659526', datetime.datetime(2008, 3, 18, 0, 0)), appears 2 times
Unique well CTU date combo: ('122194', '122194_2393954_659447', datetime.datetime(2024, 6, 5, 0, 0)), appears 2 times
Unique well CTU date combo: ('1000029195', '1000029195_2395589_659526', datetime.datetime(2022, 8, 22, 0, 0)), appears 2 times
Unique well CTU date combo: ('542025', '542025_2393887_659526', datetime.datetime(2019, 5, 16, 0, 0)), appears 2 times
Unique well CTU date combo: ('175188', '175188_2395589_659526', datetime.datetime(2006, 5, 25, 0, 0)), appears 3 times
Unique well CTU date combo: ('457693', '457693_665966_659526', datetime.datetime(2022, 4, 21, 0, 0)), appears 2 times
Unique well CTU date combo: ('257614', '257614_2393644_659526', datetime.datetime(2022, 1, 25, 0, 0)), appears 2 times
Unique well CTU date combo: ('710101', '710101_2395589_659526', datetime.datetime(2024, 8, 15, 0, 0)), appears 2 times
Unique well CTU date combo: ('257704', '257704_664354_659526', datetime.datetime(2019, 5, 28, 0, 0)), appears 2 times
Unique well CTU date combo: ('593137', '593137_665966_659526', datetime.datetime(2022, 6, 2, 0, 0)), appears 3 times
Unique well CTU date combo: ('645631', '645631_665966_659526', datetime.datetime(2021, 10, 25, 0, 0)), appears 2 times
Unique well CTU date combo: ('761657', '761657_665966_659526', datetime.datetime(2019, 9, 24, 0, 0)), appears 2 times
Unique well CTU date combo: ('340455', '340455_2395345_659472', datetime.datetime(2019, 10, 31, 0, 0)), appears 2 times
Unique well CTU date combo: ('257279', '257279_2395589_659526', datetime.datetime(2006, 4, 27, 0, 0)), appears 2 times
Unique well CTU date combo: ('257929', '257929_2395589_659526', datetime.datetime(2020, 6, 18, 0, 0)), appears 2 times
Unique well CTU date combo: ('142338', '142338_663529_659526', datetime.datetime(2019, 6, 3, 0, 0)), appears 2 times
Unique well CTU date combo: ('437819', '437819_2395589_659526', datetime.datetime(2018, 9, 14, 0, 0)), appears 3 times
Unique well CTU date combo: ('609276', '609276_2394130_659449', datetime.datetime(2020, 9, 30, 0, 0)), appears 3 times
Unique well CTU date combo: ('142313', '142313_2393887_659526', datetime.datetime(2023, 2, 28, 0, 0)), appears 2 times

In [ ]:

# check the number of records we expect in the new table
cat.showFieldinfo(fc)

Summary of Fields in Feature Class: MN_WaterContamination_SpatialJoin

Name                 Type         Length   Unique_Values
 OBJECTID             | OID        | 4     | 17443     
 Shape                | Geometry   | 0     | 12914     
 Join_Count           | Integer    | 4     | 2         
 TARGET_FID           | Integer    | 4     | 17443     
 facility_c           | String     | 20    | 303       
 facility_i           | Double     | 8     | 303       
 facility_t           | String     | 20    | 7         
 facility_n           | String     | 60    | 303       
 sys_loc_co           | String     | 20    | 11246     
 loc_name             | String     | 80    | 3795      
 loc_type             | String     | 20    | 2         
 loc_type_2           | String     | 50    | 18        
 well_statu           | String     | 20    | 6         
 geologic_u           | String     | 254   | 73        
 depth_of_w           | Double     | 8     | 1145      
 bottom_of_           | Double     | 8     | 1202      
 surf_elev            | String     | 20    | 6199      
 latest_sam           | Date       | 8     | 2305      
 number_of_           | Double     | 8     | 106       
 result_cou           | Double     | 8     | 2167      
 non_detect           | Double     | 8     | 313       
 detect_bel           | Double     | 8     | 54        
 exceedance           | Double     | 8     | 34        
 exceedan_1           | Double     | 8     | 2098      
 exceedan_2           | String     | 10    | 1806      
 exceedan_3           | String     | 254   | 945       
 exceedan_4           | String     | 254   | 1156      
 exceedan_5           | String     | 254   | 4778      
 exceedan_6           | String     | 254   | 4742      
 detect_b_1           | String     | 254   | 2971      
 detect_b_2           | String     | 254   | 3190      
 detect_b_3           | String     | 254   | 6499      
 detect_b_4           | String     | 254   | 5210      
 detect_b_5           | String     | 10    | 2086      
 latitude             | Double     | 8     | 12313     
 longitude            | Double     | 8     | 12595     
 eda_url              | String     | 80    | 11246     
 data_downl           | String     | 121   | 11246     
 GNIS_FEATU           | Integer    | 4     | 276       
 FEATURE_NA           | String     | 254   | 269       
 CTU_CLASS            | String     | 25    | 3         
 COUNTY_GNI           | Integer    | 4     | 73        
 COUNTY_COD           | String     | 2     | 73        
 COUNTY_NAM           | String     | 100   | 73        
 POPULATION           | Integer    | 4     | 262       
 SHAPE_Leng           | Double     | 8     | 277       
 Combined_GNIS_IDs    | String     | 50    | 277       
 JoinKey              | String     | 50    | 11275     

Total Record Count: 17443, Total Geometry Count(by WKT): 12914

In [ ]:

### What happened?
### 11246 unique records exist by well ID alone but we know there are more unique wells if we combine well ID + CTU location ID
# We expected 11275 records because their are 11275 unique joinkey values
    # suggests duplicate well IDs

In [ ]:

from collections import defaultdict

well_ctu = defaultdict(set) # stores each well id with a unique CTU ID

fc= "MN_WaterContamination_SpatialJoin"


with arcpy.da.SearchCursor(fc, ["sys_loc_co", "Combined_GNIS_IDS"])as cursor:
    for row in cursor:
        well_id = row[0]
        ctu_id = row[1]
        # use the well id as key
        well_ctu[well_id].add(ctu_id) # store the CTU ids , if well_id exist with different CTUs they will be added 

print("Wells that appear in multiple CTUS\n")

well_ctu_counter = 0

for well_id, ctu in well_ctu.items():
    if len(ctu)>1:
        print(f"{well_id} found in {len(ctu)} CTU list here\n:{sorted(ctu)}")
        well_ctu_counter += 1

print(f"\n{well_ctu_counter} wells found associated with multiple CTUs")

Wells that appear in multiple CTUS

1000026217 found in 2 CTU list here
:['665784_659454', '665935_659523']
1000026215 found in 2 CTU list here
:['663946_659476', '665784_659454']
122002 found in 2 CTU list here
:['663529_659526', '665966_659526']
649680 found in 2 CTU list here
:['2393887_659526', '2397369_659526']
1000026216 found in 2 CTU list here
:['2393954_659447', '665784_659454']
480304 found in 2 CTU list here
:['2393644_659526', '663965_659526']
1000026235 found in 2 CTU list here
:['2393954_659447', '665784_659454']
159537 found in 2 CTU list here
:['2393887_659526', '2397369_659526']
1000026236 found in 2 CTU list here
:['2393954_659447', '665784_659454']
1000022122 found in 2 CTU list here
:['2394568_662850', '665966_659526']
544438 found in 2 CTU list here
:['2393644_659526', '2397369_659526']
707094 found in 2 CTU list here
:['2394963_659526', '2395589_659526']
617693 found in 2 CTU list here
:['663529_659526', '665966_659526']
257205 found in 2 CTU list here
:['2393644_659526', '665966_659526']
1000022134 found in 2 CTU list here
:['2393472_659464', '665784_659454']
1000023198 found in 2 CTU list here
:['2395609_659526', '2395610_659526']
649697 found in 2 CTU list here
:['2393644_659526', '664354_659526']
162973 found in 2 CTU list here
:['663529_659526', '665966_659526']
257250 found in 2 CTU list here
:['2393644_659526', '2395854_659514']
1000026213 found in 2 CTU list here
:['664979_659493', '665784_659454']
569775 found in 2 CTU list here
:['2395589_659526', '665966_659526']
1000020749 found in 2 CTU list here
:['2393954_659447', '2394747_659464']
809769 found in 2 CTU list here
:['2395589_659526', '664354_659526']
245657 found in 2 CTU list here
:['2395282_659447', '664701_663198']
1000026219 found in 2 CTU list here
:['665235_659525', '665784_659454']
501969 found in 2 CTU list here
:['2394417_659472', '2395350_659472']
1000026210 found in 2 CTU list here
:['664601_659485', '665784_659454']
430924 found in 2 CTU list here
:['2394963_659526', '2395589_659526']
430920 found in 2 CTU list here
:['2397369_659526', '664386_659464']

29 wells found associated with multiple CTUs

In [45]:

# check why there are less unique sys loc codes then unique well id+ CTU id
fc= "MN_WaterContamination_SpatialJoin"

row_counter = 0

with arcpy.da.SearchCursor(fc, ["JoinKey", "latest_sam"], sql_clause=["DISTINCT" , "ORDER BY SYS_LOC_CO"]) as cursor:
    for row in cursor:
            row_counter +=1

print(f"{row_counter} total distinct rows")

11275 total distinct rows

In [ ]:

### We found a few things that indicate the need to deduplicate the data by combination of well_id + CTU_id
# 1. Same well ids are in different CTU polygons
# 2. Same well id has same date repeated multiple tmes
# 3. The Joinkey retains all the distinct well data by CTU including latest sampling date

Deduplicate Dataset and Create a new Feature Class¶

Situation:
- Our goal here is to map the maximum reported contamintant exceedance level across all well for each CTU polygon
  - Map the distribution of contaminant levels exceeding guidelines as a color gradient for ease of comparison across all CTUs
- The MN Ground Water Contamination Altas Summary dataset contains duplicate well sampling sites with the same values for: well ID, collection date, contaminants, etc
Task:
- Reduce this dataset down to the only the critical data to display max contaminants for unique wells
- Create a new feature class where one row = a unique well within a unique CTU for the most recent sampling date #### CreateFeatureClass()
Creates a new feature class (fc) in a geodatabase or shape file in a folder
arcpy.management.CreateFeatureclass(out_path, out_name, {geometry_type}, {template}, {has_m}, {has_z}, {spatial_reference}, {config_keyword}, {spatial_grid_1}, {spatial_grid_2}, {spatial_grid_3}, {out_alias}, {oid_type})
Key fields
- out_path: the enterprise or file geodatabase or folder to store the output
- out_name: name of the fc to be created
- geometry_type: specify the geometry of the output fc (e.g., POINT, MULTIPOINT, POLYGON, POLYLINE, MULTIPATCH (use for 3D features))
- template: an existing feature class or list of fc to define the attribute fields of the new feature class (see example below)
- spatial_reference: can specify with .prj file, path to a feature class, or defining a spatial reference (optional but if no spatial ref provided will be undefined)

In [ ]:

# List the fields in the fc
fc=  "MN_WaterContamination_SpatialJoin" # input layer

print([f.name for f in arcpy.ListFields(fc)]) # list comprehension extracting all field names

['OBJECTID', 'Shape', 'Join_Count', 'TARGET_FID', 'facility_c', 'facility_i', 'facility_t', 'facility_n', 'sys_loc_co', 'loc_name', 'loc_type', 'loc_type_2', 'well_statu', 'geologic_u', 'depth_of_w', 'bottom_of_', 'surf_elev', 'latest_sam', 'number_of_', 'result_cou', 'non_detect', 'detect_bel', 'exceedance', 'exceedan_1', 'exceedan_2', 'exceedan_3', 'exceedan_4', 'exceedan_5', 'exceedan_6', 'detect_b_1', 'detect_b_2', 'detect_b_3', 'detect_b_4', 'detect_b_5', 'latitude', 'longitude', 'eda_url', 'data_downl', 'GNIS_FEATU', 'FEATURE_NA', 'CTU_CLASS', 'COUNTY_GNI', 'COUNTY_COD', 'COUNTY_NAM', 'POPULATION', 'SHAPE_Leng', 'Combined_GNIS_IDs', 'JoinKey']

In [ ]:

## We can use this approach to acess fields, add fields and to update the order
fc=  "MN_WaterContamination_SpatialJoin" # input layer

# extracrt all the fields in the fc into  a list
orig_fields = [f.name for f in arcpy.ListFields(fc)]

# Move a field to the front of a new list
reorg_list = ["facility_c"]+ [f for f in orig_fields if f != 'facility_c'] 

print([reorg_list])

[['facility_c', 'OBJECTID', 'Shape', 'Join_Count', 'TARGET_FID', 'facility_i', 'facility_t', 'facility_n', 'sys_loc_co', 'loc_name', 'loc_type', 'loc_type_2', 'well_statu', 'geologic_u', 'depth_of_w', 'bottom_of_', 'surf_elev', 'latest_sam', 'number_of_', 'result_cou', 'non_detect', 'detect_bel', 'exceedance', 'exceedan_1', 'exceedan_2', 'exceedan_3', 'exceedan_4', 'exceedan_5', 'exceedan_6', 'detect_b_1', 'detect_b_2', 'detect_b_3', 'detect_b_4', 'detect_b_5', 'latitude', 'longitude', 'eda_url', 'data_downl', 'GNIS_FEATU', 'FEATURE_NA', 'CTU_CLASS', 'COUNTY_GNI', 'COUNTY_COD', 'COUNTY_NAM', 'POPULATION', 'SHAPE_Leng', 'Combined_GNIS_IDs', 'JoinKey']]

In [ ]:

### OPTIONAL STEP
# Create a  layer where we keep rows where the joinkey is unique 
    # this is more for analytic and QC purposes that we didn't lose data
# one row = one unique well id within one unique CTU

fc=  "MN_WaterContamination_SpatialJoin" # input layer
outfc= "wells_per_CTU_unique"


# select fields without having to hardcode eachone
fields = [f.name for f in arcpy.ListFields(fc) ]
fields = ["JoinKey"] + [f for f in fields if f != "JoinKey"] # place the joinkey upfront

# Order by joinkey and lates_sam date descending
sql_c= (None, "ORDER BY JoinKey, latest_sam DESC")


# will store unique records as a dictionary with joinkey (well_id, CTU_ID) as the key and all row data as the values
unique_records = {}

with arcpy.da.SearchCursor(fc,fields, sql_clause=sql_c) as cursor:
    for row in cursor:
        join_id = row[0]

        # check unique records dictionary, keeps the first joinkey to avoid duplicates
        if join_id not in unique_records:
            # store the entire row
            unique_records[join_id]= row  # store corresponding row and join key

sr= arcpy.Describe(fc).spatialReference

### TOOL: arcpy.management.CreateFeatureclass(out_path, out_name, {geometry_type}, {template}, {has_m}, {has_z}, {spatial_reference}, {config_keyword}, {spatial_grid_1}, {spatial_grid_2}, {spatial_grid_3}, {out_alias}, {oid_type})
arcpy.management.CreateFeatureclass(wksp, outfc, "POINT", template=fc, spatial_reference= sr) # template copies the fc schema-field types, lengh, order

try:
    with arcpy.da.InsertCursor(outfc, fields) as icursor:
        for row in unique_records.values():
            icursor.insertRow(row)


except Exception as e:
    print("Error occured", e, type(e).__name__)

except arcpy.ExecuteError:
    print("Arcpy Error", arcpy.GetMessages(2))

In [40]:

cat.showFieldinfo(outfc)

*Summary*
 Fields in the Point Feature Class: wells_per_CTU_unique

Name                 Type         Length   Unique_Values
 OBJECTID             | OID        | 4     | 11275     
 Shape                | Geometry   | 0     | 11168     
 Join_Count           | Integer    | 4     | 2         
 TARGET_FID           | Integer    | 4     | 11275     
 facility_c           | String     | 20    | 297       
 facility_i           | Double     | 8     | 297       
 facility_t           | String     | 20    | 7         
 facility_n           | String     | 60    | 297       
 sys_loc_co           | String     | 20    | 11246     
 loc_name             | String     | 80    | 3366      
 loc_type             | String     | 20    | 2         
 loc_type_2           | String     | 50    | 18        
 well_statu           | String     | 20    | 6         
 geologic_u           | String     | 254   | 68        
 depth_of_w           | Double     | 8     | 1106      
 bottom_of_           | Double     | 8     | 1186      
 surf_elev            | String     | 20    | 5589      
 latest_sam           | Date       | 8     | 2305      
 number_of_           | Double     | 8     | 106       
 result_cou           | Double     | 8     | 2167      
 non_detect           | Double     | 8     | 313       
 detect_bel           | Double     | 8     | 54        
 exceedance           | Double     | 8     | 34        
 exceedan_1           | Double     | 8     | 2098      
 exceedan_2           | String     | 10    | 1806      
 exceedan_3           | String     | 254   | 945       
 exceedan_4           | String     | 254   | 1156      
 exceedan_5           | String     | 254   | 4778      
 exceedan_6           | String     | 254   | 4742      
 detect_b_1           | String     | 254   | 2971      
 detect_b_2           | String     | 254   | 3190      
 detect_b_3           | String     | 254   | 6499      
 detect_b_4           | String     | 254   | 5210      
 detect_b_5           | String     | 10    | 2086      
 latitude             | Double     | 8     | 10922     
 longitude            | Double     | 8     | 10966     
 eda_url              | String     | 80    | 11246     
 data_downl           | String     | 121   | 11246     
 GNIS_FEATU           | Integer    | 4     | 276       
 FEATURE_NA           | String     | 254   | 269       
 CTU_CLASS            | String     | 25    | 3         
 COUNTY_GNI           | Integer    | 4     | 73        
 COUNTY_COD           | String     | 2     | 73        
 COUNTY_NAM           | String     | 100   | 73        
 POPULATION           | Integer    | 4     | 262       
 SHAPE_Leng           | Double     | 8     | 277       
 Combined_GNIS_IDs    | String     | 50    | 277       
 JoinKey              | String     | 50    | 11275     

Total Record Count: 11275, Total Geometry Count(by WKT): 11168
Warning 107 Potential Duplicate Features Detected in the Feature Class

In [63]:

## Now we focus on getting the latest sample reading across all points in each CTU

# Input feature class
fc = "MN_WaterContamination_SpatialJoin"
outfc = "wells_per_CTU_latest"

# Get all field names (preserving full schema)
fields = [f.name for f in arcpy.ListFields(fc)]
fields = ["Combined_GNIS_IDs"] + [f for f in fields if f != "Combined_GNIS_IDs"]  # Ensure CTU ID is first

# Order by CTU ID and latest sampling date descending
sql_c = (None, "ORDER BY Combined_GNIS_IDs, latest_sam DESC, exceedance DESC")

# Store the most recent record for each CTU
unique_records = {}

with arcpy.da.SearchCursor(fc, fields, sql_clause=sql_c) as cursor:
    for row in cursor:
        ctu_id = row[0]  # Grouping by Combined_GNIS_IDs
        if ctu_id not in unique_records:
            unique_records[ctu_id] = row  # Keep only the first/latest well sample record per CTU

# Use original spatial reference and schema
sr = arcpy.Describe(fc).spatialReference

# Create the output feature class using a schema template
arcpy.management.CreateFeatureclass(wksp,outfc,"POINT",template=fc,spatial_reference=sr)

# Insert the filtered rows into the output feature class
try:
    with arcpy.da.InsertCursor(outfc, fields) as icursor:
        for row in unique_records.values():
            icursor.insertRow(row)
    print("Successfully created a new Feature Class", outfc, arcpy.GetMessages(0))
except Exception as e:
    print("Error occured", e, type(e).__name__)

except arcpy.ExecuteError:
    print("Arcpy Error", arcpy.GetMessages(2))

Successfully created a new Feature Class wells_per_CTU_latest Start Time: Thursday, April 3, 2025 11:11:14 PM
Succeeded at Thursday, April 3, 2025 11:11:14 PM (Elapsed Time: 0.19 seconds)

Spatial Join Analysis¶

`arcpy.analysis.SpatialJoin()`¶

arcpy.analysis.SpatialJoin(target_features, join_features, out_feature_class, {join_operation}, {join_type}, {field_mapping}, {match_option}, {search_radius}, {distance_field_name}, {match_fields})
- Target features : features to be enriched with attributes
- Join features: features contributing attribute data
Key optional arguments:
- join operation: Type of join either JOIN_ONE_TO_ONE or JOIN_ONE_TO_MANY. default one-to-one
- join type : Number of target records to keep either 'KEEP_ALL'(default, keep all records) or 'KEEP_COMMON' (only matched records are kept). Default keeps all fields
- field mapping: Use this to customize output layer fields, add, delete, rename, re-order,change properties, and combine fields
- match option: Intersect, Within, Contains, Closest, Closest_Geodesic,Are_identical_to, Largest_overlap...and more!
- search_radius: distance to search for matches
- match fields: pairs of fields to restrict the join to only those fields with matching values

In [64]:

### Goal: Create the final layer with the latest well data for the max exceedance value
    # Will show the max contaminanat levels for each CTU based on the latest well data
    # Spatial join to enrich the CTU layer with the latest max contaminant data


                                    
ctu_polygons = "city_township_unorg" # target features
latest_points = "wells_per_CTU_latest"# join features
out_fc = "CTU_latest_sample_summary"

try:
                         ## target         # join features
    arcpy.analysis.SpatialJoin(ctu_polygons,latest_points,out_fc,join_operation="JOIN_ONE_TO_ONE",join_type="KEEP_ALL",match_option="INTERSECT")
    print("Successfully joined the features", arcpy.GetMessages(0))
    
except Exception as e:
    print("Error occured", e, type(e).__name__)

except arcpy.ExecuteError:
    print("Arcpy Error", arcpy.GetMessages(2))

Successfully joined the features Start Time: Thursday, April 3, 2025 11:11:17 PM
Succeeded at Thursday, April 3, 2025 11:11:17 PM (Elapsed Time: 0.51 seconds)

Savage, MN GIS

Other articles

QC and Join data in Arcpy

What is Arcpy?:¶

What is ArcPY? Part3

What is Arcpy?: Part 3¶

What is ArcPY? Part2

What is Arcpy?: Part 2¶

What is ArcPY? Part1

Other articles

What is Arcpy?:¶

What is Arcpy?: Part 3¶

What is Arcpy?: Part 2¶

links

social