作為一名程式設計師兼職業餘空間資料分析師，空間資料的處理一直是一抹揮之不去的烏雲。雖然GIS軟體視覺化的操作介面、包羅永珍的功能，已經能解決日常工作中幾乎所有問題，但對於身為程式設計師的我來說，一切不能用鍵盤上26個字母解決的問題，都是不科學的~所以這些年我一直致力于思考一個問題，如何把Arcgis的功能，用程式碼的方式實現。

我瞭解到，Geopandas是一個極好的工具，今天我們就以一個案例來管窺一二。

 資料準備&任務明確

我有一張上海的街鎮地圖，一份上海的房源資料和一份上海餐廳資料。

 

今天，我們需要製作每一個房源的1KM緩衝區，計算緩衝區內有多少家餐廳

，並且計算上海每個街鎮內有多少家餐廳，然後繪製在地圖上展示。

這是一個理論上簡單但操作繁瑣的過程，有些同學可能已經知道GIS的操作：

匯入地圖資料——匯入兩份EXCEL——並關聯經緯度——獲得點資料——並儲存點資料圖層——然後對房源點資料進行緩衝區製作——儲存緩衝區圖層——分別使用空間關聯把餐廳資料關聯到上海街鎮和房源緩衝區圖層中——並分別儲存圖層……

好了好了，我已經說不下去了，我們來看看程式碼怎麼解決吧~

 

 匯入資料 

先匯入各種包

import pandas as

 pd
import geopandas as gpd
from shapely import geos
from shapely.geometry import Point
import fiona
import matplotlib.pyplot as plt
from fiona.crs import from_epsg,from_string

我們來觀察下python中用到的地理處理工具包的名字：Geopandas. 顧名思義，這個包可以讓我們像使用pandas一樣處理地理資料，大家可以想象一下這過程，一定如絲般順滑！

Geopandas其實是各種地理資料分析包的集大成者，包括shapely,Fiona等,當然還有資料分析相關的包numpy,pandas，所以，在這些提及的包中的功能，都可以混合使用，切換不留痕跡。這裡就一併匯入了。

 

匯入地圖shp

shanghai_map = gpd.GeoDataFrame.from_file('./上海街鎮/上海街鎮.shp') #讀取shapfile資料為geodataframe格式

Geopandas提供了一種資料格式叫GeoDataFrame，用直白的話概括就是DataFrame加了一列資料，表達地理資訊。匯入功能Geopandas底層呼叫的是Fiona包，所以，一些基本引數和可以匯入的資料格式，可以參考Fiona的說明文件。

 

檢視匯入後的資料

和DataFrame如出一轍，但多了一列Geometry來存放地理資訊，繪圖看一下：

嗯！是長這樣！

 

接下來，匯入房源和餐廳資料

 *這兩份資料是CSV格式，匯入成DataFrame，我們發現數據中含有經緯度欄位，我們可以根據這兩個欄位，也把資料轉換成GeoDataFrame格式

通過經緯度轉換點資料：

def point_to_geo(df,lon,lat):
    df['geometry'] = gpd.GeoSeries(list(zip(df[lon],df[lat]))).apply(Point) #識別經緯度，轉換點資料
    df = gpd.GeoDataFrame(df) #轉換Geodataframe格式
    df.crs = {'init':'epsg:4326'} #定義座標系WGS84
    del df[lon]
    del df[lat]
    return df


house_data = point_to_geo(house_data,'lon_WGS','lat_WGS') #轉換Geodataframe格式
restaurant_data = point_to_geo(restaurant_data,'lon_WGS','lat_WGS') #轉換Geodataframe格式

 

*官方文件請看這裡：

繪圖，看看成功與否——以上海街鎮為底圖，兩份點資料疊加在底圖上

base = shanghai_map.plot(color='lightyellow',edgecolor='black',figsize=(15, 15)) #畫底圖
restaurant_data.plot(ax=base,marker='o', color='green', markersize=5) #在底圖上疊加餐廳點資料
house_data.plot(ax=base,marker='o', color='red', markersize=5) #在底圖上疊加房源點資料
plt.gca().xaxis.set_major_locator(plt.NullLocator())#去掉x軸刻度
plt.gca().yaxis.set_major_locator(plt.NullLocator())#去掉y軸刻度
plt.savefig('./map.png',dpi=400) #儲存圖片

*這裡設定下圖形的基本引數，畫圖功能是基於matplotlib，所以它的一些寫法和功能是通用的。

 

 顯示OK，那麼就進入核心環節~

 分析過程 

在製作緩衝區和空間關聯之前，先需要對圖形進行投影變換，這裡仍然寫一個函式：

def wgs84_to_CGCS2000(df,code):
    result = df.to_crs(from_epsg(code))
    return result

shanghai_map_pcs = wgs84_to_CGCS2000(shanghai_map,4549)
house_data_pcs = wgs84_to_CGCS2000(house_data,4549)
restaurant_data_pcs = wgs84_to_CGCS2000(restaurant_data,4549)

 *只要知道投影座標系的ESPG程式碼，就可以任意轉換，我們使用GSC2000座標系。

 

