需求:
基于属性表不同字段(例如:砷、镉、铜、铅、汞、镍)和属性表中分层字段中的不同分层(第1层、第2层.........、第6层)的要素记录进行空间插值。
例如:在空间插值之前,将砷字段基于分层字段细分为6个单独的砷图层,即第1层的砷字段对应的要素数据、第2层的砷字段对应的要素数据、第3层的砷字段对应的要素数据、第4层的砷字段对应的要素数据、第5层的砷字段对应的要素数据、第6层的砷字段对应的要素数据;其它字段(镉、铜、铅、汞、镍)均遵循此数据选择方式,当然了这些过程都是在代码中完成的。
代码(pycharm中编译):
# -*-coding:utf-8-*-
import arcpy
from arcpy import env
# Set workspace
env.coincidentPoints = "MAX"
env.workspace = "D:/test/data.mdb"
# 定义污染指标
pollutionParameter = ['砷', '镉', '铜','铅', '汞', '镍']
# 遍历污染指标
for indexwrzb in pollutionParameter:
count = 1
print("污染指标: "+indexwrzb)
# 遍历每个污染指标分层类型,例如:第1层、第2层、第3层、第4层、第5层、第6层
for layerId in range(6):
in_features = "土壤污染数据全区"
out_feature_class = "土壤污染%s数据全区第%s层" %(indexwrzb, (layerId + 1))
currentWrzb = '[%s]' % indexwrzb
where_clause = '%s>0 and [分层]="第%s层" ' % (currentWrzb, (layerId + 1))
# Execute Select
arcpy.Select_analysis(in_features, out_feature_class, where_clause)
inPointFeatures = out_feature_class
#判断图层要素记录的条数,如果少于3条则提示不能进行空间差值操作;
#如果多余或者等于3条,那么可以进行空间插值操作
featureRecord = arcpy.GetCount_management(inPointFeatures)
featureRecordCount = int(featureRecord.getOutput(0))
print (str(indexwrzb)+"的图层"+str(count)+"的要素记录条数为: "+str(featureRecordCount))
if(featureRecordCount>2):
#将当前遍历出的污染指标赋值给zField,用此作为空间插值的Z值
zField = indexwrzb
outLayer = "outLPI"+str(indexwrzb)+str(count)
outRaster = "D:/test/" + str(indexwrzb) + "数据第" + str(layerId+1) + "层.tif"
cellSize = 3
power = 2
# Set variables for search neighborhood
majSemiaxis = 250
minSemiaxis = 250
angle = 0
maxNeighbors = 15
minNeighbors = 10
sectorType = "ONE_SECTOR"
searchNeighbourhood = arcpy.SearchNeighborhoodStandard(majSemiaxis, minSemiaxis,
angle, maxNeighbors,
minNeighbors, sectorType)
# Check out the ArcGIS Geostatistical Analyst extension license
arcpy.CheckOutExtension("GeoStats")
# Execute IDW
arcpy.IDW_ga(inPointFeatures, zField, outLayer, outRaster, cellSize,power, searchNeighbourhood,"")
count = count + 1
else:
print (str(indexwrzb)+"图层"+str(count)+"的要素记录少于3个,不能进行空间差值......")
continue
print ("ok.......")
注意事项:
(1)# -*-coding:utf-8-*- 如果脚本中有中文字符,则需要在脚本第一行加上utf-8的编码
(2)env.coincidentPoints = "MAX",多个点重合在一起,用其中值最大的点来进行插值。
(3)len(pollutionParameter) 获取字段wrzb的长度,此示例中wrzb的长度是6
(4)range(6) 输出记录[0,1,2,3,4,5],从0开始到6,参考链接:http://www.runoob.com/python/python-func-range.html
(5)where_clause = '[砷]>0 and [分层]="第%s层" ' % (indexwrzb + 1) 条件语句,获取砷字段下属性值大于0且分层字段属性值为第%s层的要素记录,其中%s是一个变量值。
(6)在代码中判断图层要素记录的条数,如果少于3条则提示不能进行空间差值操作; 如果多余或者等于3条,那么可以进行空间插值操作。
(7)获取要素记录的数量,用下述代码实现:
featureRecord = arcpy.GetCount_management(inPointFeatures)
featureRecordCount = int(featureRecord.getOutput(0))
(8)arcpy.CheckOutExtension("GeoStats")是用于检查地统计分析扩展的许可,如果是在独立的python编译器中需要检查扩展许可,如果在arcgis自带的python window中则不要检查扩展许可。
日志输出:
污染指标: 砷
砷的图层1的要素记录条数为: 455
砷的图层2的要素记录条数为: 683
砷的图层3的要素记录条数为: 272
砷的图层4的要素记录条数为: 16
砷的图层5的要素记录条数为: 6
砷的图层6的要素记录条数为: 1
砷图层6的要素记录少于3个,不能进行空间差值......
污染指标: 镉
镉的图层1的要素记录条数为: 455
镉的图层2的要素记录条数为: 683
镉的图层3的要素记录条数为: 273
镉的图层4的要素记录条数为: 16
镉的图层5的要素记录条数为: 6
镉的图层6的要素记录条数为: 1
镉图层6的要素记录少于3个,不能进行空间差值......
污染指标: 铜
铜的图层1的要素记录条数为: 455
铜的图层2的要素记录条数为: 683
铜的图层3的要素记录条数为: 273
铜的图层4的要素记录条数为: 16
铜的图层5的要素记录条数为: 6
铜的图层6的要素记录条数为: 1
铜图层6的要素记录少于3个,不能进行空间差值......
污染指标: 铅
铅的图层1的要素记录条数为: 455
铅的图层2的要素记录条数为: 683
铅的图层3的要素记录条数为: 273
铅的图层4的要素记录条数为: 16
铅的图层5的要素记录条数为: 6
铅的图层6的要素记录条数为: 1
铅图层6的要素记录少于3个,不能进行空间差值......
污染指标: 汞
汞的图层1的要素记录条数为: 455
汞的图层2的要素记录条数为: 682
汞的图层3的要素记录条数为: 273
汞的图层4的要素记录条数为: 16
汞的图层5的要素记录条数为: 6
汞的图层6的要素记录条数为: 1
汞图层6的要素记录少于3个,不能进行空间差值......
污染指标: 镍
镍的图层1的要素记录条数为: 455
镍的图层2的要素记录条数为: 683
镍的图层3的要素记录条数为: 273
镍的图层4的要素记录条数为: 16
镍的图层5的要素记录条数为: 6
镍的图层6的要素记录条数为: 1
镍图层6的要素记录少于3个,不能进行空间差值......
ok.......
输出结果:
参考资料:
http://www.runoob.com/python/python-func-range.html range函数的用法
https://blog.csdn.net/u010608964/article/details/87428910 利用arcpy遍历shp文件并获取要素记录数