Post 2: Python Scripting (3 Posts)

Pyscripting 1: March 15, 2015                                                     

                                                Introduction

The purpose of this exercise is to expose the student to Python Scripting (PS) as it pertains to tools found in ArcGIS.

According to ArcGIS Help, PS is a free programming language that can be integrated into a number of computer software programs, like ArcGIS products. Python can be used to create batch-processing of several datasets at one time. For example, in this Python scripting (PS) exercise several raster datasets were processed in succession using PS in conjunction with ArcMap tools, like projection and extract.

In addition to its capability to process numerous datasets at once, PS is a powerful tool because it is a used by numerous users of varying skill levels. As such, PS tutorials and information are found at several websites. In fact, ArcGIS Help pages for Geodatabase toolsets usually have a section that gives an example of PS for each tool (fig. 1).

Following are several websites that are useful for those interested in learning more about Python Scripting:

• Python.org: https://www.python.org/

                              https://wiki.python.org/moin/BeginnersGuide/NonProgrammers

• PennState.edu: https://www.e-education.psu.edu/geog485/node/17

• After Hours Programming:http://www.afterhoursprogramming.com/tutorial/Python/Overview/

Figure 1. Python Scripting code sample for loading raster sets to a geodatabase. From AcrGIS Help 10.1 (http://resources.arcgis.com/en/help/main/10.1/index.html#//001200000025000000).

                                                    Methods

In order to create the PS for this exercise (fig. 2), several steps were followed:

• Setting  the environment defined the workspace (i.e. folders) where the data to be processed originated from:

#set environment

env.workspace ="Q:\StudentCoursework\CHupy\g2155.GEOG.337.001\XXX\Ex_5\Tremp_raster2"

arcpy.env.overwriteOutput = True

• ListRasters returns a list of rasters available in the workspace (defined above) based on name and/or raster type (e.g. tiff, jpeg, etc.):

#Get and print a list all rasters in workspace

listofrasters = arcpy.ListRasters()

for raster in listofrasters:

    print(raster)

• Before projecting the rasters, they needed to be defined:

#set up output variable

    rasterOut = '{}_Out.tif'.format(raster)

    rasterExtract = '{}_Extract.tif'.format(raster)

• Rasters were projected into the appropriate coordinate system:

#project raster

    arcpy.ProjectRaster_management(raster, rasterOut, "Q:\StudentCoursework\CHupy\g2155.GEOG.337.001\XXX\Lab\Ex_5\Tremp_raster2\TMP.gdb\Fema\S_BASE_INDEX")

• Rasters were extracted to the appropriate extents and saved:

#extract rasters

    outExtractByMask = ExtractByMask(rasterOut, "Q:\StudentCoursework\CHupy\g2155.GEOG.337.001XXX\Lab\Ex_5\Tremp_raster2\TMP.gdb\Boundaries\County_Boundary")

  #save extract

    outExtractByMask.save(rasterExtract)

• The projected and extracted rasters were sent to the appropriate database:

#onvert rasters to TMP geodatabase

    arcpy.RasterToGeodatabase_conversion(outExtractByMask, "Q:\StudentCoursework\CHupy\g2155.GEOG.337.001\XXXLab\Ex_5\Tremp_raster2\TMP.gdb")

Figure 2. Screen-shot of the Python Script created for the exercise.

                                                      Conclusion

Once the rasters were projected, extracted, and saved to the appropriate geodatabase, they were used as part of the semester project for my geography 337 course (see Post 3: 'Data Collection').

 

Pyscripting 2: April 21, 2015

                                                   Introduction

 The Python Scripting assignment in this section correlates to Post 5: Network Analysis. Post 5 deals with the routing of frac sand from mining facilities in Wisconsin to railroad depots for shipment across the country.

  However, not all mining facilities that were geocoded in Post 4 are active. As such, a Python Script (fig. 1) will be written to only select mines that are active and save them as a layer. The Python Script will also create a layer which is comprised of feature classes that are categorized as mines, as opposed to facilities that are classified as processing facilities, for example. Finally, the script will create a feature layer of all the active mine facilities that have no railroads within 1.5 km of them.



Figure 1. Python Script for creating feature classes for active mining facilities at a distance of greater than 1.5 km from railroad depots for network analysis activity in Post 5.

                                                      Methods

• System modules were imported and the environment was set in the Python Scripting program:

 

#import system modules
import arcpy
from arcpy import env

#set environment
env.workspace= "Q:\StudentCoursework\CHupy\g2155.GEOG.337.001\LUCZAKJN\Lab\ex7\ex7.gdb"

arcpy.env.overwriteOutput = True

• Variables used in the script were defined:

 

#set variables
all = "all_mines"
active = "active_mines"
status = "status_mines"
nrmines = "mines_norail"
wi = "wi"
rails = "rails_wtm"
finalnr = "mines_norail_final"

• SQL statements were established to select:

Active Mines:
#write SQL statements to select all active mines
activeSQL = "Status" + " = " + "'Active'"

Mining Facilities:
#write SQL statement to select "mines"
minesSQL = "Facility_T" + " LIKE " + "'%Mine%'"

Mines without a rail depot:
#write SQL statement to select all mines without a rail depot
nrminesSQL = " NOT " + "Facility_T" + " LIKE " + "'%Rail%'"

• Temporary feature layer is created for each SQL above (i.e. activeSQL, minesSQL, and nrminesSQL:

Active:
#select mines that are active
arcpy.MakeFeatureLayer_management(all, active, activeSQL)

Facility:
#select mines that are active with "mine" facility field
arcpy.MakeFeatureLayer_management(active, status, minesSQL)

No Rails:
#select mines that are active with no rail depot
arcpy.MakeFeatureLayer_management(status, nrmines, nrminesSQL)

• Mines within 1.5 Km of a rail depot were removed:

 #remove mines within 1.5 miles of a rail line
arcpy.SelectLayerByLocation_management(nrmines, "WITHIN_A_DISTANCE", rails, "1.5 KILOMETER", "REMOVE_FROM_SELECTION")

• New layers selected above were saved to the ex7 geodatabase:

#save new layers of selected information
arcpy.CopyFeatures_management(nrmines, finalnr)

 

Pyscripting 2: May 13, 2015

Introduction

This python script corresponds to blog Post 6: Raster Analysis of Sand Mine Locations (Blog Post 6 ) The purpose of this python script is to create a weighted index model based on the five risk assessment models generated in the last post (Figure 1). One of the models is given a weighted value and a new risk model was generated (Figure 2).

 

Figure 1. Screenshot of the PyScripter interface.

 

Figure 2. Weighted risk assessment map generated using PyScripter.



 

 

Methods

 

• System modules were imported into the PyScripter interface, the working environment was set, and the spatial analyst extension was activated:

 

#import system modules
import arcpy
from arcpy import env

 

#set environment
env.workspace = "Q:\StudentCoursework\CHupy\g2155.GEOG.337.001\LUCZAKJN\Lab\Ex_5\Tremp_raster2\NewTMP.gdb"
arcpy.env.overwriteOutput = True

 

#Activate License
arcpy.CheckOutExtension ("spatial")

 

 

• Variables for the five rasters were specified:

 

#Set Variables
streams  =  arcpy.Raster("streamsreclassuse1")
park  =  arcpy.Raster("parkreclassuse1")
school  =  arcpy.Raster("school_reclass1")
prime  =  arcpy.Raster("farmlandreclassuse")
zone  =  arcpy.Raster("zoning_reclassuse1")

 

 

• Raster math (multiplication) was performed on the variable of interest (i.e. school):

 

#Raster Math
outweight  =  (school*1.5)

 

 

• The weighted index was created by adding cell values of each risk assessment model

 

#Weighted Index
weighted  =  (outweight + park + streams + prime + zone)

 

• The resulting index model was saved in PyScripter since rasters file are temporary by default

#Save
weighted.save("weighted_result")

 

print "the script has completed"