利用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!)"  ## 不能用別名
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)