利用python繪製分析路易斯安那州巴吞魯日市的人口密度格局

有我之境發表於2022-01-09
Python

前言

資料來源於王法輝教授的GIS和數量方法，以後有空，我會利用python來實現裡面的案例，這裡向王法輝教授致敬。

繪製普查人口密度格局

使用屬性查詢提取區邊界

import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import arcpy
from arcpy import env
plt.style.use('ggplot')#使用ggplot樣式
%matplotlib inline#輸出線上圖片

plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.sans-serif'] = ['SimHei']# 替換sans-serif字型為黑體
plt.rcParams['axes.unicode_minus'] = False   # 解決座標軸負數的負號顯示問題

regions = gpd.GeoDataFrame.from_file('../Census.gdb',layer='County')
regions

資料預覽

BRTrt = regions[regions.NAMELSAD10=='East Baton Rouge Parish']

投影

BRTrt = BRTrt.to_crs('EPSG:26915')
BRTrt.crs

投影資訊

BRTrt.to_file('BRTrt.shp')

裁剪資料

Tract = gpd.GeoDataFrame.from_file('../Census.gdb',layer='Tract')
Tract = Tract.to_crs('EPSG:26915')

TractUtm = gpd.GeoDataFrame.from_file('TractUtm.shp')
BRTrtUtm = gpd.GeoDataFrame.from_file('BRTrt.shp')

# Set workspace
env.workspace = r"MyProject"

# Set local variables
in_features = "TractUtm.shp"
clip_features = "BRTrt.shp"
out_feature_class = "BRTrtUtm.shp"
xy_tolerance = ""

# Execute Clip
arcpy.Clip_analysis(in_features, clip_features, out_feature_class, xy_tolerance)

裁剪

計算面積和人口密度

BRTrtUtm = gpd.GeoDataFrame.from_file('BRTrtUtm.shp')
BRTrtUtm['area'] = BRTrtUtm.area/1000000

## 計算人口密度
BRTrtUtm['PopuDen'] = BRTrtUtm['DP0010001']/BRTrtUtm['area']
BRTrtUtm.to_file('BRTrtUtm.shp')

描述統計

BRTrtUtm['PopuDen'].describe()

人口密度描述統計

人口密度圖

fig = plt.figure(figsize=(12,12)) #設定畫布大小
ax = plt.gca()
ax.set_title("巴吞魯日市2010年人口密度模式",fontsize=24,loc='center')
BRTrtUtm.plot(ax=ax,column='PopuDen',linewidth=0.5,cmap='Reds'
              ,edgecolor='k',legend=True,)
# plt.savefig('巴吞魯日市2010年人口密度模式.jpg',dpi=300)

plt.show()

巴吞魯日市2010年人口密度模式


分析同心環區的人口密度格式

生成同心環

## 兩種方法生成多重緩衝區的閾值
dis = list(np.arange(2000,26001,2000))
dis
dis = list(range(2000,26001,2000))
dis

多重緩衝區數值

## 真的特別神奇distances只有這樣寫列表才可以執行

 
# Set local variables
inFeatures = "BRCenter"
outFeatureClass = "rings.shp"
distances = [2000, 4000, 6000, 8000, 10000, 
             12000, 14000, 16000, 18000, 
             20000, 22000, 24000, 26000]
bufferUnit = "meters"

# Execute MultipleRingBuffer
arcpy.MultipleRingBuffer_analysis(inFeatures, outFeatureClass, distances, bufferUnit, "", "ALL")

相交

try:
    # Set the workspace (to avoid having to type in the full path to the data every time)
    arcpy.env.workspace = "MyProject"    
    
    # Process: Find all stream crossings (points)
    inFeatures = ["rings", "BRTrtUtm"]
    intersectOutput = "TrtRings.shp"   
    arcpy.Intersect_analysis(inFeatures, intersectOutput,)
 
except Exception as err:
    print(err.args[0])

TrtRings = gpd.GeoDataFrame.from_file('TrtRings.shp')
TrtRings['area'] = TrtRings.area/1000000

