- Flow chart of a workflow to prepare two datasets for spatial join
- The analysis includes checks for unique ID fields and duplicate geometry
What is Arcpy?:¶
Reminder: 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 ## need this to process and extract downloaded zipfiles later on
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)
In [2]:
### BE SURE TO INSTALL MAGIC LIBRARY FIRST for filetype validation example below :
# step1 activate arcgispro-py3 environment: conda activate "C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3"
# step2 install library: pip install python-magic-bin on windows, for unix system use conda install -c conda-forge python-magic
import magic ## to validate file type requested from web
What we did last time¶
- Used requests libary to pull a shapefile dataset of ground water contamination from MN Geospatial Commons
- Used arcpy.conversion.FeatureClassToGeodatabase() to load a list of shapefiles into a geodatabase
- Created a function that would simplify checking the existance and characteristics of our feature classes
What we will do this time:¶
- Revisit the custom function to check our gdb
- Create a new function using the request library
- Pull some new data into our project
- Perform several QC checks on attributes and geometry
- Perform Spatial Analysis within Arcpy
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()
In [3]:
## Let's load in our custom tool
import custom_arcpy_tools as cat
In [ ]:
# note: When you import the custom arcpy tools module it sets a variable called initial_dir to store the path where your script is located
# this is intentionally done to support the modules functionality i.e., creating folders locally from a known reference point to store data
In [4]:
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
In [ ]:
# See the ValueError?
### Don't remember the arguments for this function?
### Use this:
# help(function_name)
In [5]:
help(cat.listFC_dataset_Attributes) ## use help(function_name) to get the documentation for the function
In [6]:
# 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
In [7]:
# list the files
%ls
In [ ]:
## should see a file ending with .gdb. Copy the folder path and add its name below
In [8]:
## 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
In [9]:
arcpy.Exists(wksp) ## same as using print statement following if not os.path.exists(gdb_path)
Out[9]:
In [10]:
# Now. Retry the function that failed earlier
cat.listFC_dataset_Attributes()
In [11]:
# Confirm we have those feature classes
arcpy.ListFeatureClasses()
Out[11]:
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 three walkthroughs
- Now will validate uniqueness of the data and join the datasets
In [12]:
# quick reminder of where our current workspace is set
print("Data Folder is here" , wksp)
# let's call it gdb for clarity
gdb=wksp
Preparing for Spatial Join¶
When joining any datasets, spatial or otherwise, it is a good practice to:
- Validate the record counts of the initial and joined datasets
Understand the grain of the input datasets
- Example: one dataset is organized by City,Townships, and Unorganized territories (CTU) assigned a unique ID called GNIS_Feature_ID
- GNIS_Feature_ID = Geographic Names Information System (GNIS) Feature ID
- A unique ID from the GNIS federal database developed by USGS
- GNIS database contains official names and GNIS_Feature_ID for geographic features (cities, townships, lakes, etc)
- Why does this matter?:
- A unique standard ID defines the grain or resolution of a dataset
- Use that ID for analysis of discrete entities (e.g., cities)
- e.g., Each GNIS_Feature_ID corresponds to a single unique location (City,Township or Unorganized Territory)
Look for unique standard ID fields to track discrete real-world geographical entities:
- GNIS_Feature_ID (USGS)
- GEOID (Census)
- FIPS Codes
- UUID, PLACE_ID
MN City,Township, Unorganized Territory Dataset¶
- Data represents boundaries of cities, townships and unorganized territories(CTU) in Minnesota
- key fields
- GNIS_Feature_ID: official ID for geographic features that relates them to their Name and location in the GNIS Database
- County_GNIS_Feature_ID: official county ID that relates counties to their Name and Location in the GNIS Database
- FEATURE_NAME: Name of the city, township or unorganized territory location
- Source : MN Geospatial Commons:
In [13]:
## 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])}")
In [14]:
## 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])}")
In [15]:
## List and summarize all fields
print("Summary of Fields")
for field in arcpy.ListFields("city_township_unorg"):
# For better display left align, pad to same length. e.g, field.name:20, means left align text and pad to 20 characters
print(f" - {field.name:15} | {field.type:10} | {field.length}")
In [16]:
### lets wrap this in a function and add total record counts
# A function to get the fields and record counts
def showFieldinfo(fc):
record_count = arcpy.management.GetCount(fc)
print("Summary of Fields")
for field in arcpy.ListFields(fc):
print( f" {field.name:15} | {field.type:15} | {field.length}")
print(f"\nRecord Count:{record_count} in Feature Class {fc}")
return
In [17]:
fc = "city_township_unorg"
showFieldinfo(fc)
In [ ]:
## We now have a way to summarize all fields in a more visual way
QC Overlapping Polygons: Do CTU polygons overlap?¶
- This is to ensure that we do have a unique polygon for each CTU
- If we do, bucketing well contamination sites by unique CTU will make sense
- otherwise, some wells would be grouped into multiple CTU making analysis more complex
Add a positive control for overlap¶
- we will grab a random polygon from CTU layer using the SearchCursor and the SHAPE@WKT field(i.e., string of polygon shape vertices coordinates)
- Then insert the polygon into the same CTU layer to create a duplicate
- This way we know we have an overlap and that our method can detect overlap
In [ ]:
# We can even access the geometry of the layer
fc = "city_township_unorg"
with arcpy.da.SearchCursor(fc, ["SHAPE@WKT"]) as cursor:
for row in cursor:
shape= row[0]
print(shape)