ArcGIS
ArcGIS
因为某种神的号召，我要用 Python 做这次实习。调用 ArcGIS 的 ArcPy 之前的准备工作
解压数据
打开地图
打开 Python Shell
然后就是把相应代码粘贴到 Python Shell 中执行（回车即可）
代码放在各 DEMO 文件夹的 x-py 文件夹内，有些我试验的 Python 脚本放在 more-scripts 文件夹
path_base = 'X://DEMO/';
path_gdb = path_base + 'demo.gdb/';
path_data = path_base + 'Spatial/';
# import modules
import sys, os, arcpy;
import arcpy.mapping as mapping;
from arcpy.sa import *
from arcpy.mapping import *# set workspace
arcpy.env.workspace = path_gdb;# source data
path_elevation = path_data + 'elevation'
path_landuse = path_data + 'landuse'
path_school_site = path_data + 'School_site.shp'
path_destination = path_data + 'Stowe.gdb/destination'
path_rec_sites = path_data + 'Stowe.gdb/rec_sites'
path_roads = path_data + 'Stowe.gdb/roads'
path_schools = path_data + 'Stowe.gdb/schools'# get document
mxd = mapping.MapDocument("CURRENT");
df = mapping.ListDataFrames(mxd)[0]; # Get First Data frame, a.k.a. "Layers" (empty by far)# Add Layers
layer = mapping.Layer(path_elevation);
mapping.AddLayer(df, layer, "AUTO_ARRANGE")
layer = mapping.Layer(path_landuse);
mapping.AddLayer(df, layer, "AUTO_ARRANGE")
layer = mapping.Layer(path_school_site);
mapping.AddLayer(df, layer, "AUTO_ARRANGE")
layer = mapping.Layer(path_destination);
mapping.AddLayer(df, layer, "AUTO_ARRANGE")
layer = mapping.Layer(path_rec_sites);
mapping.AddLayer(df, layer, "AUTO_ARRANGE")
layer = mapping.Layer(path_roads);
mapping.AddLayer(df, layer, "AUTO_ARRANGE")
layer = mapping.Layer(path_schools);
mapping.AddLayer(df, layer, "AUTO_ARRANGE")
# Adding Layers Done, Fantastic~# set relative path
mxd.relativePaths = True;for df in mapping.ListDataFrames(mxd):
if (df.name == 'Layers'): # that's for sure
layers = mapping.ListLayers(mxd, 'elevation', df);
# set extent to whole map, that's what we want
arcpy.env.extent = layers[0].getExtent();arcpy.CheckOutExtension("Spatial");如果你没打开的话，大概会看到个这个
spatial extension error
根据地形生成阴影、坡度栅格图，根据欧式距离生成 DisTanceRec 和 DistanceSch 图
# Execute HillShade
outHillShade = arcpy.sa.Hillshade("elevation", 180, 75, "SHADOWS", 1); # done
outHillShade.save(path_gdb + 'HillShade');
# fine
Slope = arcpy.sa.Slope("elevation", "DEGREE");
Slope.save(path_gdb + 'Slope');
# Distance to Rec Site
DistanceRec = arcpy.sa.EucDistance("rec_sites", cell_size = 30);
DistanceRec.save(path_gdb + 'DistanceRec');
# Distance to School
DistanceSch = arcpy.sa.EucDistance("schools", cell_size = 30);
DistanceSch.save(path_gdb + 'DistanceSch');降低距离的分类，重分类
# this is it, Reclassify
ReclassSlope = arcpy.sa.Reclassify("Slope", "Value",
RemapRange([[0,8,1],[8,15,2],[15,23,3],[23,31,4],[31,39,5],
[39,46,6],[46,54,7],[54,62,8],[62,69,9],[69,78,10]]), "NODATA"); # this is what we want~
ReclassSlope.save(path_gdb + 'ReclassSlope');
ReclassRec = arcpy.sa.Reclassify("DistanceRec", "Value",
RemapRange([[0,1349,1],[1349,2698,2],[2698,4046,3],[4046,5395,4],[5395,6744,5],
[6744,8093,6],[8093,9441,7],[12139,13488,8],[10790,12139,9],[12139,13488,10]]), "NODATA");
ReclassRec.save(path_gdb + 'ReclassRec');手工太累，用一个函数生成重分类的 RemapRange 矩阵，再重分类
# seem we need a function to help us out
def partition (beg, end, num):
pace = ( end - beg ) / num + 1;
return [[beg+pace*i, beg+pace*i+pace] for i in xrange(num)]
range = partition (0, 16928, 10);
ReclassSch = arcpy.sa.Reclassify("DistanceSch", "Value",
RemapRange(range), "NODATA");
ReclassSch.save(path_gdb + 'ReclassSch');按照实习指导书的要求，输出
with arcpy.da.SearchCursor("landuse", ("LANDUSE")) as cursor:
for row in sorted(cursor):
print row[0]
粘贴并运行
运行框的输出信息
桌面右下角的输出信息
代码输出结果
最终输出结果
DEMO2, test1
打开 mxd 空地图文件，用相应脚本加载数据
通过属性选定数据，用 Arcpy 中的 SelectLayerByAttribute_management。
arcpy.SelectLayerByAttribute_management("Streets", "NEW_SELECTION", '"STR" = \'I40\'');ArcPy 里同一个函数有两种调用方法，一种是这里的方式：功能（SelectLayerByAttribute）+ 分类（management）；另一种是按照 ArcToolBox 的结构分层次，如
arcpy.sa.Slope()，其中 sa 指 Spatial Analysis，是相应 ArcToolBox 的别名（alias），可以通过右键相应 ToolBox 查看。后一种的层次性和条理性更清晰，但前一种，更直接得面向功能，即你想干什么。根据具体情况，选择适合的方式。
arcpy.SelectLayerByLocation_management("Stations", "WITHIN_A_DISTANCE", "Streets", "1000 Feet", "NEW_SELECTION");
count = arcpy.GetCount_management("Stations");
# Ques: How many sites, and their name?
# Ans:
print 'There are ' + str(count) + ' Stations are ... within 140 feets from 140 road'; # ===> 2, ...
with arcpy.da.SearchCursor("Stations", ("NAME")) as cursor:
for row in sorted(cursor):
print row[0] # have a better looking than `print row'arcpy.SelectLayerByLocation_management("Business", "WITHIN_A_DISTANCE", "Stations", "1320 Feet", "NEW_SELECTION");
count = arcpy.GetCount_management( "Business"); # ==> 19
print "有 " + str(count) + " 个商业站点位于已选中的加油站 1320 英尺以内"
arcpy.SelectLayerByAttribute_management("Zoning", "NEW_SELECTION", '"DESCRIPTIO" = \'DRAIN\'');
arcpy.SelectLayerByLocation_management("Zoning", "BOUNDARY_TOUCHES", "Zoning"); # ==> 37
juris= [row[0] for row in arcpy.da.SearchCursor("Zoning", ("JURISDICTI"))];
JURIS=set(juris);
# >>> print JURIS
# set([u'COUNTY', u'CITY'])
acres= [row[0] for row in arcpy.da.SearchCursor("Zoning", ("ACRES"))];
key-val= [(row[0], row[1]) for row in arcpy.da.SearchCursor("Zoning", ("JURISDICTI", "ACRES"))]
# how to use dict?
DOME2, other tests
其他的几个处理都大同小异，这里只看关键的几个代码
根据属性（非空间的）和位置（空间的）来选择 Feature（表中的行，数据项）
arcpy.SelectLayerByAttribute_management("Centract", "NEW_SELECTION", '"POP_90" < 5786.054054');
# arcpy.SelectLayerByLocation_management("Centract", "COMPLETELY_CONTAINS", "MidSchol");
arcpy.SelectLayerByLocation_management("Centract", "COMPLETELY_CONTAINS", "MidSchol", selection_type="SUBSET_SELECTION");
arcpy.SelectLayerByAttribute_management("Counties", "NEW_SELECTION", '"Name" = \'Cobb County\'');
arcpy.SelectLayerByLocation_management("Centract", "COMPLETELY_WITHIN", "Counties", selection_type="SUBSET_SELECTION");
count = arcpy.GetCount_management( "Centract");
print count
percap= [row[0] for row in arcpy.da.SearchCursor("Centract", ("PER_CAPINC"))];
print percap # ===> [19867.69921875, 12793.2998046875, 17364.80078125]输出效果
# arcpy.analysis.Buffer("vernal_pools", "vernal_pools_Buffer.shp", "2000 Feets", "FULL", "ROUND", "ALL");
# Not Feets, it's Feet
arcpy.analysis.Buffer("vernal_pools", "vernal_pools_Buffer.shp", "2000 Feet", "FULL", "ROUND", "ALL");
arcpy.SelectLayerByLocation_management("parcels", "INTERSECT", "vernal_pools_Buffer");
arcpy.MakeFeatureLayer_management("parcels", "parcels_out") # 靠，直接导出，就是导出 selected features，查了我半天需要注意的是
2000 Feet不能写成2000 Feets。有时候 ArcGIS 会无故卡死，不知道什么情况
# Processing
arcpy.analysis.Buffer("Roads", "Roads_BUFFER", "500 Meter", "FULL", "ROUND", "ALL");
arcpy.analysis.Buffer("Water", "Water_BUFFER", "500 Meter", "FULL", "ROUND", "ALL");
arcpy.analysis.Union(("Water_BUFFER", "Roads_BUFFER"), "water_road", "ALL", gaps="GAPS")
