ofo数据获取坐标转换生成shp数据

随着共享单车的不断增加以及其重要性，很多人开始通过共享单车大数据对城市生活进行研究与分析。前段时间做ofo数据分析时一直被数据所困扰，通过学习与借鉴其他学者的研究，利用python对ofo单车数据进行爬取。

相对于其他爬取程序，该程序有以下特点：

本程序利用python2.7进行编写，考虑到arcpy模块使用版本问题；
由于爬取的是高德地图上面的ofo定位数据，所以又将火星坐标系下的定位数据转换为WGS84坐标系下的定位数据；
同时将转换后的数据导出为shp点数据，方便操作与研究。

首先呢，该ofo定位数据爬取程序是通过登录ofo网页爬取数据的，登录地址：https://common.ofo.so/newdist/?Journey。进行爬取前必须获取认证信息token，获取方式如下图所示：

在下面的程序代码中用的是我的token，为了长久使用，大家可以自己登录获取token。

同时该程序是通过定义矩形区域进行爬取的，所以要事先查询要爬取区域的左上角经纬度与右下角经纬度，调用start函数时填写该数据。

详细代码如下：

#!coding=utf-8
from __future__ import division
import requests
import datetime
import threading
import json
import os
import pandas as pd
import numpy as np
import time
import sqlite3
import math
import arcpy
from requests.packages.urllib3.exceptions import InsecureRequestWarning
from requests_toolbelt.multipart.encoder import MultipartEncoder
from concurrent.futures import ThreadPoolExecutor
requests.packages.urllib3.disable_warnings(InsecureRequestWarning)class Crawler:def __init__(self):self.start_time = datetime.datetime.now()self.db_name = "file:database?mode=memory&cache=shared"self.save_path = "./data/" + datetime.datetime.now().strftime("%Y%m%d")self.file_name = datetime.datetime.now().strftime("%Y%m%d_%H%M") + "_ofo"self.lock = threading.Lock()self.total = 0self.done = 0self.bikes_count = 0self.x_pi = 3.14159265358979324 * 3000.0 / 180.0self.pi = 3.1415926535897932384626  # πself.a = 6378245.0  # 长半轴self.ee = 0.00669342162296594323  # 偏心率平方self.message = os.path.isdir(self.save_path)def get_nearby_bikes(self, args):try:url = "https://san.ofo.so/ofo/Api/nearbyofoCar"headers = {'charset': "utf-8",'Accept': '*/*','Accept-Encoding': 'gzip, deflate','Accept-Language': 'zh-CN','Content-Length': '516','Content-Type': 'multipart/form-data; boundary=----ofo-boundary-MC4wOTk1ODUy','Host': 'san.ofo.so','Origin': 'https://common.ofo.so','Referer': 'https://common.ofo.so/newdist/','User-Agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/51.0.2704.79 Safari/537.36 Edge/14.14393'}self.request(headers, args, url)except Exception as ex:print(ex)def request(self, headers, args, url):multipart_encoder = MultipartEncoder(fields={"token": "e923c290-0d27-11e7-b9d2-b5e857d3318f","source": "0","source-version": "9999",# "lat": "36.103235","lat": str(args[0]),# "lng":"103.709681""lng": str(args[1])# file为路径},boundary='----ofo-boundary-MC4wOTk1ODUy')response = requests.request("POST", url, headers=headers,timeout=30,verify=False,data=multipart_encoder)with self.lock:with self.connect_db() as c:try:decoded = json.loads(response.text)['values']['info']['cars']self.done += 1for x in decoded:self.bikes_count += 1at = self.gcj02_to_wgs84_lat(x['lng'], x['lat'])#将火星坐标转换为WGS84--纬度on = self.gcj02_to_wgs84_lng(x['lng'], x['lat'])#将火星坐标转换为WGS84--经度c.execute("INSERT OR IGNORE INTO ofo VALUES (%d,'%s',%f,%f)" % (int(time.time()) * 1000, x['carno'], at, on))timespent = datetime.datetime.now() - self.start_timepercent = self.done / self.totaltotal = timespent.total_seconds() / percentprint("位置 %s, 未去重单车数量 %s, 进度 %0.2f%%, 速度 %0.2f个/分钟, 总时间 %s, 剩余时间 %s" % (args, self.bikes_count, percent * 100, self.done / timespent.total_seconds() * 60, total,total - timespent.total_seconds()))except Exception as ex:print(ex)def connect_db(self):return sqlite3.connect(self.db_name)def generate_create_table_sql(self, brand):return '''CREATE TABLE {0}("bikeId" VARCHAR(12),lat DOUBLE,lng DOUBLE,CONSTRAINT "{0}_bikeId_lat_lon_pk"PRIMARY KEY (bikeId, lat, lng));'''.format(brand)#创建shp点数据def CreateFeaturclass(self, savepath, featurename, spatial):if arcpy.Exists(savepath + '\\' + featurename + '.shp') == False:arcpy.CreateFeatureclass_management(savepath, featurename, 'POINT', '', '', '', spatial)else:pass#添加字段def AddField(self, savepath, featurename):arcpy.AddField_management(savepath + '\\' + featurename + '.shp', 'bikeid', 'TEXT')arcpy.AddField_management(savepath + '\\' + featurename + '.shp', 'lon', 'TEXT')arcpy.AddField_management(savepath + '\\' + featurename + '.shp', 'lat', 'TEXT')#遍历点并添加字段值def InsertRow(self, savepath, featurename, data):Insercur = arcpy.InsertCursor(savepath + '\\' + featurename + '.shp')for value in range(1,len(data)+1):point = arcpy.Point()newrow = Insercur.newRow()point.X = float(data.head(value)['lng'][value-1])point.Y = float(data.head(value)['lat'][value-1])newrow.setValue('Id', value)newrow.setValue('bikeid',data.head(value)['bikeId'][value-1])newrow.setValue('lon',data.head(value)['lng'][value-1])newrow.setValue('lat',data.head(value)['lat'][value-1])pointGeo = arcpy.PointGeometry(point)newrow.shape = pointGeoInsercur.insertRow(newrow)def gcj02_to_wgs84_lng(self, lng1, lat1):if self.out_of_china(lng1, lat1):  # 判断是否在国内return lng1, lat1dlng = self.transformlng(lng1 - 105.0, lat1 - 35.0)radlat = lat1 / 180.0 * self.pimagic = math.sin(radlat)magic = 1 - self.ee * magic * magicsqrtmagic = math.sqrt(magic)dlng = (dlng * 180.0) / (self.a / sqrtmagic * math.cos(radlat) * self.pi)mglng = lng1 + dlngreturn lng1 * 2 - mglngdef gcj02_to_wgs84_lat(self, lng1, lat1):if self.out_of_china(lng1, lat1):  # 判断是否在国内return lng1, lat1dlat = self.transformlat(lng1 - 105.0, lat1 - 35.0)radlat = lat1 / 180.0 * self.pimagic = math.sin(radlat)magic = 1 - self.ee * magic * magicsqrtmagic = math.sqrt(magic)dlat = (dlat * 180.0) / ((self.a * (1 - self.ee)) / (magic * sqrtmagic) * self.pi)mglat = lat1 + dlatreturn lat1 * 2 - mglatdef transformlat(self, lng1, lat1):ret = -100.0 + 2.0 * lng1 + 3.0 * lat1 + 0.2 * lat1 * lat1 + \0.1 * lng1 * lat1 + 0.2 * math.sqrt(math.fabs(lng1))ret += (20.0 * math.sin(6.0 * lng1 * self.pi) + 20.0 *math.sin(2.0 * lng1 * self.pi)) * 2.0 / 3.0ret += (20.0 * math.sin(lat1 * self.pi) + 40.0 *math.sin(lat1 / 3.0 * self.pi)) * 2.0 / 3.0ret += (160.0 * math.sin(lat1 / 12.0 * self.pi) + 320 *math.sin(lat1 * self.pi / 30.0)) * 2.0 / 3.0return retdef transformlng(self, lng1, lat1):ret = 300.0 + lng1 + 2.0 * lat1 + 0.1 * lng1 * lng1 + \0.1 * lng1 * lat1 + 0.1 * math.sqrt(math.fabs(lng1))ret += (20.0 * math.sin(6.0 * lng1 * self.pi) + 20.0 *math.sin(2.0 * lng1 * self.pi)) * 2.0 / 3.0ret += (20.0 * math.sin(lng1 * self.pi) + 40.0 *math.sin(lng1 / 3.0 * self.pi)) * 2.0 / 3.0ret += (150.0 * math.sin(lng1 / 12.0 * self.pi) + 300.0 *math.sin(lng1 / 30.0 * self.pi)) * 2.0 / 3.0return retdef out_of_china(self, lng1, lat1):return not (lng1 > 73.66 and lng1 < 135.05 and lat1 > 3.86 and lat1 < 53.55)def group_data(self):print("正在导出数据")conn = self.connect_db()self.export_to_shp(conn, "ofo")def export_to_shp(self, conn, brand):spRef = arcpy.SpatialReference(4326)df = pd.read_sql_query("SELECT * FROM %s" % brand, conn, parse_dates=True)self.CreateFeaturclass(self.save_path, self.file_name, spRef)self.AddField(self.save_path, self.file_name)self.InsertRow(self.save_path, self.file_name, df)print(brand)print ("去重后数量")print (len(df))def start(self, top_lng, top_lat, bottom_lng, bottom_lat):while True:self.__init__()if self.message == False:os.makedirs(self.save_path)#创建路径try:with self.connect_db() as c:c.execute(self.generate_create_table_sql('ofo'))except Exception as ex:print(ex)passexecutor = ThreadPoolExecutor(max_workers=100)print("Start")self.total = 0offset = 0.002lat_range = np.arange(float(top_lat), float(bottom_lat), -offset)for lat in lat_range:lng_range = np.arange(float(top_lng), float(bottom_lng), offset)for lon in lng_range:self.total += 1executor.submit(self.get_nearby_bikes, (lat, lon))executor.shutdown()self.group_data()#是否继续运行always_run = Falseif not always_run:breakwaittime = 1print("等待%s分钟后继续运行" % waittime)time.sleep(waittime * 60)if __name__ == '__main__':c = Crawler()c.start(103.686592, 36.114191, 103.741781, 36.091515)#爬取范围（左上角，右下角经纬度）print("完成")

