空间分析实习报告

实习目的

ArcGIS

ArcGIS

实习过程

过程综述

  1. 实习按照实习指导书进行,分为栅格、矢量、三维等几大版块;
  2. 按照实习指导书要求,回答相应问题,并提交几份地图;
  3. 最后按要求撰写实习报告。

一点说明

因为某种神的号召,我要用 Python 做这次实习。调用 ArcGIS 的 ArcPy 之前的准备工作

解压数据

解压数据

打开地图

打开地图

打开 Python Shell

打开 Python Shell

然后就是把相应代码粘贴到 Python Shell 中执行(回车即可)

代码放在各 DEMO 文件夹的 x-py 文件夹内,有些我试验的 Python 脚本放在 more-scripts 文件夹

栅格部分(一个完整的例子, DEMO)

  1. 设置工作路径,加入必要的 Python Module
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 *
  1. 设置工作区
# set workspace
arcpy.env.workspace = path_gdb;
  1. 指定数据来源
# 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'
  1. 获取当前地图的第一个 dataFrame,以添加数据
# get document
mxd = mapping.MapDocument("CURRENT");
df = mapping.ListDataFrames(mxd)[0]; # Get First Data frame, a.k.a. "Layers" (empty by far)
  1. 添加数据
# 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~
  1. 设置数据的路径为相对路径便于分享拷贝后使用
# set relative path
mxd.relativePaths = True;
  1. 设置处理的边界(包围盒?),这样处理就是对全图,而不是局部,便于
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();
  1. 确定下,你的空间分析组件有没有打上勾
arcpy.CheckOutExtension("Spatial");

如果你没打开的话,大概会看到个这个

spatial extension error

spatial extension error

  1. 开始处理

根据地形生成阴影、坡度栅格图,根据欧式距离生成 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]
  1. 处理结果截图
粘贴并运行

粘贴并运行

运行框的输出信息

运行框的输出信息

桌面右下角的输出信息

桌面右下角的输出信息

代码输出结果

代码输出结果

最终输出结果

最终输出结果

矢量部分

DEMO2, test1

  1. 打开 mxd 空地图文件,用相应脚本加载数据

  2. 通过属性选定数据,用 Arcpy 中的 SelectLayerByAttribute_management

arcpy.SelectLayerByAttribute_management("Streets", "NEW_SELECTION", '"STR" = \'I40\'');
注:

ArcPy 里同一个函数有两种调用方法,一种是这里的方式:功能(SelectLayerByAttribute)+ 分类(management);另一种是按照 ArcToolBox 的结构分层次,如 arcpy.sa.Slope(),其中 sa 指 Spatial Analysis,是相应 ArcToolBox 的别名(alias),可以通过右键相应 ToolBox 查看。后一种的层次性和条理性更清晰,但前一种,更直接得面向功能,即你想干什么。根据具体情况,选择适合的方式。

  1. 选定数据,统计数量,按照指导书要求回答问题
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'
  1. 一些别的处理
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?
  1. 输出

DOME2, other tests

  1. 其他的几个处理都大同小异,这里只看关键的几个代码

  2. 根据属性(非空间的)和位置(空间的)来选择 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]

输出效果

  1. 缓冲区操作,最后导出到 shp 文件
# 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 会无故卡死,不知道什么情况 arcgis-died

  1. 缓冲区操作,并合并
# 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")

其他

  1. 其他的部分用 ArcPy 没法搞定,首先就是无法打开网络分析的 “Streets_ND” 数据,从我查到的资料表明 ArcPy 做不了。
  2. 可能是我的 ArcGIS 权限不够的问题,我没法新建 ToolBox。
  3. 手工做了一点其他的内容,实在太烂,不敢多发

实习感悟

  1. 通过实习,巩固了 ArcGIS 的使用,巩固了 Python 编程;
  2. 锻炼了查资料,查文档的能力,英文搜索投放关键字的水平也大大提高了;
  3. 熟悉了 Windows 下 notepad++ 的使用, Linux 下 Bash 脚本处理文档的小技巧;
  4. 这是一个快乐而充实的五一国际劳动节;
  5. 还是 Linux 好用,写报告多方便多了,see:

「Load Disqus | 加载评论」