利用Arcpy和Geopandas进行多规比对
在ArcGIS中可以很容易进行土规和城规的比对，这里，利用Arcpy和Geopandas来进行多规比对。
# encoding: utf-8
# @Author : hanqi
数据预处理
导入相关模块
## 导入相关模块
import arcpy
from arcpy import env
import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
防止乱码
plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.sans-serif'] = ['SimHei']#替换sans-serif字体为黑体
plt.rcParams['axes.unicode_minus'] = False # 解决坐标轴负数的负号显示问题
设置arcpy工作环境和路径
env.workspace = r"F:\ArcGIS\多规对比"
path = "F:\ArcGIS\多规对比\图斑对比.gdb"
tg_path = "F:\ArcGIS\多规对比\图斑对比.gdb\土规用地"
cg_path = "F:\ArcGIS\多规对比\图斑对比.gdb\城规用地"
土规和城规用地分类
读取数据
land = gpd.GeoDataFrame.from_file(path,layer='土规用地')
city = gpd.GeoDataFrame.from_file(path,layer='城规用地')
land.head()
city.head()
arcpy数据处理
## 删除字段
## 因为我做过一遍，所以先删除
arcpy.DeleteField_management(tg_path, "tg_yongdi")
arcpy.DeleteField_management(cg_path, "cg_yongdi")
## 添加字段
fieldName1 = "tg_yongdi"
fieldLength = 30
# Execute AddField twice for two new fields
arcpy.AddField_management(tg_path, fieldName1, "TEXT", field_length=fieldLength)
fieldName2 = "cg_yongdi"
fieldLength = 30
arcpy.AddField_management(cg_path, fieldName2, "TEXT", field_length=fieldLength)
## 字段计算器
## 城规字段计算器
## 只有一个水域不是建设用地
# Set local variables
inTable = cg_path
fieldName = "cg_yongdi"
expression = "getClass((!cg_fenlei!))"
codeblock = """
def getClass(cg_fenlei):
if cg_fenlei in "E11":
return "城规_非建设用地"
else:
return "城规_建设用地" """
# Execute CalculateField
arcpy.CalculateField_management(inTable, fieldName, expression, "PYTHON_9.3", codeblock)
## 土规字段计算器
# Set local variables
inTable = tg_path
fieldName = "tg_yongdi"
expression = "getClass((!Layer!))"
codeblock = """
def getClass(Layer):
if Layer in ("TG-一般农田", "TG-基本农田", "TG-林地", "TG-水系", "TG-林地"):
return "土规_非建设用地"
if Layer is None:
return None
else:
return "土规_建设用地" """
# Execute CalculateField
arcpy.CalculateField_management(inTable, fieldName, expression, "PYTHON_9.3", codeblock)
结果就是上图，土规和城规的用地分为两类了，这里不做演示。
多规比对
## 联合工具
inFeatures = [cg_path, tg_path]
outFeatures = outFeatures = path + "\\多规比对"
arcpy.Union_analysis (inFeatures, outFeatures, "ALL")
dg_path = "F:\ArcGIS\多归对比\图斑对比.gdb\多规比对"
dg = gpd.GeoDataFrame.from_file(path, layer="多规比对")
dg.head()
# Set local variables
## 多规字段计算器
inTable = dg_path
fieldName = "cg_yongdi"
expression = "getClass3((!cg_yongdi!))"
codeblock = """
def getClass3(cg_yongdi):
if cg_yongdi is None:
return "城规_非建设用地"
else:
return cg_yongdi
"""
# Execute CalculateField
arcpy.CalculateField_management(inTable, fieldName, expression, "PYTHON_9.3", codeblock)
## 删除字段
arcpy.DeleteField_management(dg_path, ["FID_城规用地","Layer","cg_fenlei","FID_土规用地"])
## 添加合并字段
fieldName3 = "dg_hebing"
fieldLength = 30
# Execute AddField twice for two new fields
arcpy.AddField_management(dg_path, fieldName3, "TEXT", field_length=fieldLength)
## 城规字段合并
# Set local variables
inTable = dg_path
fieldName = "dg_hebing"
expression = "!tg_yongdi!+'_'+!cg_yongdi!"
# Execute CalculateField
arcpy.CalculateField_management(inTable, fieldName, expression, "PYTHON_9.3")
dg = gpd.GeoDataFrame.from_file(path, layer="多规比对")
dg_path = path+'\\多规比对'
dg
冲突对比
## 添加索引
rec = 0
inTable = dg_path
fieldName = "Identify"
expression = "autoIncrement()"
codeblock = """
def autoIncrement():
global rec
pStart = 1
pInterval = 1
if (rec == 0):
rec = pStart
else:
rec = rec + pInterval
return rec """
# Execute AddField
arcpy.AddField_management(inTable, fieldName, "LONG")
# Execute CalculateField
arcpy.CalculateField_management(inTable, fieldName, expression, "PYTHON_9.3", codeblock)
dg = gpd.GeoDataFrame.from_file(path, layer="多规比对")
dg.head()
## 按属性分割
# Set local variables
in_feature_class = dg_path
target_workspace = path
fields = ['dg_hebing']
arcpy.SplitByAttributes_analysis(in_feature_class, target_workspace, fields)
cgfjs_tgjs = gpd.GeoDataFrame.from_file(path, layer="土规_非建设用地_城规_建设用地")
cgfjs_tgjs.plot()
cgjs_tgfjs = gpd.GeoDataFrame.from_file(path, layer="土规_建设用地_城规_非建设用地")
cgjs_tgfjs.plot()
区划调整
土规非建设用地_城规建设用地调整
## 定义路径
tgfjs_cgjs_path = path+"\\土规_非建设用地_城规_建设用地"
# Set local variables
## 土地面积大于10000依城规
inTable = tgfjs_cgjs_path
fieldName = "gz_1"
expression = "getClass((!Shape_Area!))"
codeblock = """
def getClass(Shape_Area):
if Shape_Area>10000:
return "依城规用地"
else:
return "依土规用地" """
# Execute CalculateField
arcpy.CalculateField_management(inTable, fieldName, expression, "PYTHON_9.3", codeblock)
土规建设用地_城规非建设用地调整
## 限制1000以内不隐藏
pd.set_option('max_columns',1000)
pd.set_option('max_row',1000)
tgjs_cgfjs_path = path+"\\土规_建设用地_城规_非建设用地"
tgjs_cgfjs = gpd.GeoDataFrame.from_file(path, layer="土规_建设用地_城规_非建设用地")
tgjs_cgfjs.sort_values('Shape_Area',ascending=False) ##观察面积
# Set local variables
## 土地面积大于10000以土规
inTable = tgjs_cgfjs
fieldName = "gz_2"
expression = "getClass5((!Shape_Area!))"
codeblock = """
def getClass5(Shape_Area):
if Shape_Area>10000:
return "依土规用地"
else:
return "依城规用地" """
# Execute CalculateField
arcpy.CalculateField_management(inTable, fieldName, expression, "PYTHON_9.3", codeblock)
tgjs_cgfjs = gpd.GeoDataFrame.from_file(path, layer="土规_建设用地_城规_非建设用地")
tgjs_cgfjs.plot(column='gz_2',legend=True)
数据连接
layerName = dg_path
joinTable = cgjs_tgfjs_path
joinField = "Identify"
# Join the feature layer to a table
arcpy.AddJoin_management(layerName, joinField, joinTable, joinField)
layerName = '多规比对_Layer1'
joinTable = cgfjs_tgjs_path
joinField = "Identify"
# Join the feature layer to a table
arcpy.AddJoin_management(layerName, joinField, joinTable, joinField)
## 导出数据
# 输入
in_features = ['多规比对_Layer1']
# 输出
out_location = path
# 执行 FeatureClassToGeodatabase
arcpy.FeatureClassToGeodatabase_conversion(in_features, out_location)
dg_3 = gpd.GeoDataFrame.from_file(path, layer="多规比对_Layer3")
dg_3.head()
## 重命名
# Set local variables
in_data = "多规比对_Layer1"
out_data = "多规比对_最终"
# Execute Rename
arcpy.Rename_management(in_data, out_data)
# Set local variables
inTable = path+"\\多规对比_最终"
fieldName = "gz_all"
expression = "getClass5(!土规_非建设用地_城规_建设用地_gz_1! , !土规_建设用地_城规_非建设用地_gz_2!)" ## 不能用别名 枣庄人流医院哪家好 http://mobile.0632-3679999.com/
codeblock = """
def getClass6(gz_1,gz_2):
if gz_1 is not None:
return gz_1
if gz_2 is not None:
return gz_2
else:
return "无冲突" """
arcpy.AddField_management(inTable, fieldName, "TEXT", 30)
# Execute CalculateField
arcpy.CalculateField_management(inTable, fieldName, expression, "PYTHON_9.3", codeblock)
dg_zz = gpd.GeoDataFrame.from_file(path, layer="多规比对_最终")
dg_zz.plot(column='gz_all',legend=True)
用地判定
就是根据冲突比对来返回想对呀的地块用地，这里我用字段计算pyhton做不出来，用vb可以，但是可以用pandas来处理。
vb脚本
if [gz_all] = "依土规用地" then
value = [多规比对_tg_yongdi]
if [gz_all] = "依成规用地" then
value = [多规比对_cg_yongdi]
if [gz_all] = "无冲突" then
value = [多规比对_cg_yongdi]
end if
value
pandas
## 最终用地
for i in dg_zz.index:
if dg_zz.at[i,'gz_all']=="无冲突":
dg_zz["最终用地"].at[i] = dg_zz["多规比对_tg_yongdi"].at[i]
if dg_zz.at[i,'gz_all'] == "依城规用地":
dg_zz["最终用地"].at[i] = dg_zz["多规比对_cg_yongdi"].at[i]
else:
dg_zz["最终用地"].at[i] = dg_zz["多规比对_tg_yongdi"].at[i]
dg_zz.head()
## 分列
dg_zz["用地"] = dg_zz['最终用地'].str.split('_',expand=True)[1]
dg_zz.head()
#dg_zz["用地"] = s.split('_')[1]
结果，这里我进行过很多次计算了
## 写入excel文件
dg_zz.to_excel("成果.xlsx",index=None,encoding="utf8")
然后连接回shp文件就行了，有些ArcGIS不支持xlsx，保存的时候可以注意点，我的可以
字段计算器分列
## 分割字段
# Set local variables
inTable = path + "\\多规比对_最终"
fieldName = "yongdi"
expression = "getClass6(!zz_yongdi!)" ## 不能用别名
codeblock = """
def getClass6(zz_yongdi):
value = zz_yongdi.split('_')[1]
return value """
arcpy.AddField_management(inTable, fieldName, "TEXT", 30)
# Execute CalculateField
arcpy.CalculateField_management(inTable, fieldName, expression, "PYTHON_9.3", codeblock)
dg_zz = gpd.GeoDataFrame.from_file(path, layer="多规比对_最终")
dg_zz.plot(column='zz_yongdi',legend=True)
成果出图
fig, ax = plt.subplots(figsize=(10, 10))
ax = dg_zz.plot(ax=ax,column='yongdi',cmap='Set1',legend=True,legend_kwds={
'loc': 'lower left',
'title': '图例',
'shadow': True})
plt.suptitle('土地利用建设分布图', fontsize=24) # 添加最高级别标题
plt.tight_layout(pad=4.5) # 调整不同标题之间间距
plt.grid(True, alpha=0.3)
fig.savefig('土地利用建设分布图.png', dpi=300)
fig, ax = plt.subplots(figsize=(10, 10))
ax = dg_zz.plot(ax=ax,column='gz_all',cmap='tab10',legend=True,legend_kwds={
'loc': 'lower left',
'title': '图例',
'shadow': True})
plt.suptitle('规划冲突调整图', fontsize=24) # 添加最高级别标题
plt.tight_layout(pad=4.5) # 调整不同标题之间间距
plt.grid(True, alpha=0.3)
fig.savefig('规划冲突调整图.png', dpi=300)