Automation of geoprocessing makes the manipulation of data easier, faster, and reduces user errors. A great way to automate geoprocessing workflows is through the open-source programming language Python. Python can be run inside ArcMap or in an integrated development environment (IDE) such as PythonWin or PyScripter. Data management, display, and other aspects of GIS can also be automated with Python scripts, however, the automation of geoprocessing will be the focus of this GIS II sand mining analysis. PyScripter will be used to write and execute scripts because it allows for more extensive debugging and error checking, making it a great environment for learning proper Python logic and syntax.
As each script is written, it will be posted here along with updates to the notes section below.
Notes on Python Use and Syntax:
Script 1
- For help with Python syntax and sample code, ArcGIS desktop help can be consulted.
- CODE: import arcpy
- When scripting outside of ArcMap(in an IDE), the 'arcpy' package needs to be imported to use geoprocessing tools from the ArcToolbox. Depending on the desired tools, extra extensions may also need to be imported (e.g. the spatial analysis extension)
- CODE: from arcpy import env
env.workspace = "file_path"
arcpy.env.overwriteOutput = True - A workspace sets the default location (file_path) the script will use to find files. The last line of code sets the overwrite output option to true allowing for the script to be run multiple times without having errors.
- A # indicates a comment line which will not be executed when the script is run. This allows for descriptions of the code to be written directly in the script to make it more understandable and user friendly. Comments appear in a green color in PyScripter.
- A function preforms an action and returns a value. Functions follow this basic syntax: FunctionName(parameter, parameter, etc). The value and number of parameters depend on the specific function being preformed.
- The establishment of variables allows a script to be run in a variety of scenarios. A simple example of a variable is: x = 3. In this example 'x' is the variable and can be referenced in functions as a parameter. When x is referenced, the value 3 is returned. If the value of 'x' needs to be changed, only the line of code which established the variable needs to be modified.
- When referencing variables, they are not put in quotes because they are defined in the script. For parameters that are not defined in the script (e.g. file paths), quotes are needed.
- Lists of data can be established through use of multiple functions all beginning with the word 'List' (e.g. ListDatasets() or ListFields()). After a list is established, a for-loop can be used to cycle through the list and process each element individually. Lines of code that are indented after the for-loop are included in the loop. To end the loop, the next line of code should not be indented.
- The character r can be placed in front of a file path to specify that the string of text is indeed a file path. Python could become confused with the large number of \ characters which are used for differing functions.
- When writing SQL statements via script, spaces are necessary before and/or after Boolean operators. (e.g. "Facility_T" + " LIKE " + "'%Mine%'" reads as: Facility_T LIKE '%Mine%'. whereas "Facility_T" + "LIKE" + "'%Mine%'" reads Facility_TLIKE'%Mine%'. The later cannot be understood and an error will be generated.)
Python Script 1:
For the first script of the semester, three rasters (a DEM, land cover, and cropland) will be projected, clipped (extract by mask), and loaded into a geodatabase. This geoprocessing is automated through use of a for-in-loop.
#-------------------------------------------------------------------------------
# Name: Ex5
# Purpose: Project, clip, and load rasters into a geodatabase
#
# Author: foxp
#
# Created: 10/10/2014
#-------------------------------------------------------------------------------
#import system module
import arcpy
#import spatial analyst extension
from arcpy.sa import *
#turn on the spatial analyst license
arcpy.CheckOutExtension("spatial")
#establish workspace
from arcpy import env
env.workspace = "Q:\StudentCoursework\CHupy\Geog337_f14\FOXP\Ex5\pythBase"
arcpy.env.overwriteOutput = True
#Generate list of rasters
ListOfRasters = arcpy.ListRasters()
#Cycle through rasters in the workspace, one by one, and preform each tool
for raster in ListOfRasters:
print(raster)
#Assign names for output rasters
rasterOut = '{}_Out.tif'.format(raster)
#project all rasters to the coordinate system of the geodatabase through
#referencing one of it's feature classes
arcpy.ProjectRaster_management(raster, rasterOut,
"Q:\StudentCoursework\CHupy\Geog337_f14\FOXP\Ex5\pythBase\TMP.gdb\Land\Dams")
#extract the rasters by the desired boundary
outExtractByMask = ExtractByMask(rasterOut,
"Q:\StudentCoursework\CHupy\Geog337_f14\FOXP\Ex5\pythBase\TMP.gdb\Boundaries\County_Boundary")
#Load rasters into geodatabase
arcpy.RasterToGeodatabase_conversion(outExtractByMask,
"Q:\StudentCoursework\CHupy\Geog337_f14\FOXP\Ex5\pythBase\TMP.gdb")
print("{} converted successfully:".format(raster))
print("The script has completed")
Python Script 2:
For the second script of the semester, all active mines without any relation to railways will be selected and made a new feature class. The goal of this script is to identify all mines that will need to use road networks to haul sand for a subsequent network analysis.
#-------------------------------------------------------------------------------
# Name: Ex7
# Purpose: Select active mines, exclude inactive mines, exclude mines within
# 1.5 km of a rail line
#
# Author: foxp
#
# Created: 05/11/2014
# Copyright: (c) foxp 2014
#-------------------------------------------------------------------------------
#import system modules
import arcpy
from arcpy import env
#establish enviornment settings
env.workspace = r"Q:\StudentCoursework\CHupy\Geog337_f14\FOXP\Ex7\ExSeven.gdb"
arcpy.env.overwriteOutput = True
#establish variables for the feature classes
mines = "all_mines"
studyarea = "wi"
rails = "rails_wtm"
#Establish variable to select active mines
activeSQL = "Status" + "=" + "'Active'"
#Establish variable to select all mines with 'mine' in the facility type
minesSQL = "Facility_T" + " LIKE " + "'%Mine%'"
#Establish variable to select all mines without the word 'rail' in the facility type
norailSQL = "NOT " + "Facility_T" + " LIKE " + "'%Rail%'"
#Perform the SQL variables in succession
arcpy.MakeFeatureLayer_management(mines, "active_mines", activeSQL)
arcpy.MakeFeatureLayer_management("active_mines", "status_mines", minesSQL)
arcpy.MakeFeatureLayer_management("status_mines", "norail_mines", norailSQL)
#Select all queried mines within WI
allnorailSel = arcpy.SelectLayerByLocation_management("norail_mines", "WITHIN", studyarea)
#delete all mines within 1.5 km of a rail line
arcpy.SelectLayerByLocation_management(allnorailSel, "WITHIN_A_DISTANCE", rails, 1500, "REMOVE_FROM_SELECTION" )
#Save the selection as a feature
arcpy.CopyFeatures_management("norail_mines", "norail_mines_final")
print "Script Completed"
Python Script 3:
The last script of the semester takes the reclassified rasters from the impact model in the Raster Analysis blog post and creates a weighted index. The output of the script is the same as the output of the Impact model, however one variable will be given emphasis via establishing a weight.
#-------------------------------------------------------------------------------
# Name: Ex8 Python
# Purpose: Take Ex8 Impact Results and create a weighted index
#
# Author: foxp
#
# Created: 15/12/2014
#-------------------------------------------------------------------------------
#import system modules
import arcpy
from arcpy import env
from arcpy.sa import *
#turn on the spatial analyst license
arcpy.CheckOutExtension("spatial")
#establish workspace
env.workspace = "Q:\StudentCoursework\CHupy\Geog337_f14\FOXP\Ex5\pythBase\TMP.gdb"
arcpy.env.overwriteOutput = True
#establish raster variables
streams = arcpy.Raster("streams_reclass")
primefarm = arcpy.Raster("PrimeFarm_reclass")
popzone = arcpy.Raster("popzone_reclass")
school = arcpy.Raster("school_reclass")
wildlife = arcpy.Raster("wildlife_reclass")
#establish weighted variable
outweight = (streams * 1.5)
#create weighted index by adding the variables together (same as raster calculator)
weighted = (outweight + primefarm + popzone + school + wildlife)
weighted.save("weighted_index")
print "The Script has completed"
No comments:
Post a Comment