接著，構建緩衝區圖層，使用buffer()方法，引數為緩衝區的半徑。

house_data_buffer = house_data_pcs.buffer(1000) #建立一公里緩衝區

base = shanghai_map_pcs.plot(color='lightyellow',edgecolor='black',figsize=(15, 15))
house_data_buffer.plot(ax=base,color='gray', markersize=5,alpha=0.5)
plt.gca().xaxis.set_major_locator(plt.NullLocator()) #去掉x軸刻度
plt.gca().yaxis.set_major_locator(plt.NullLocator()) #去年y軸刻度
plt.savefig('./map2.png',dpi=400) #儲存圖片

構建完成，同樣輸出成地圖檢視效果：

增加了透明度，我們可以看到市中心的房源密度較高。生成的資料為GeoSeries，資料為緩衝區的地理資訊。這裡需要使用任意方法把地理資訊和其他欄位匹配在一起。

buffer_temp = house_data_pcs[['name','geometry']]
buffer_temp['geometry'] = house_data_buffer
house_data_buffer = buffer_temp

 

下面，我們使用空間關聯，連線餐廳點資料，並獲取餐廳的個數。

spacial_join_restaurant = gpd.sjoin(house_data_buffer,restaurant_data_pcs,how='left',op='contains') #空間連線
spacial_join_restaurant = spacial_join_restaurant.groupby(['name']).count()['title'].to_frame().reset_index() #聚合計算個數
spacial_join_restaurant.columns = ['name','restaurant_count'] #更改列名，方便操作
buffer_result = pd.merge(house_data_pcs,spacial_join_restaurant,left_on='name',right_on='name',how='left') #欄位匹配

 *方法是sjoin()，功能和GIS相同，分相交，包含和被包含三種。我們使用包含進行連線，並使用pandas的groupby()方法分組計數。

這樣我們就完成了空間關聯和計算，我們可以使用地圖顯示，並新增圖例：

base = shanghai_result.plot(column='restaurant_count', cmap='Oranges',scheme = 'fisher_jenks'
                             ,legend=True,edgecolor='black',figsize=(15, 15)) #按個數多少疊加底色
plt.gca().xaxis.set_major_locator(plt.NullLocator()) #去掉x軸刻度
plt.gca().yaxis.set_major_locator(plt.NullLocator()) #去年y軸刻度
plt.savefig('./map3.png',dpi=400) #儲存圖片

 

這樣就可以繪製出以街鎮為單位的點密度圖：

最後，我們需要把資料匯出儲存。我們可以儲存為shapfile格式，這裡推薦大家直接儲存為csv格式，方便分享。也可以根據需要，轉成任意格式。

shanghai_result = wgs84_to_CGCS2000(shanghai_result,4326) #轉換成地理座標系
shanghai_result.columns = ['town','region','geometry','restaurant_count'] #更改列名為英文，因為資料儲存對中文支援不佳
buffer_result = wgs84_to_CGCS2000(buffer_result,4326) #轉換成地理座標系

shanghai_result.to_csv('./result/shanghai_result.csv') #儲存成csv
shanghai_result.to_file('./result/shanghai_result.shp') #儲存成shapfile
buffer_result.to_csv('./result/buffer_result.csv') #儲存成csv
buffer_result.to_file('./result/buffer_result.shp') #儲存成shapfile

 

 結果檢查 

我們把儲存的資料匯入Arcgis，可以正常使用~太棒了！

 

是不是很實用呢？所以，對於空間資料分析中一些重複操作，工程化的作業，建議大家可以考慮用程式碼實現，高效，便捷。

作者簡介：本文作者資料俠倪家禹，城市資料團特約撰稿人，資料分析師（Python）微專業學員，喜歡用資料探勘生活中的小祕密。對資料研究有著敏銳的洞察力，善於把複雜的問題簡單化，簡單的問題流程化。希望大家通過資料感受生活的魅力。

宣告：本文為作者投稿，版權歸其個人所有。

---

 熱 文 推 薦  

☞ 「傻瓜」才能寫出好程式碼！

☞ 漫畫 | 從搬家到容器技術 Docker 應用場景解析

☞ Hacker News 12 月招聘趨勢：React 已連續霸榜 19 個月

☞ 從傾家蕩產到身價百億，這個85後只用了8年

☞ 難逃寒冬裁員的“大追殺”，30 歲女碼農該何去何從？

☞ OpenStack 2018 年終盤點

☞ 拼多多黃崢給陸奇“兼職”，欲挖掘這類AI人才

☞ 老程式設計師肺腑忠告：千萬別一輩子靠技術生存！

print_r('點個好看吧！');
var_dump('點個好看吧！');
NSLog(@"點個好看吧！");
System.out.println("點個好看吧！");
console.log("點個好看吧！");
print("點個好看吧！");
printf("點個好看吧！\n");
cout << "點個好看吧！" << endl;
Console.WriteLine("點個好看吧！");
fmt.Println("點個好看吧！");
Response.Write("點個好看吧！");
alert("點個好看吧！")
echo "點個好看吧！"

喜歡就點選“好看”吧！