(ArcGIS 10 for Economics Research)
mkudamatsu.github.io/gis_lecture2.html
Masayuki Kudamatsu
12 October, 2018
Press SPACE to proceed.
To go back to the previous slide, press SHIFT+SPACE.
Match features from two vector datasets by location
Henderson et al. (2018) find these ecoregions predict where economic activities take place
But we cannot easily match this data to administrative areas without spatial join
Weather data: increasingly popular among economists
But grid cells are the units of observation in weather data
Evenly spaced meridians and parallels divide the earth into grid cells
(Image source: ArcGIS Help "What are grids and graticules?")
$\Rightarrow$ Spatial Join matches grid cells with socio-economic units
1. Alsan (2015)
2. Replicate Alsan (2015)
3. Introduction to Python for ArcGIS
Did Tsetse flies cause Africa's poverty?
Important?
Original?
Feasible?
Ethnographic Atlas + Murdock (1959)'s map
Weather at 2°x2° grid cells in 1871 from 20th Century Reanalysis Project
(Figure 2 of Alsan 2015)
Feed temperature and humidity into the these curves
$\Rightarrow$ Tsetse Suitability Index (steady-state population)
We will learn how to match these two data
\begin{align*} Y_j = \alpha + \delta TSI_j + \boldsymbol{X}'_{j}\boldsymbol{\omega} + \varepsilon_{j} \end{align*}
$Y_j$ | Precolonial outcome for ethnic group $j$ |
$TSI_j$ | Tsetse fly suitability index |
$\boldsymbol{X}_j$ | Controls (incl.\ temperature, humidity in 1871) |
(Table 1 of Alsan 2015)
1. Use the Spatial Join tool (Exercise 2)
2. Export attribute tables to Stata (Exercise 3)
3. Create grid cell polygons (Exercise 1)
In natural science data (e.g. weather data),
the unit of observations is often a grid cell.
Examples:
Grid cell polygons help merge these data with other data
We will use 3 geo-processing tools for this exercise:
As in Lecture 1, we use these tools in Model Builder.
code/lecture2.tbx/exercises1-3
Spatial resolution: 2° x 2°
Africa: roughly spans within
$\Rightarrow$ Create square polygons whose centroid is from (-20°, -40°) to (60°, 40°) at the interval of 2°
Output Feature Class: ...\Lecture2\temporary\fishnet.shp
temporary/
folder
Fishnet Origin Coordinate
X Coordinate: -21
Y Coordinate: -41
$\Rightarrow$ Bottom-left polygon's centroid will be (-20°, -40°)
Y-Axis Coordinate
X Coordinate: -21
Y Coordinate: -10 (or any value other than -41)
$\Rightarrow$ Y-Axis will be perpendicular to the horizontal X-Axis
Number of Rows: 41
Number of Columns: 41
Uncheck "Create Label Points"
Geometry Type: POLYGON
The Create Fishnet tool doesn't assign any coordination system to the output file
Ethnographic Atlas data uses geographic coordinate system
$\Rightarrow$ Assign WGS 1984
How? We already learned this in Lecture 1 Exercise 6.
Input Dataset or Feature Class: fishnet.shp
Coordinate System: GCS_WGS_1984
Now save and run the Model.
Browse the output.
Overlay the Ethnographic Atlas data
.../Lecture2/input/tribemap.shp
20th Century Reanalysis data cell polygons should cover the whole mainland Africa.
Browse the attribute table of the output (cf. Lecture 1)
We don't see centroid coordinates to match with weather data...
$\Rightarrow$ Use the Add Geometry Attributes tool
To add the centroid coordinate values in attribute table
Input Features: fishnet.shp (2)
fishnet.shp
(which is the output from Create Fishnet)
Geometry Properties: check CENTROID
Note: this tool overwrites the input file
Now save and run the Model.
Browse the output and its attribute table.
Is everything as expected?
Ethnic homelands: polygons
How to match polygons with polygons?
One way to go: match polygon centroids with polygons
Mean position of all points on the polygon boundary
$\Rightarrow$ Can be matched with a single weather data grid cell
1. Create ethnic homeland centroid point features
.../input/tribalmap.shp
)
2. Use Spatial Join to match ethnic homeland centroids with 20th Century Reanalysis grid cells
We now add these processes to the Model for Excercise 1
XY Table: ...\input\tribalmap.shp
X Field: tribelon
Y Field: tribelat
Spatial Reference: GCS_WGS_1984
Leave the other options as they are.
Input Features: tribalmap_Layer
Output Feature Class: ...\temporary\ethno_centroids.shp
Now save and run the Model.
Browse the output and its attribute table.
Is everything as expected?
Target Feature Class: ethno_centroids.shp
Join Features: fishnet.shp
Output Feature Class: ...\temporary\ethno_centroids_with_20cr2.shp
Join Operation: JOIN_ONE_TO_ONE
Check "Keep All Target Features"
Field Map of Join Features: leave as it is
Match Option: INTERSECT
Now save and run the Model.
Browse the output and its attribute table.
Is everything as expected?
Now we want to export the attribute table of ethno_centroids_with_20cr2.shp
to Stata
3 ways to export the attribute table to the format readable by Stata
Converts the attribute table into an Excel file
Then use import excel in Stata
But Excel cannot handle more than 65535 rows...
Converts the attribute table into an ASCII text file
Then use import delimited in Stata
This tool will automatically add centroid coordinates to the output ASCII file
Unlike Table To Excel, must specify which fields to export
This Stata ado directly reads a shapefile's attribute table
But it works only with polygons
1. Table To Excel
...\output\ethno_centroids_with_20cr2.xls
2. Export Feature Attribute to ASCII
...\output\ethno_centroids_with_20cr2.txt
Browse the exported tables
ethno_centroids_with_20cr2.txt
has two additional columns
Export Feature Attribute to ASCII is preferrable for
But Table To Excel is simpler to execute
Look at lec2model.tbx/exercises1-3
in the downloaded data folder solutions4exercises
Essential for replication
Convenient to repeat the same geo-processing
Use Model Builder to write a draft script
Then edit the draft script
Throughout the course, we will learn various scripting tips
1. If the model is not opened yet, right-click the model in Catalogue Window and click "Edit"
2. Click in the menu bar "Model > Export > To Python Script"
3. Save as ***.py
lecture2.py
and save it in .../Desktop/Lecture2/code/
Use Python's default editor IDLE
Or use your favourite text editor
1. Click the Windows icon at the bottom-left
2. Click ArcGIS > IDLE (Python GUI)
3. Click File > Open to read a script
1. Click the Atom icon (either on Desktop or Taskbar)
2. Click File > Add Project Folder... and choose ...Desktop\Lecture2
Lecture2/
folder
3. Click code > lecture2.py to read a script
Have your own template script for executing ArcGIS geo-processing tools
Then copy-and-paste command lines from the exported script to your template
For this lecture, I've prepared the template script (.../Lecture2/code/template4L2.py
)
Using this template script, we now learn the basics of Python scripting for ArcGIS
cf. To learn Python, visit Codecademy (free of charge)
1. Object, Property, Method
2. Try-Except statement
3. String variables
4. Print
5. Comments
An object is a set of commands
arcpy
is an object containing all the ArcGIS geoprocessing tools
The import
command reads this object (i.e. launch ArcGIS)
An object may contain smaller objects
arcpy.env
is an object contained in arcpy
Two types of commands in an object: property and method
Contains a value (numbers, strings, True/False)
To assign a new value of the property, type:
object.property = value
overwriteOutput
: whether the overwriting of output files is allowed
workspace
: working directory name
By default, overwriting the output file is not allowed.
The following command line (case-sensitive!)
arcpy.env.overwriteOutput = True
allows you to overwrite output files.
arcpy.env.workspace = "D:/temp"
Never use \
\
is used for line continuation
\
is used for directory path
Use /
or \\
\\
is a function (i.e. transform inputs into outputs)
e.g. Geo-processing tool in ArcGIS: methods of the arcpy
object
To use a method, type:
object.method(arg1, arg2, ...)
The arcpy
method names mix uppercase and lowercase letters as well as underscores
Arguments for these methods are often very long with lots of [ ] " ' ; ,
etc.
$\Rightarrow$ Better to first use Model Builder to export the draft script (until you get used to Python programming)
Visit ArcGIS for Desktop website and type the method name in the search field at the top row
1. Object, Property, Method
2. Try-Except statement
3. String variables
4. Print
5. Comments
If you simply run commands and get an error, you won't get a message on what went wrong
Try-except statement is essential to know what the error is about
Python first attempts to execute indented commands below try:
If there is no error:
Python skips all the indented commands after except:
If there is an error:
Python executes the indented commands below except:
$\Rightarrow$ Write commands for error messages after except:
1. Copy the geoprocessing commands in the exported script
2. Paste them between try:
and except:
in the template script
3. Select commands to be indented
4. Click the tab key (Atom) or Format > Indent Region (IDLE)
############################## Geoprocessing ##########################
try:
# Process: Create Fishnet
arcpy.CreateFishnet_management(fishnet_shp, "-21 -41", "-21 -10", "2", "2", "41", "41", "", "NO_LABELS", "DEFAULT", "POLYGON")
# Process: Make XY Event Layer
arcpy.MakeXYEventLayer_management(tribalmap_shp, "lon", "lat", tribalmap_Layer, "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119522E-09;0.001;0.001;IsHighPrecision", "")
# Process: Copy Features
arcpy.CopyFeatures_management(tribalmap_Layer, ethno_centroids_shp, "", "0", "0", "0")
# Process: Define Projection
arcpy.DefineProjection_management(fishnet_shp, "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
# Process: Add Geometry Attributes
arcpy.AddGeometryAttributes_management(fishnet_shp__2_, "CENTROID", "", "", "")
# Process: Spatial Join
arcpy.SpatialJoin_analysis(ethno_centroids_shp, fishnet_shp__4_, ethno_centroid_with_20cr_shp, "JOIN_ONE_TO_ONE", "KEEP_ALL", "", "INTERSECT", "", "")
# Process: Table To Excel
arcpy.TableToExcel_conversion(ethno_centroid_with_20cr_shp, ethno_centroid_with_20cr_xls, "NAME", "CODE")
# Process: Export Feature Attribute to ASCII
arcpy.ExportXYv_stats(ethno_centroid_with_20cr_shp, "CENTROID_X;CENTROID_Y;NAME;TRIBE_CODE", "SPACE", ethno_centroid_with_20cr_txt, "ADD_FIELD_NAMES")
# Return geoprocessing specific errors
except arcpy.ExecuteError:
print arcpy.GetMessages()
1. Object, Property, Method
2. Try-Except statement
3. String variables
4. Print
5. Comments
It's similar to local macro in Stata
To create a variable called file1
which contains a string data.shp
, type
file1 = "data.shp"
Recommended to assign a string variable to each data file name in the Python script
1. Copy local variables in the exported script
2. Paste them before try:
in the template script
3. Sort these names by inputs, intermediates, and outputs
4. Assign the directory path for outputs to arcpy.env.workspace
5. Delete the directory path for outputs and intermediates
Now notice
fishnet_shp__2_ = fishnet_shp
fishnet_shp__4_ = fishnet_shp__2_
This happens when we use overwriting geoprocessing tools (Define Projection, Add Geometry Attributes)
It's confusing. So:
6. Press ctrl + f (or click Edit > Replace for IDLE) to replace:
fishnet_shp__2_
with fishnet_shp
fishnet_shp__4_
with fishnet_shp
7. Delete the duplicated lines
import arcpy
### Set environment ###
# Allow the overwriting of the output files
arcpy.env.overwriteOutput = True # This command is CASE-SENSITIVE
# Set the working directory.
arcpy.env.workspace = "C:\\Users\\your_username\\Desktop\\Lecture2" # NEVER USE single backslash (\).
### Local variables ###
tribalmap_shp = "input\\tribalmap.shp"
tribalmap_Layer = "tribalmap_Layer"
ethno_centroids_shp = "temporary\\ethno_centroids.shp"
fishnet_shp = "temporary\\fishnet.shp"
ethno_centroid_with_20cr_shp = "temporary\\ethno_centroid_with_20cr.shp"
ethno_centroid_with_20cr_xls = "output\\ethno_centroid_with_20cr.xls"
ethno_centroid_with_20cr_txt = "output\\ethno_centroid_with_20cr.txt"
############################## Geoprocessing ##########################
try:
# Process: Create Fishnet
1. Object, Property, Method
2. Try-Except statement
3. String variables
4. Print
5. Comments
It's a Python's equivalent of display
in Stata
To show a message on the IDLE screen while running the script, type:
print "message"
print arcpy.GetMessages()
For each geo-processing step, add the print
command to explain what's being processsed.
Also add the print commands at the beginning of the script
############################## Preamble #################################
print "Launch ArcGIS"
import arcpy
### Set environment ###
print "Enable the overwriting"
arcpy.env.overwriteOutput = True # This command is CASE-SENSITIVE
print "Set the working directory"
arcpy.env.workspace = "C:\\Users\\Masayuki Kudamatsu\\Desktop\\Lecture2" # NEVER USE single backslash (\).
### Local variables ###
print "Input file names being set"
tribalmap_shp = "input\\tribalmap.shp"
print "Intermediate file names being set"
tribalmap_Layer = "tribalmap_Layer"
ethno_centroids_shp = "temporary\\ethno_centroids.shp"
fishnet_shp = "temporary\\fishnet.shp"
ethno_centroid_with_20cr_shp = "temporary\\ethno_centroid_with_20cr.shp"
print "Output files being set"
ethno_centroid_with_20cr_xls = "output\\ethno_centroid_with_20cr.xls"
ethno_centroid_with_20cr_txt = "output\\ethno_centroid_with_20cr.txt"
############################## Geoprocessing ##########################
try:
print "Create 20th Century Reanalysis data cell polygons step 1/3"
arcpy.CreateFishnet_management(fishnet_shp, "-21 -41", "-21 -10", "2", "2", "41", "41", "", "NO_LABELS", "DEFAULT", "POLYGON")
print "Create 20th Century Reanalysis data cell polygons step 2/3"
arcpy.DefineProjection_management(fishnet_shp, "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
print "Create 20th Century Reanalysis data cell polygons step 3/3"
arcpy.AddGeometryAttributes_management(fishnet_shp, "CENTROID", "", "", "")
print "Create ethnic group centroids step 1/2"
arcpy.MakeXYEventLayer_management(tribalmap_shp, "lon", "lat", tribalmap_Layer, "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119522E-09;0.001;0.001;IsHighPrecision", "")
print "Create ethnic group centroids step 2/2"
arcpy.CopyFeatures_management(tribalmap_Layer, ethno_centroids_shp, "", "0", "0", "0")
print "Match Ethnographic Atlas ethnic groups with 20th Century Reanalysis data cells"
arcpy.SpatialJoin_analysis(ethno_centroids_shp, fishnet_shp, ethno_centroid_with_20cr_shp, "JOIN_ONE_TO_ONE", "KEEP_ALL", "", "INTERSECT", "", "")
print "Export output data as Excel sheet"
arcpy.TableToExcel_conversion(ethno_centroid_with_20cr_shp, ethno_centroid_with_20cr_xls, "NAME", "CODE")
print "Export output data as csv file"
arcpy.ExportXYv_stats(ethno_centroid_with_20cr_shp, "CENTROID_X;CENTROID_Y;NAME;TRIBE_CODE", "SPACE", ethno_centroid_with_20cr_txt, "ADD_FIELD_NAMES")
# Return geoprocessing specific errors
except arcpy.ExecuteError:
1. Object, Property, Method
2. Try-Except statement
3. String variables
4. Print
5. Comments
To insert comments, type #
at the beginning of a comment
Useful for others to understand (and for you to remember) what each command does
Also useful if you want to skip some commands for the time being
Suppose you're not sure if this script will properly create ethnic group centroids.
So you want to temporarily skip the other parts of the script
1. Select command lines to be skipped
except:
part
2. Press ctrl + /
$\Rightarrow$ #
will appear at the beginning of selected commands
############################## Geoprocessing ##########################
try:
# print "Create 20th Century Reanalysis data cell polygons step 1/3"
# arcpy.CreateFishnet_management(fishnet_shp, "-21 -41", "-21 -10", "2", "2", "41", "41", "", "NO_LABELS", "DEFAULT", "POLYGON")
#
# print "Create 20th Century Reanalysis data cell polygons step 2/3"
# arcpy.DefineProjection_management(fishnet_shp, "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
#
# print "Create 20th Century Reanalysis data cell polygons step 3/3"
# arcpy.AddGeometryAttributes_management(fishnet_shp, "CENTROID", "", "", "")
arcpy.MakeXYEventLayer_management(tribalmap_shp, "lon", "lat", tribalmap_Layer, "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119522E-09;0.001;0.001;IsHighPrecision", "")
print "Create ethnic group centroids step 1/2"
print "Create ethnic group centroids step 2/2"
arcpy.CopyFeatures_management(tribalmap_Layer, ethno_centroids_shp, "", "0", "0", "0")
# print "Match Ethnographic Atlas ethnic groups with 20th Century Reanalysis data cells"
# arcpy.SpatialJoin_analysis(ethno_centroids_shp, fishnet_shp, ethno_centroid_with_20cr_shp, "JOIN_ONE_TO_ONE", "KEEP_ALL", "", "INTERSECT", "", "")
#
# print "Export output data as Excel sheet"
# arcpy.TableToExcel_conversion(ethno_centroid_with_20cr_shp, ethno_centroid_with_20cr_xls, "NAME", "CODE")
#
# print "Export output data as csv file"
# arcpy.ExportXYv_stats(ethno_centroid_with_20cr_shp, "CENTROID_X;CENTROID_Y;NAME;TRIBE_CODE", "SPACE", ethno_centroid_with_20cr_txt, "ADD_FIELD_NAMES")
# Return geoprocessing specific errors
except arcpy.ExecuteError:
3. Save and run the script
4. Cross your fingers :-)
5. Check if ethno_centroids.shp
is created properly.
6. Now uncomment those skipped commands
7. Save and run the entire script
8. Cross your fingers :-)
9. Check if everything is as expected
Most likely, you get an error message:
ERROR 000464: Cannot get exclusive schema lock. Either being edited or in use by another application.
This error occurs when the script is trying to overwrite the output that's been shown in ArcMap
Before you run the Python script, always close all the output files in ArcMap
Sometimes, the output file is not what you intended.
e.g., New fields in the attribute table are all zero
Revising the script may not solve the issue / You cannot find anything wrong with the script
If so, reboot ArcMap/IDLE and even Windows. This often solves the issue.
Still unsolved, ask Google.
Last resort: ask at Stack Exchange or ArcGIS Forum.
Look at
solutions4exercises/lec2script.py
in the downloaded data folder for Lecture 2
Do you remember which geo-processing tools you used for each of these tasks?