最后将导出为shp数据在ArcGIS中展示如下图所示：

![将导出的数据在ArcGIS中展示](https://img-blog.csdn.net/2018082118581859?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzI0NjU1NzAx/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70) ![属性表](https://img-blog.csdn.net/20180821190051209?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzI0NjU1NzAx/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70)

参考文献
https://github.com/SilverBooker/ofoSpider

本文仅供参考学习，有不到之处，望大家谅解

ofo数据获取坐标转换生成shp数据相关推荐

【GlobalMapper精品教程】053：打开dbf文件并生成有坐标系的shp数据
本文讲解在globalmapper汇总打开dbf文件并生成有坐标系的shp数据. 文章目录一.dbf文件解读二.打开dbf文件二.另存为shp文件一.dbf文件解读我们可以通过Excel或F ...
arcgis分隔图层重复出文件_已知坐标点txt文件在ArcGIS中转换成shp数据的两种方法...
在平时工作中,经常会遇到只有txt坐标的勘测定界图等数据,通过以下操作步骤可将txt数据转换成shp数据. 方法一 1.打开txt数据,如下图所示,该数据为面数据,坐标系为西安80,投影方式是3度分带 ...
arcmap shp导出cad无反应_如何使用ArcMap将Excel数据转换为shp数据
1. 概述对ArcMap而言,除了shapefile等数据源,还可以将包含地理位置的表格数据以 XY 坐标的形式添加到地图中,比如Excel格式的数据,如果包含有坐标数据在里面,就可以通过导入XY数 ...
OSM OpenStreetMap 获取城市路网数据及转为ESRI shp数据的方法
目录一.OSM简介二.获取OSM数据方法方法1:官网下载方法2:通过ArcGIS Editor for OpenStreetMap插件下载方法3:通过QGIS下载三.获取OSM城市路网数据 ...
C#不用ArcEngine，生成Shp文件(五)---------读取.shx文件和生成.shx文件
这一篇来写一下.shx文件的读取跟生成.测试数据下载地址为:http://download.csdn.net/detail/gis0911178/9650967 在第一篇时候有介绍. 索引文件(.sh ...
FME将用SHP数据对栅格影像数据进行裁剪
提示:文章写完后,目录可以自动生成,如何生成可参考右边的帮助文档文章目录前言一.需求二.实操 1.数据准备 2.FME的流程设置 3.裁剪的结果总结前言利用FME对数据进行操作与处理的过 ...
数据处理（1）：十万加SHP数据重载导入标准字段图层
使用工具:Arcgis 时常我们会从客户.网络得到数据,但是数据里面字段名于公司规定得字段名有出入,此时就需要进行数据调整,首先我们打开软件"ArcCatalog" 在任意文件夹下 ...
JavaWeb和WebGIS学习笔记（三）——GeoServer 发布shp数据地图
JavaWeb和WebGIS学习笔记(三)--GeoServer 发布shp数据地图系列链接: Java web与web gis学习笔记(一)--Tomcat环境搭建 Java web与web gi ...
MATLAB导入任意省市地图（Shp数据）并从EXCEL中导入经纬度数据在图上手动加点
目录导入任意各省市Shp数据从EXCEL导入经纬度数据并在图上手动加点导入任意各省市Shp数据 1.首先下载任意省市Shp数据,详细步骤请看: 转载:三步教你免费下载省,市,区县行政区Shp数据 ...

ofo数据获取坐标转换生成shp数据

ofo数据获取坐标转换生成shp数据相关推荐

最新文章

热门文章