一、簡介
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
放大后可以看到位於中山四路附近的三峽博物館，跟高德地圖上的對比一下，還是我們的底圖比較素雅~：
以上就是本文的全部內容，如有疏漏之處望指出。