智慧海洋建设-Task1 地理数据分析常用工具
在地理空间数据分析中，常会用到许多地理分析的工具，在本模块中主要是针对常用的shapely、geopandas、folium、kepler.gl、geohash等工具进行简单介绍。其中shapely和geopandas是做地理空间数据的分析很好的工具，而folium和kepler.gl是进行地理数据可视化的工具，geohash是将经纬度坐标进行数据编码的方式。通过了解不同的方法将有助于我们去思考如何在现有的工具下去做数据的分析和特征的提取功能
学习目标
1.了解和学习shapely和geopandas的基本功能，掌握用python中的这两个库实现几何对象之间的空间操作方法。
2.掌握folium和kepler.gl的数据可视化工具的使用。
3.学习与掌握geohash编码方法。
内容介绍
shapely
空间数据模型
几何对象的一些功能特性
Point
LineStrings
LineRings
Polygon
几何对象之间的关系
geopandas
Folium
Kepler.gl
GeoHash
注意事项
Shapely是python中开源的空间几何对象库，支持Point、Curve和Surface等基本几何对象类型以及相关空间操作。另外，几何对象类型的特征分别有interior、boundary和exterior。
空间数据模型
1.point类型对应的方法在Point类中。curve类型对应的方法在LineString和LinearRing类中。surface类型对应的方法在Polygon类中。 2.point集合对应的方法在MultiPoint类中，curves集合对应的反方在MultiLineString类中，surface集合对应的方法在MultiPolygon类中。
几何对象的一些功能特性
Point、LineString和LinearRing有一些功能非常有用。
几何对象可以和numpy.array互相转换。
可以求线的长度(length)，面的面积（area)，对象之间的距离(distance),最小最大距离(hausdorff_distance),对象的bounds数组(minx, miny, maxx, maxy)
可以求几何对象之间的关系：相交(intersect)，包含(contain)，求相交区域(intersection)等。
可以对几何对象求几何中心(centroid),缓冲区(buffer),最小旋转外接矩形(minimum_rotated_rectangle)等。
可以求线的插值点(interpolate),可以求点投影到线的距离(project),可以求几何对象之间对应的最近点(nearestPoint)
可以对几何对象进行旋转(rotate)和缩放(scale)
!conda install shapely
^C
from shapely import geometry as geo
from shapely import wkt
from shapely import ops
import numpy as np
# point有三种赋值方式，具体如下
point = geo.Point(0.5,0.5)
point_2 = geo.Point((0,0))
point_3 = geo.Point(point)
# 其坐标可以通过coords或x，y，z得到
print(list(point_3.coords))
print(point_3.x)
print(point_3.y)
#批量进行可视化
geo.GeometryCollection([point,point_2])
print(np.array(point))#可以和np.array进行互相转换
[(0.5, 0.5)]
0.5
0.5
[0.5 0.5]
#代码示例
arr=np.array([(0,0), (1,1), (1,0)])
line = geo.LineString(arr) #等同于 line = geo.LineString([(0,0), (1,1), (1,0)])
print ('两个几何对象之间的距离:'+str(geo.Point(2,2).distance(line)))#该方法即可求线线距离也可以求线点距离
print ('两个几何对象之间的hausdorff_distance距离:'+str(geo.Point(2,2).hausdorff_distance(line)))#该方法求得是点与线的最长距离
print('该几何对象的面积:'+str(line.area))
print('该几何对象的坐标范围:'+str(line.bounds))
print('该几何对象的长度:'+str(line.length))
print('该几何对象的几何类型:'+str(line.geom_type))
print('该几何对象的坐标系:'+str(list(line.coords)))
center = line.centroid #几何中心
geo.GeometryCollection([line,center])
两个几何对象之间的距离:1.4142135623730951
两个几何对象之间的hausdorff_distance距离:2.8284271247461903
该几何对象的面积:0.0
该几何对象的坐标范围:(0.0, 0.0, 1.0, 1.0)
该几何对象的长度:2.414213562373095
该几何对象的几何类型:LineString
该几何对象的坐标系:[(0.0, 0.0), (1.0, 1.0), (1.0, 0.0)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-kOeIIKcZ-1618329331199)(output_4_1.svg)]
bbox = line.envelope #envelope可以求几何对象的最小外接矩形
geo.GeometryCollection([line,bbox])
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-MoXzMd43-1618329331209)(output_5_0.svg)]
rect = line.minimum_rotated_rectangle #最小旋转外接矩形
geo.GeometryCollection([line,rect])
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-FjZSapXQ-1618329331211)(output_6_0.svg)]
pt_half = line.interpolate(0.5,normalized=True) #插值
geo.GeometryCollection([line,pt_half])
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-HUi35tkR-1618329331214)(output_7_0.svg)]
pt_half = line.interpolate(0.5,normalized=True) #插值
geo.GeometryCollection([line,pt_half])
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-86X8HLkX-1618329331217)(output_8_0.svg)]
line1 = geo.LineString([(0,0),(1,-0.2),(2,0.3),(3,-0.5),(5,0.2),(7,0)])
line1_simplify = line1.simplify(0.4, preserve_topology=False) #Douglas-Pucker算法
print(line1)
print(line1_simplify)
line1_simplify
LINESTRING (0 0, 1 -0.2, 2 0.3, 3 -0.5, 5 0.2, 7 0)
LINESTRING (0 0, 2 0.3, 3 -0.5, 5 0.2, 7 0)
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-pvpeMn3l-1618329331218)(output_9_1.svg)]
buffer_with_circle = line1.buffer(0.2) #端点按照半圆扩展
geo.GeometryCollection([line1,buffer_with_circle])
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-3fLDpdwT-1618329331219)(output_10_0.svg)]
LinearRings
class LinearRing(coordinates) LineStrings构造函数传入参数是2个或多个点元组
元组序列可以通过在第一个和最后一个索引中传递相同的值来显式关闭。否则，将第一个元组复制到最后一个索引，从而隐式关闭序列。 与LineString一样，元组序列中的重复点是允许的，但可能会导致性能上的损失，应该避免在序列中设置重复点。
# from shapely.geometry.polygon import LinearRing
ring = geo.polygon.LinearRing([(0, 0), (1, 1), (1, 0)])
print(ring.length)#相比于刚才的LineString的代码示例，其长度现在是3.41，是因为其序列是闭合的
print(ring.area)
geo.GeometryCollection([ring])
3.414213562373095
0.0
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-sA7w2UY0-1618329331222)(output_12_1.svg)]
from shapely.geometry.polygon import LinearRing
ring = geo.polygon.LinearRing([(0, 0), (1, 1), (1, 0)])
print(ring.length)#相比于刚才的LineString的代码示例，其长度现在是3.41，是因为其序列是闭合的
print(ring.area)
geo.GeometryCollection([ring])
from shapely.geometry import Polygon
polygon1 = Polygon([(0, 0), (1, 1), (1, 0)])
ext = [(0, 0), (0, 2), (2, 2), (2, 0), (0, 0)]
int = [(1, 0), (0.5, 0.5), (1, 1), (1.5, 0.5), (1, 0)]
polygon2 = Polygon(ext, [int])
print(polygon1.area)
print(polygon1.length)
print(polygon2.area)#其面积是ext的面积减去int的面积
print(polygon2.length)#其长度是ext的长度加上int的长度
print(np.array(polygon2.exterior)) #外围坐标点
geo.GeometryCollection([polygon2])
0.5
3.414213562373095
3.5
10.82842712474619
[[0. 0.]
[0. 2.]
[2. 2.]
[2. 0.]
[0. 0.]]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-i6Mktkpb-1618329331223)(output_14_1.svg)]
几何对象关系¶
一个几何对象特征分别有interior、boundary和exterior。下面的叙述直接用内部、边界和外部等名词概述
1.object.contains(other)
如果object的外部没有其他点，或者至少有一个点在该object的内部，则返回True
a.contains(b)与 b.within(a)的表达是等价的
from shapely.geometry import LineString
Point要有geo
from shapely.geometry import LineString
coords = [(0, 0), (1, 1)]
print(LineString(coords).contains(geo.Point(0.5, 0.5)))#线与点的关系
print(LineString(coords).contains(geo.Point(1.0, 1.0)))#因为line的边界不是属于在该对象的内部，所以返回是False
polygon1 = Polygon( [(0, 0), (0, 2), (2, 2), (2, 0), (0, 0)])
print(polygon1.contains(geo.Point(1.0, 1.0)))#面与点的关系
#同理这个contains方法也可以扩展到面与线的关系以及面与面的关系
geo.GeometryCollection([polygon1,geo.Point(1.0, 1.0)])
True
False
True
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-tmKrYdqi-1618329331225)(output_17_1.svg)]
print( LineString(coords).crosses(LineString([(0, 1), (1, 0)])))
print(geo.Point(0, 0).disjoint(geo.Point(1, 1)))
print( LineString(coords).intersects(LineString([(0, 1), (1, 0)])))
True
True
True
# 在下图中即为在给定6个point之后求其凸包，并绘制出来的凸包图形
points1 = geo.MultiPoint([(0, 0), (1, 1), (0, 2), (2, 2), (3, 1), (1, 0)])
hull1 = points1.convex_hull
geo.GeometryCollection([hull1,points1])
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-oSuiJICr-1618329331226)(output_19_0.svg)]
# object.intersection 返回对象与对象之间的交集
polygon1 = Polygon( [(0, 0), (0, 2), (2, 2), (2, 0), (0, 0)])
hull1.intersection(polygon1)
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-FDsSY8V2-1618329331231)(output_20_0.svg)]
#返回对象与对象之间的并集
hull1.union(polygon1)
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-yLdmCOEo-1618329331234)(output_21_0.svg)]
hull1.difference(polygon1) #面面补集
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Q8YrJ9tV-1618329331235)(output_22_0.svg)]
from shapely.geometry import asPoint,asLineString,asMultiPoint,asPolygon
import numpy as np
pa = asPoint(np.array([0.0, 0.0]))#将numpy数组转换成point格式
la = asLineString(np.array([[1.0, 2.0], [3.0, 4.0]]))#将numpy数组转换成LineString格式
ma = asMultiPoint(np.array([[1.1, 2.2], [3.3, 4.4], [5.5, 6.6]]))#将numpy数组转换成multipoint集合
pg = asPolygon(np.array([[1.1, 2.2], [3.3, 4.4], [5.5, 6.6]]))#将numpy数组转换成polygon
print(np.array(pa))#将Point转换成numpy格式
'''
from shapely.geometry import asPoint,asLineString,asMultiPoint,asPolygon
import numpy as np
pa = asPoint(np.array([0.0, 0.0]))#将numpy数组转换成point格式
la = asLineString(np.array([[1.0, 2.0], [3.0, 4.0]]))#将numpy数组转换成LineString格式
ma = asMultiPoint(np.array([[1.1, 2.2], [3.3, 4.4], [5.5, 6.6]]))#将numpy数组转换成multipoint集合
pg = asPolygon(np.array([[1.1, 2.2], [3.3, 4.4], [5.5, 6.6]]))#将numpy数组转换成polygon
print(np.array(pa))#将Point转换成numpy格式
'''
[0. 0.]
'\nfrom shapely.geometry import asPoint,asLineString,asMultiPoint,asPolygon\nimport numpy as np\npa = asPoint(np.array([0.0, 0.0]))#将numpy数组转换成point格式\nla = asLineString(np.array([[1.0, 2.0], [3.0, 4.0]]))#将numpy数组转换成LineString格式\nma = asMultiPoint(np.array([[1.1, 2.2], [3.3, 4.4], [5.5, 6.6]]))#将numpy数组转换成multipoint集合\npg = asPolygon(np.array([[1.1, 2.2], [3.3, 4.4], [5.5, 6.6]]))#将numpy数组转换成polygon\nprint(np.array(pa))#将Point转换成numpy格式\n'
geopandas
GeoPandas提供了地理空间数据的高级接口，它让使用python处理地理空间数据变得更容易。GeoPandas扩展了pandas。)使用的数据类型，允许对几何类型进行空间操作。几何运算由shapely执行。Geopandas进一步依赖fiona进行文件访问，依赖matplotlib进行绘图。
geopandas和pandas一样，一共有两种数据类型：
GeoSeries
GeoDataFrame 它们继承了pandas数据结构的大部分方法。这两个数据结构可以当做地理空间数据的存储器，shapefile文件的pandas呈现。
Shapefile文件用于描述几何体对象：点，折线与多边形。例如，Shapefile文件可以存储井、河流、湖泊等空间对象的几何位置。除了几何位置，shp文件也可以存储这些空间对象的属性，例如一条河流的名字，一个城市的温度等等。
例如，当安装geopandas库后，便可通过matplotlib直接画出当安装geopandas数据集中的世界地图
!pip install geopandas
^C
import pandas as pd
import geopandas
import matplotlib.pyplot as plt
world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))#read_file方法可以读取shape文件，转化为GeoSeries和GeoDataFrame数据类型。
world.plot()#将GeoDataFrame变成图形展示出来，得到世界地图
plt.show()
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-V90TJlIO-1618329331237)(output_26_0.png)]
!pip install folium
Requirement already satisfied: folium in c:\programdata\anaconda3\lib\site-packages (0.12.1)
Requirement already satisfied: jinja2>=2.9 in c:\programdata\anaconda3\lib\site-packages (from folium) (2.11.3)
Requirement already satisfied: requests in c:\programdata\anaconda3\lib\site-packages (from folium) (2.25.1)
Requirement already satisfied: branca>=0.3.0 in c:\programdata\anaconda3\lib\site-packages (from folium) (0.4.2)
Requirement already satisfied: numpy in c:\programdata\anaconda3\lib\site-packages (from folium) (1.19.2)
Requirement already satisfied: MarkupSafe>=0.23 in c:\programdata\anaconda3\lib\site-packages (from jinja2>=2.9->folium) (1.1.1)
Requirement already satisfied: chardet<5,>=3.0.2 in c:\programdata\anaconda3\lib\site-packages (from requests->folium) (4.0.0)
Requirement already satisfied: certifi>=2017.4.17 in c:\programdata\anaconda3\lib\site-packages (from requests->folium) (2020.12.5)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in c:\programdata\anaconda3\lib\site-packages (from requests->folium) (1.26.3)
Requirement already satisfied: idna<3,>=2.5 in c:\programdata\anaconda3\lib\site-packages (from requests->folium) (2.10)
import folium
import os
#首先，创建一张指定中心坐标的地图，这里将其中心坐标设置为北京。zoom_start表示初始地图的缩放尺寸，数值越大放大程度越大
m=folium.Map(location=[39.9,116.4],zoom_start=10)
m
Make this Notebook Trusted to load map: File -> Trust Notebook
import folium
import numpy as np
from folium.plugins import HeatMap
#先手动生成data数据，该数据格式由[纬度，经度，数值]构成
data=(np.random.normal(size=(100,3))*np.array([[1,1,1]])+np.array([[48,5,1]])).tolist()
# data
m=folium.Map([48,5],tiles='stamentoner',zoom_start=6)
HeatMap(data).add_to(m)
m
Make this Notebook Trusted to load map: File -> Trust Notebook
m=folium.Map([90,5],tiles='stamentoner',zoom_start=6)
HeatMap(data).add_to(m)
m
Make this Notebook Trusted to load map: File -> Trust Notebook
m=folium.Map([48,5],tiles='stamentoner',zoom_start=16)
HeatMap(data).add_to(m)
m
Make this Notebook Trusted to load map: File -> Trust Notebook
import pandas as pd
import geopandas as gpd
from pyproj import Proj
from keplergl import KeplerGl
from tqdm import tqdm
import os
import matplotlib.pyplot as plt
import shapely
import numpy as np
from datetime import datetime
import warnings
warnings.filterwarnings('ignore')
plt.rcParams['font.sans-serif'] = ['SimSun'] # 指定默认字体为新宋体。
#获取文件夹中的数据
def get_data(file_path,model):
assert model in ['train', 'test'], '{} Not Support this type of file'.format(model)
paths = os.listdir(file_path)
# print(len(paths))
tmp = []
for t in tqdm(range(len(paths))):
p = paths[t]
with open('{}/{}'.format(file_path, p), encoding='utf-8') as f:
next(f)
for line in f.readlines():
tmp.append(line.strip().split(','))
tmp_df = pd.DataFrame(tmp)
if model == 'train':
tmp_df.columns = ['ID', 'lat', 'lon', 'speed', 'direction', 'time', 'type']
else:
tmp_df['type'] = 'unknown'
tmp_df.columns = ['ID', 'lat', 'lon', 'speed', 'direction', 'time', 'type']
tmp_df['lat'] = tmp_df['lat'].astype(float)
tmp_df['lon'] = tmp_df['lon'].astype(float)
tmp_df['speed'] = tmp_df['speed'].astype(float)
tmp_df['direction'] = tmp_df['direction'].astype(int)#如果该行代码运行失败，请尝试更新pandas的版本
return tmp_df
# 平面坐标转经纬度，供初赛数据使用
# 选择标准为NAD83 / California zone 6 (ftUS) (EPSG:2230)，查询链接：https://mygeodata.cloud/cs2cs/
def transform_xy2lonlat(df):
x = df['lat'].values
y = df['lon'].values
p=Proj('+proj=lcc +lat_1=33.88333333333333 +lat_2=32.78333333333333 +lat_0=32.16666666666666 +lon_0=-116.25 +x_0=2000000.0001016 +y_0=500000.0001016001 +datum=NAD83 +units=us-ft +no_defs ')
df['lon'], df['lat'] = p(y, x, inverse=True)
return df
#修改数据的时间格式
def reformat_strtime(time_str=None, START_YEAR="2019"):
"""Reformat the strtime with the form '08 14' to 'START_YEAR-08-14' """
time_str_split = time_str.split(" ")
time_str_reformat = START_YEAR + "-" + time_str_split[0][:2] + "-" + time_str_split[0][2:4]
time_str_reformat = time_str_reformat + " " + time_str_split[1]
# time_reformat=datetime.strptime(time_str_reformat,'%Y-%m-%d %H:%M:%S')
return time_str_reformat
#计算两个点的距离
def haversine_np(lon1, lat1, lon2, lat2):
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2
c = 2 * np.arcsin(np.sqrt(a))
km = 6367 * c
return km * 1000
def compute_traj_diff_time_distance(traj=None):
"""Compute the sampling time and the coordinate distance."""
# 计算时间的差值
time_diff_array = (traj["time"].iloc[1:].reset_index(drop=True) - traj[
"time"].iloc[:-1].reset_index(drop=True)).dt.total_seconds() / 60
# 计算坐标之间的距离
dist_diff_array = haversine_np(traj["lon"].values[1:], # lon_0
traj["lat"].values[1:], # lat_0
traj["lon"].values[:-1], # lon_1
traj["lat"].values[:-1] # lat_1
)
# 填充第一个值
time_diff_array = [time_diff_array.mean()] + time_diff_array.tolist()
dist_diff_array = [dist_diff_array.mean()] + dist_diff_array.tolist()
traj.loc[list(traj.index),'time_array'] = time_diff_array
traj.loc[list(traj.index),'dist_array'] = dist_diff_array
return traj
#对轨迹进行异常点的剔除
def assign_traj_anomaly_points_nan(traj=None, speed_maximum=23,
time_interval_maximum=200,
coord_speed_maximum=700):
"""Assign the anomaly points in traj to np.nan."""
def thigma_data(data_y,n):
data_x =[i for i in range(len(data_y))]
ymean = np.mean(data_y)
ystd = np.std(data_y)
threshold1 = ymean - n * ystd
threshold2 = ymean + n * ystd
judge=[]
for data in data_y:
if (data < threshold1)|(data> threshold2):
judge.append(True)
else:
judge.append(False)
return judge
# Step 1: The speed anomaly repairing
is_speed_anomaly = (traj["speed"] > speed_maximum) | (traj["speed"] < 0)
traj["speed"][is_speed_anomaly] = np.nan
# Step 2: 根据距离和时间计算速度
is_anomaly = np.array([False] * len(traj))
traj["coord_speed"] = traj["dist_array"] / traj["time_array"]
# Condition 1: 根据3-sigma算法剔除coord speed以及较大时间间隔的点
is_anomaly_tmp = pd.Series(thigma_data(traj["time_array"],3)) | pd.Series(thigma_data(traj["coord_speed"],3))
is_anomaly = is_anomaly | is_anomaly_tmp
is_anomaly.index=traj.index
# Condition 2: 轨迹点的3-sigma异常处理
traj = traj[~is_anomaly].reset_index(drop=True)
is_anomaly = np.array([False] * len(traj))
if len(traj) != 0:
lon_std, lon_mean = traj["lon"].std(), traj["lon"].mean()
lat_std, lat_mean = traj["lat"].std(), traj["lat"].mean()
lon_low, lon_high = lon_mean - 3 * lon_std, lon_mean + 3 * lon_std
lat_low, lat_high = lat_mean - 3 * lat_std, lat_mean + 3 * lat_std
is_anomaly = is_anomaly | (traj["lon"] > lon_high) | ((traj["lon"] < lon_low))
is_anomaly = is_anomaly | (traj["lat"] > lat_high) | ((traj["lat"] < lat_low))
traj = traj[~is_anomaly].reset_index(drop=True)
return traj, [len(is_speed_anomaly) - len(traj)]
df=get_data(r'hy_round1_train_20200102','train')
100%|██████████████████████████████████████████████████████████████████████████████| 7000/7000 [01:24<00:00, 83.00it/s]
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-120-0513d9f9f524> in <module>
----> 1 df=get_data(r'hy_round1_train_20200102','train')
<ipython-input-118-7d9c7956a205> in get_data(file_path, model)
20 tmp_df['lon'] = tmp_df['lon'].astype(float)
21 tmp_df['speed'] = tmp_df['speed'].astype(float)
---> 22 tmp_df['direction'] = tmp_df['direction'].astype(int)#如果该行代码运行失败，请尝试更新pandas的版本
23 return tmp_df
24 # 平面坐标转经纬度，供初赛数据使用
C:\ProgramData\Anaconda3\lib\site-packages\pandas\core\generic.py in astype(self, dtype, copy, errors)
5875 else:
5876 # else, only a single dtype is given
-> 5877 new_data = self._mgr.astype(dtype=dtype, copy=copy, errors=errors)
5878 return self._constructor(new_data).__finalize__(self, method="astype")
5879
C:\ProgramData\Anaconda3\lib\site-packages\pandas\core\internals\managers.py in astype(self, dtype, copy, errors)
629 self, dtype, copy: bool = False, errors: str = "raise"
630 ) -> "BlockManager":
--> 631 return self.apply("astype", dtype=dtype, copy=copy, errors=errors)
632
633 def convert(
C:\ProgramData\Anaconda3\lib\site-packages\pandas\core\internals\managers.py in apply(self, f, align_keys, ignore_failures, **kwargs)
425 applied = b.apply(f, **kwargs)
426 else:
--> 427 applied = getattr(b, f)(**kwargs)
428 except (TypeError, NotImplementedError):
429 if not ignore_failures:
C:\ProgramData\Anaconda3\lib\site-packages\pandas\core\internals\blocks.py in astype(self, dtype, copy, errors)
624
625 if dtype is not None:
--> 626 dtype = pandas_dtype(dtype)
627
628 # may need to convert to categorical
C:\ProgramData\Anaconda3\lib\site-packages\pandas\core\dtypes\common.py in pandas_dtype(dtype)
1797 # raise a consistent TypeError if failed
1798 try:
-> 1799 npdtype = np.dtype(dtype)
1800 except SyntaxError as err:
1801 # np.dtype uses `eval` which can raise SyntaxError
TypeError: First element of field tuple is neither a tuple nor str
df=transform_xy2lonlat(df)
df['time']=df['time'].apply(reformat_strtime)
df['time']=df['time'].apply(lambda x: datetime.strptime(x,'%Y-%m-%d %H:%M:%S'))
#这一个cell的代码不用运行，DF.csv该数据已经放到了github上面给出的附件数据里面
#对轨迹进行异常点剔除，对nan值进行线性插值
ID_list=list(pd.DataFrame(df['ID'].value_counts()).index)
DF_NEW=[]
Anomaly_count=[]
for ID in tqdm(ID_list):
df_id=compute_traj_diff_time_distance(df[df['ID']==ID])
df_new,count=assign_traj_anomaly_points_nan(df_id)
df_new["speed"] = df_new["speed"].interpolate(method="linear", axis=0)
df_new = df_new.fillna(method="bfill")
df_new = df_new.fillna(method="ffill")
df_new["speed"] = df_new["speed"].clip(0, 23)
Anomaly_count.append(count)#统计每个id异常点的数量有多少
DF_NEW.append(df_new)
DF=pd.concat(DF_NEW)
#读取github的数据
DF=pd.read_csv('DF.csv')
line= shapely.geometry.LineString(np.array(df[df['ID']=='11'][['lon','lat']]))
ax=gpd.GeoSeries([line]).plot(color='red')
ax = gpd.GeoSeries([line]).simplify(tolerance=0.000000001).plot(color='blue',
ax=ax,
linestyle='--')
LegendElement = [plt.Line2D([], [], color='red', label='简化前'),
plt.Line2D([], [], color='blue', linestyle='--', label='简化后')]
# 将制作好的图例映射对象列表导入legend()中，并配置相关参数
ax.legend(handles = LegendElement,
loc='upper left',
fontsize=10)
# ax.set_ylim((-2.1, 1))
# ax.axis('off')
print('化简前数据长度：'+str(len(np.array(gpd.GeoSeries([line])[0]))))
print('化简后数据长度：'+str(len(np.array(gpd.GeoSeries([line]).simplify(tolerance=0.000000001)[0]))))
#定义数据简化函数。即通过shapely库将经纬度转换成LineString格式
#然后放入GeoSeries数据结构中并进行简化，最后再将所有数据放入GeoDataFrame中
def simplify_dataframe(df):
line_list=[]
for i in tqdm(dict(list(df.groupby('ID')))):
line_dict={}
lat_lon=dict(list(df.groupby('ID')))[i][['lon','lat']]
line=shapely.geometry.LineString(np.array(lat_lon))
line_dict['ID']=dict(list(df.groupby('ID')))[i].iloc[0]['ID']
line_dict['type']=dict(list(df.groupby('ID')))[i].iloc[0]['type']
line_dict['geometry']=gpd.GeoSeries([line]).simplify(tolerance=0.000000001)[0]
line_list.append(line_dict)
return gpd.GeoDataFrame(line_list)
df_gpd_change=simplify_dataframe(DF)
df_gpd_change=pd.read_pickle('df_gpd_change.pkl')
map1=KeplerGl(height=800)#zoom_start与这个height类似，表示地图的缩放程度
map1.add_data(data=df_gpd_change,name='data')
#当运行该代码后，下面会有一个kepler.gl使用说明的链接，可以根据该链接进行学习参考
map1
5791
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
