1、簡介web
shp格式的文件是地理信息領域最多見的文件格式之一，很好的結合了矢量數據與對應的標量數據，而在Python中咱們可使用pyshp來完成建立shp文件的過程，本文將從如何從高德地圖獲取矢量信息開始，最終構造出相應的shp文件，並利用R中的leaflet進行可視化；json
2、數據獲取及清洗瀏覽器
2.1 數據獲取cookie
首先咱們須要從高德地圖獲取所關注對象的矢量信息，這裏點數據咱們選擇重慶軌道交通站點，線咱們選擇重慶軌道交通線路，面咱們選擇重慶市三峽博物館，考慮到只是簡單演示小規模採集數據，所以選擇selenium做爲數據爬取的工具，首先咱們須要操縱模擬瀏覽器打開高德地圖查找內容的頁面（即query帶有關鍵詞），這樣作的目的是讓咱們的瀏覽器加載所需接口對應的cookies，方便以後直接進行矢量信息的採集，以下面這個頁面：app
運行下面的代碼啓動瀏覽器，接着觀察會不會出現滑塊驗證碼（筆者親測有很大的機率觸發驗證碼）：工具
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')
這時若出現下列驗證碼則手動接觸便可（考慮到爬蟲並非本文重點所以沒有花費時間編寫模擬滑動滑塊的代碼）：url
在滑塊解除後，咱們就能夠批量獲取軌道線路矢量信息，代碼以下，注意每輪運行間隔調久一些防止被ban：spa
'''這個字典存放全部原始的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數據，接下來相似上面的作法，獲取中國三峽博物館矢量信息：.net
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文件相關的用法；
3、寫出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
放大後能夠看到位於中山四路附近的三峽博物館，跟高德地圖上的對比一下，仍是咱們的底圖比較素雅~：
以上就是本文的所有內容，若有疏漏之處望指出。