一、簡介
shp格式的檔案是地理資訊領域最常見的檔案格式之一，很好的結合了向量資料與對應的標量資料，而在Python中我們可以使用pyshp來完成建立shp檔案的過程，本文將從如何從高德地圖獲取向量資訊開始，最終構造出相應的shp檔案，並利用R中的leaflet進行視覺化；
二、資料獲取及清洗
2.1 資料獲取
首先我們需要從高德地圖獲取所關注物件的向量資訊，這裡點資料我們選擇重慶軌道交通站點，線我們選擇重慶軌道交通線路，面我們選擇重慶市三峽博物館，考慮到只是簡單演示小規模採集資料，因此選擇selenium作為資料爬取的工具，首先我們需要操縱模擬瀏覽器開啟高德地圖查詢內容的頁面（即query帶有關鍵詞），這樣做的目的是讓我們的瀏覽器載入所需介面對應的cookies，方便之後直接進行向量資訊的採集，如下面這個頁面：
執行下面的程式碼啟動瀏覽器，接著觀察會不會出現滑塊驗證碼（筆者親測有很大的概率觸發驗證碼）：
from selenium import webdriver from tqdm import tqdm import time option = webdriver.ChromeOptions() '''這個實驗引數用於避免被高德檢測到為driver驅動的瀏覽器''' option.add_experimental_option('excludeSwitches', ['enable-automation']) browser = webdriver.Chrome(options=option) '''訪問指定網址拿到cookies''' browser.get('https://www.amap.com/search?query=軌道交通3號線&city=500000&geoobj=106.477496%7C29.407019%7C106.642291%7C29.665101&zoom=12')
這時若出現下列驗證碼則手動接觸即可（考慮到爬蟲並不是本文重點因此沒有花費時間編寫模擬滑動滑塊的程式碼）：
在滑塊解除後，我們就可以批量獲取軌道線路向量資訊，程式碼如下，注意每輪執行間隔調久一些防止被ban：
'''這個字典存放所有原始的json資料''' rawSHP = {} crtLines = ['軌道交通1號線','軌道交通2號線','軌道交通3號線','軌道交通4號線','軌道交通5號線','軌道交通6號線','軌道交通10號線', '軌道交通3號線北延伸段(空港線)','軌道交通6號線支線','軌道交通環線'] for line in tqdm(crtLines): browser.get(f'https://www.amap.com/service/poiInfo?query_type=TQUERY&pagesize=20&pagenum=1&qii=true&cluster_state=5&need_utd=true&utd_sceneid=1000&div=PC1000&addr_poi_merge=true&is_classify=true&zoom=12&city=500000&geoobj=106.477496%7C29.394307%7C106.642291%7C29.677779&keywords={line}') '''這裡從網頁內容標籤中抽取json部分內容''' rawSHP[line] = eval(browser.find_elements_by_xpath("//pre")[0].text) time.sleep(8)
這樣我們就得到對應重慶軌道交通線路和站點的原始json資料，接下來類似上面的做法，獲取中國三峽博物館向量資訊：
browser.get('https://www.amap.com/service/poiInfo?query_type=TQUERY&pagesize=20&pagenum=1&qii=true&cluster_state=5&need_utd=true&utd_sceneid=1000&div=PC1000&addr_poi_merge=true&is_classify=true&zoom=12&city=500000&geoobj=106.477496%7C29.394307%7C106.642291%7C29.677779&keywords=中國三峽博物館') '''這裡從網頁內容標籤中抽取json部分內容''' museumSX = eval(browser.find_elements_by_xpath("//pre")[0].text)
經過上面的步驟我們就得到了所需內容的原始格式，接下來進行清洗；
2.2 資料清洗
首先提取點資料，rawSHP為字典，鍵為線路名稱，值為所對應包含的全部內容，我們需要的經緯度資訊就包含在其中，以環線為例：
按照上圖箭頭所指的路徑便可找到對應的站點名稱name和經緯度xy_coords，而對於線資料，如下圖：
同樣可以找到對應每個折點的經度xs與緯度ys，對於面資料，在museumSX變數下data->poi_list->domain_list中name屬性為'aoi'的元素中可以找到其對應的面向量資訊：
獲悉所需資料的位置之後，接下來我們在寫入shp檔案的過程中同時完成清洗過程，在此之間首先需要介紹pyshp中寫出shp檔案相關的用法；
三、寫出shp檔案
3.1 用pyshp寫出shp檔案
pyshp是以純Python程式碼的方式對ESRI shapefiles檔案進行讀寫、編輯等操作的模組，以用法方便快捷功能高效強大著稱，寫出時使用到其中的Writer類，其主要有三個引數：
target：檔案最終存出的具體位置及檔名稱
shapeType：int型，用於決定檔案型別，型別與數字對應關係如下：
autoBalance：int型，建議傳入1，即定義的屬性有秩序的自動跟隨定義的要素之後，避免出現錯亂；
而pyshp中的Writer物件有如下常用方法：
field：用於建立跟隨向量要素的屬性表欄位，其name引數用於定義欄位名；fieldType引數用於控制資料型別，'C'代表字串，‘N’代表數值型，‘F’代表浮點型，‘L’代表bool型，‘D’代表日期；引數size為字元型，用於控制資料長度，最大限制為‘2046’
point：傳入點的經度與緯度
line：傳入單條或多條線每個折點的經緯度
poly：傳入面中對應每個邊界點的經緯度
除了上述三種最基本的，還有很多傳入其他格式向量資訊的方法，本文未使用到不再贅述；
record：傳入屬性表對應欄位的值
close：在最後存出檔案時呼叫
因為我們爬取的資料來自高德地圖，因此如果有轉換座標系的需求，可以使用下列程式碼完成百度座標、火星座標系、wgs84之間的互轉：
import math x_pi = 3.14159265358979324 * 3000.0 / 180.0 pi = 3.1415926535897932384626 # π a = 6378245.0 # 長半軸 ee = 0.00669342162296594323 # 偏心率平方 class LngLatTransfer(): def __init__(self): pass def gcj02_to_bd09(self, lng, lat): """ 火星座標系(GCJ-02)轉百度座標系(BD-09) 谷歌、高德——>百度 :param lng:火星座標經度 :param lat:火星座標緯度 :return: """ z = math.sqrt(lng * lng + lat * lat) + 0.00002 * math.sin(lat * x_pi) theta = math.atan2(lat, lng) + 0.000003 * math.cos(lng * x_pi) bd_lng = z * math.cos(theta) + 0.0065 bd_lat = z * math.sin(theta) + 0.006 return [bd_lng, bd_lat] def bd09_to_gcj02(self, bd_lon, bd_lat): """ 百度座標系(BD-09)轉火星座標系(GCJ-02) 百度——>谷歌、高德 :param bd_lat:百度座標緯度 :param bd_lon:百度座標經度 :return:轉換後的座標列表形式 """ x = bd_lon - 0.0065 y = bd_lat - 0.006 z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * x_pi) theta = math.atan2(y, x) - 0.000003 * math.cos(x * x_pi) gg_lng = z * math.cos(theta) gg_lat = z * math.sin(theta) return [gg_lng, gg_lat] def wgs84_to_gcj02(self, lng, lat): """ WGS84轉GCJ02(火星座標系) :param lng:WGS84座標系的經度 :param lat:WGS84座標系的緯度 :return: """ if self.out_of_china(lng, lat): # 判斷是否在國內 return [lng, lat] dlat = self._transformlat(lng - 105.0, lat - 35.0) dlng = self._transformlng(lng - 105.0, lat - 35.0) radlat = lat / 180.0 * pi magic = math.sin(radlat) magic = 1 - ee * magic * magic sqrtmagic = math.sqrt(magic) dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi) dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi) mglat = lat + dlat mglng = lng + dlng return [mglng, mglat] def gcj02_to_wgs84(self, lng, lat): """ GCJ02(火星座標系)轉GPS84 :param lng:火星座標系的經度 :param lat:火星座標系緯度 :return: """ if self.out_of_china(lng, lat): return [lng, lat] dlat = self._transformlat(lng - 105.0, lat - 35.0) dlng = self._transformlng(lng - 105.0, lat - 35.0) radlat = lat / 180.0 * pi magic = math.sin(radlat) magic = 1 - ee * magic * magic sqrtmagic = math.sqrt(magic) dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi) dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi) mglat = lat + dlat mglng = lng + dlng return [lng * 2 - mglng, lat * 2 - mglat] def bd09_to_wgs84(self, bd_lon, bd_lat): lon, lat = self.bd09_to_gcj02(bd_lon, bd_lat) return self.gcj02_to_wgs84(lon, lat) def wgs84_to_bd09(self, lon, lat): lon, lat = self.wgs84_to_gcj02(lon, lat) return self.gcj02_to_bd09(lon, lat) def _transformlat(self, lng, lat): ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + \ 0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng)) ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 * math.sin(2.0 * lng * pi)) * 2.0 / 3.0 ret += (20.0 * math.sin(lat * pi) + 40.0 * math.sin(lat / 3.0 * pi)) * 2.0 / 3.0 ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 * math.sin(lat * pi / 30.0)) * 2.0 / 3.0 return ret def _transformlng(self, lng, lat): ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \ 0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng)) ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 * math.sin(2.0 * lng * pi)) * 2.0 / 3.0 ret += (20.0 * math.sin(lng * pi) + 40.0 * math.sin(lng / 3.0 * pi)) * 2.0 / 3.0 ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 * math.sin(lng / 30.0 * pi)) * 2.0 / 3.0 return ret def out_of_china(self, lng, lat): """ 判斷是否在國內，不在國內不做偏移 :param lng: :param lat: :return: """ return not (lng > 73.66 and lng < 135.05 and lat > 3.86 and lat < 53.55)
3.2 寫出shp檔案
點檔案：
思路是初始化Writer物件之後，利用迴圈從rawSHP字典中抽取所有的站點名稱、經緯度以及對應線路，因此屬性表中建立欄位name用於儲存站點名稱，route欄位用於存放線路資訊，具體程式碼如下（注意匯入名需為shapefile，即pyshp）：
'''shp檔案寫出部分'''
import shapefile
'''抽取經緯度-站點名稱-線路名稱三元組''' all_points = [] for line in tqdm(rawSHP.keys()): for line_ in rawSHP[line]['data']['busline_list']: for station in line_['stations']: '''這裡完成火星座標系向WGS84的轉換''' all_points.append([transfer.gcj02_to_wgs84(lng=float(station['xy_coords'].split(';')[0]), lat=float(station['xy_coords'].split(';')[1])), station['name'], line]) '''去重''' all_points_ = [] for item in all_points: if item not in all_points_: all_points_.append(item) '''初始化Writer物件''' w_point = shapefile.Writer(r'C:\Users\hp\Desktop\shp寫出\重慶軌道交通站點向量資料', autoBalance=shapefile.POINT) '''建立站點名稱欄位''' w_point.field('name','C') '''建立線路名稱欄位''' w_point.field('route','C') for item in all_points_: '''寫入點要素''' w_point.point(item[0][0],item[0][1]) '''追加屬性值''' w_point.record(name=item[1],route=item[2]) '''關閉存出檔案''' w_point.close()
輸出目錄中也包含了我們所需的檔案：
在arcgis中檢視：
成功~
接下來是線檔案：
'''shp檔案寫出部分''' import shapefile w_line = shapefile.Writer(r'C:\Users\hp\Desktop\shp寫出\重慶軌道交通線路向量資料', autoBalance=1) w_line.field('name','C') for key in rawSHP.keys(): lines = [] for idx in range(len(rawSHP[key]['data']['busline_list'])): lines.append([]) for lng,lat in zip(rawSHP[key]['data']['busline_list'][idx]['xs'].split(','),rawSHP[key]['data']['busline_list'][idx]['ys'].split(',')): lines[-1].append(transfer.gcj02_to_wgs84(lng=float(lng),lat=float(lat))) for l in lines: '''每段線路要素單獨寫出''' w_line.line([l]) '''追加對應的線路名稱''' w_line.record(name=key) w_line.close()
在arcgis中檢視線檔案：
面檔案
rawPoly = museumSX['data']['poi_list'][0]['domain_list'][14]['value'].split('_') rawPoly = [transfer.gcj02_to_wgs84(float(rawPoly[i].split(',')[0]),float(rawPoly[i].split(',')[1])) for i in range(len(rawPoly))] w_polygon = shapefile.Writer(r'C:\Users\hp\Desktop\shp寫出\三峽博物館面向量資料', autoBalance=shapefile.POLYGON) w_polygon.field('name','C') w_polygon.poly([rawPoly]) w_polygon.record(name='中國三峽博物館') w_polygon.close()
在arcgis中檢視：
可以與高德網頁上的形狀對比，非常吻合，至此，我們就完成了shp檔案的生成，下面我們簡單的在R中用leaflet進行視覺化，這裡選用Carto的底圖（WGS84座標系），對應的R程式碼如下：
rm(list=ls()) library(rgdal) library(leaflet) library(ggplot2) library(readxl) setwd('C:\\Users\\hp\\Desktop\\shp寫出') crt <- readOGR('重慶軌道交通線路向量資料.shp') crt_station <- readOGR('重慶軌道交通站點向量資料.shp') museum <- readOGR('三峽博物館面向量資料.shp') #用迴圈的方式疊加線 m <- leaflet() %>% addTiles(urlTemplate = '//{s}.basemaps.cartocdn.com/light_nolabels/{z}/{x}/{y}{r}.png') for(i in 1:length(crt@lines)){ m <- m %>% addPolylines(fortify(crt@lines[[i]])$long,lat=fortify(crt@lines[[i]])$lat,fillColor = sample(colors(),1),color = sample(colors(),1),weight=2) } #疊加點 m <- m %>% addCircleMarkers(lng=data.frame(crt_station@coords)$coords.x1,lat=data.frame(crt_station@coords)$coords.x2, radius=1, weight = 1, opacity = 1, color = 'red', fillOpacity = 1, label=crt_station@data ) #疊加面 m <- m %>% addPolygons(lng=fortify(museum)$long, lat=fortify(museum)$lat) m
放大後可以看到位於中山四路附近的三峽博物館，跟高德地圖上的對比一下，還是我們的底圖比較素雅~：
以上就是本文的全部內容，如有疏漏之處望指出。