利用Python統計各鄉鎮各地類面積

利用Python統計各鄉鎮各地類面積

導入相關模塊

import arcpy
from arcpy import env
import arcpy.sa as sa
import numpy as np 
import pandas as pd 
from dbfread import DBF
import geopandas as gpd
import matplotlib.pyplot as plt

env.workspace = r"F:\ArcGIS\ArcGIS文件\EX02數據準備"

建立文件地理數據庫

# Import system modules
import os
import sys


# Set local variables
out_folder_path = "F:\ArcGIS\ArcGIS文件\EX02數據準備"
out_name = "myfgdb.gdb"
# Execute CreateFileGDB
arcpy.CreateFileGDB_management(out_folder_path, out_name)
print(out_name+'建立成功!')

批量投影

咱們進行面積計算都是須要投影座標系的，這裏的地理座標統一成投影座標CGCS2000_3_Degree_GK_Zone_35。

# Local variables:
input_features = ["面狀要素_土地利用數據_Clip.shp", "面狀要素_鄉鎮行政邊界.shp", "點狀要素_鄉鎮駐地.shp"]
out_workspace = "F:\\ArcGIS\\ArcGIS文件\\EX02數據準備\\myfgdb.gdb"


# Process: 批量投影
arcpy.BatchProject_management(input_features, out_workspace, "PROJCS['CGCS2000_3_Degree_GK_Zone_35',GEOGCS['GCS_China_Geodetic_Coordinate_System_2000',DATUM['D_China_2000',SPHEROID['CGCS2000',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Gauss_Kruger'],PARAMETER['False_Easting',35500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',105.0],PARAMETER['Scale_Factor',1.0],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]", "", "")
## 我也不知道爲何這個投影名稱這麼長，應該用CGCS2000_3_Degree_GK_Zone_35就好了

空間鏈接

## 因爲個人鎮界沒有屬性信息，因此進行空間鏈接
target_features = r"F:/ArcGIS/ArcGIS文件/EX02數據準備/myfgdb.gdb/面狀要素_鄉鎮行政邊界_shp"
join_features = r"F:/ArcGIS/ArcGIS文件/EX02數據準備/myfgdb.gdb/點狀要素_鄉鎮駐地_shp"
out_feature_class = r"F:/ArcGIS/ArcGIS文件/EX02數據準備/myfgdb.gdb/鄉鎮邊界_空間鏈接_shp"

arcpy.SpatialJoin_analysis(target_features, join_features, out_feature_class)

標識

# IdentityWells.py
# Description: Simple example showing use of Identity tool
 
# Import system modules
import arcpy
from arcpy import env

# Set the workspace
env.workspace = "F:/ArcGIS/ArcGIS文件/EX02數據準備/myfgdb.gdb"

# Set local parameters
inFeatures = "面狀要素_土地利用數據_Clip_shp"
idFeatures = "鄉鎮邊界_空間鏈接_shp"
outFeatures = "鄉鎮_土地_面積"

# Process: Use the Identity function
arcpy.Identity_analysis (inFeatures, idFeatures, outFeatures)

數據透視

data = gpd.read_file(env.workspace,layer='鄉鎮_土地_面積')

data.head()

## 數據清洗
df = data[['鄉鎮名','GUIZHOU2_1','Shape_Area']]
df.head()

df.shape

## 重命名列
df = df.rename(columns={'鄉鎮名':'鄉鎮','GUIZHOU2_1':'地類','Shape_Area':'面積'})

df.head()

## 數據透視
pivot_df = pd.pivot_table(df,index=['鄉鎮','地類'],values='面積',aggfunc=np.sum)

#顯示全部列
pd.set_option('display.max_columns', None)

#顯示全部行
pd.set_option('display.max_rows', None)

#設置value的顯示長度爲100，默認爲50
pd.set_option('max_colwidth',100)

pivot_df

## 計算百分比
pivot_df['佔比'] = round(pivot_df.面積 / pivot_df.groupby(level=0).面積.transform(sum) * 100,2).astype(str) + '%'

pivot_df

保存數據

df.to_excel('各村土地利用分佈表.xlsx')

pivot_df.to_excel('各村土地利用數據透視表.xlsx')

總結

先是利用arcpy進行空間分析，而後利用pandas、geopandas等進行數據透視，注意計算面積必定用的是投影座標系，由於我數據放在gdb裏面，要是shp文件的話，還有進行計算幾何。