TrtRings['EstPopu'] = TrtRings['PopuDen'] * TrtRings['POLY_AREA']

融合

arcpy.env.workspace = "C:/data/Portland.gdb/Taxlots"

# Set local variables
inFeatures = "TrtRings"
outFeatureClass = "DissRings.shp"
dissolveFields = ["distance"]
statistics_fields = [["POLY_AREA","SUM"], ["PopuDen","SUM"]]

# Execute Dissolve using LANDUSE and TAXCODE as Dissolve Fields


arcpy.Dissolve_management(inFeatures, outFeatureClass, dissolveFields, statistics_fields,)

DissRings = gpd.GeoDataFrame.from_file('DissRings.shp')
DissRings

圈層

DissRings['PopuDen'] = DissRings['SUM_PopuDe'] / DissRings['SUM_POLY_A']
DissRings.set_index('distance',inplace=True)

DissRings['PopuDen'].plot(kind='bar',x='distance',
                          xlabel='',figsize=(8,6))
plt.savefig('同心環人口密度圖.jpg',dpi=300)
plt.show()

同心環人口密度圖

要素轉點

# Set environment settings
env.workspace = "BR.gdb"

#  Set local variables
inFeatures = "BRBlkUtm"
outFeatureClass = "BRBlkPt.shp"

# Use FeatureToPoint function to find a point inside each park
arcpy.FeatureToPoint_management(inFeatures, outFeatureClass, "INSIDE")

標識

env.workspace = "MyProject"

# Set local parameters
inFeatures = "BRBlkPt"
idFeatures = "DissRings"
outFeatures = "BRBlkPt_Identity.shp"

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

資料篩選

BRBlkPt_Identity = gpd.GeoDataFrame.from_file('BRBlkPt_Identity.shp')
BRBlkPt_Identity.shape

資料大小

BRBlkPt_Identity.tail()

需要剔除的資料

## 選取資料
BRBlkPt_Identity = BRBlkPt_Identity[~(BRBlkPt_Identity['distance']==0.0)]

資料分組

rigs_data = pd.DataFrame(BRBlkPt_Identity.groupby(by=['distance'])['POP10'].sum(),columns=['POP10'])
rigs_data.reset_index(inplace=True)
rigs_data

圈層人口資料

資料連線

EstPopu = BRBlkPt_Identity[['distance','SUM_POLY_A','SUM_PopuDe']]
PopuDen = pd.merge(rigs_data,EstPopu,how='inner',left_on='distance',right_on='distance')

## 刪除重複值，按理來說，應該沒有重複值了，可以試試外連線
PopuDen.drop_duplicates(inplace = True)

分析和比較環形區人口和密度估值

PopuDen.set_index('distance',inplace=True)
PopuDen['EstPopu'] = PopuDen['SUM_PopuDe'] / PopuDen['SUM_POLY_A']
PopuDen['PopuDen1'] = PopuDen['POP10'] / PopuDen['SUM_POLY_A']

PopuDen['EstPopu'].plot(figsize=(10,6),marker='o',xlabel='距離(米)',ylabel='密度(人/平方千米)')
PopuDen['PopuDen1'].plot(marker='s',xlabel='距離(米)',ylabel='密度(人/平方千米)')

plt.legend(['基於街道','基於普查區'])
plt.savefig('基於普查區和街區資料的人口密度模式對比.jpg',dpi=300)
plt.show()

基於普查區和街區資料的人口密度模式對比

總結

2022年的第一次寫筆記，寫的不是很好，而且發現許多問題，比如就是geopandas裡面的area和arcpy裡面的area不一樣，可能是演算法不一樣，面積要使用投影座標系，我相信這個應該沒有人不知道了吧，要對ArcGIS Pro裡面的arcpy大讚。最近感謝也比較多，比如疫情，已經有點常態化，很影響我們的生活了。心懷感恩，希望我們都有美好的未來。春燕歸，巢於林木。接下來一段時間，我要忙我的畢業論文，可能會比較忙，需要資料的可以聯絡我。

相關文章