　　利用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)

利用Arcpy和Geopandas進行多規比對

相關文章