Task1 地理数据分析常用工具

智慧海洋竞赛实践 专题一

Task1 地理数据分析常用工具

shapely和geopandas是做地理空间数据的分析

folium和kepler.gl是进行地理数据可视化的工具

geohash是将经纬度坐标进行数据编码的方式

1. 常用工具

1.1 Shapely库

库简介

Shapely是python中开源的空间几何对象库，支持point(点)、curve(线)和surface(面)等基本几何对象类型以及相关空间操作。另外，几何对象类型的特征分别有interior(内部)、boundary(边界)和exterior(外部)。

1.1.1 空间数据模型

基本几何对象类型	对应的方法(类)
point类型	Point
curve类型	LineString和LinearRing
surface类型	Polygon
point集合	MultiPoint
curves集合	MultiLineString
surface集合	MultiPolygon

1.1.2 几何对象的一些功能特性

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)

1.1.3 Python操作

1. 导入包

from shapely import geometry as geo
from shapely import wkt 
from shapely import ops
import numpy as np 

2. Point

class Point(coordinates)

点的构建

# point有三种赋值方式，具体如下
# 给x,y,z的数值
point_1 = geo.Point(0.5,0.5)
point_2 = geo.Point((0,0))  
point_3 = geo.Point(point_1) #把point_1的坐标点给point_3

# 获取点的坐标
print(list(point_3.coords))
print(point_3.x)
print(point_3.y)
# 功能特性：几何对象可以和numpy.array互相转换
print(np.array(point))

RUN:

[(0.5, 0.5)]
0.5
0.5
[0.5 0.5]

点的可视化

geo.GeometryCollection([point_1,point_2])

3. LineStrings

class LineString(coordinates)

线的构建

LineStrings构造函数传入参数是2个或多个点元组

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)))#该方法即可求线线距离也可以求线点距离
# 最长距离0
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)))

RUN:

两个几何对象之间的距离: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)]

center = line.centroid #几何中心 结果是一个点
geo.GeometryCollection([line,center])

bbox = line.envelope #envelope可以求几何对象的最小外接矩形
geo.GeometryCollection([line,bbox]) 

rect = line.minimum_rotated_rectangle #最小旋转外接矩形
geo.GeometryCollection([line,rect])

DouglasPucker算法的实现

（算法具体内容可参考：https://blog.csdn.net/wangqinghao/article/details/8207070）

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算法

4. LinearRings

class LinearRing(coordinates)

LineStrings构造函数传入参数是2个或多个点元组

ring = geo.polygon.LinearRing([(0, 0), (1, 1), (1, 0)])
print(ring.length)#相比于刚才的LineString的代码示例，其长度现在是3.41，是因为其序列是闭合的
print(ring.area)

RUN:

3.414213562373095
0.0

和3LineString相比，这是一个闭合的环，但他的面积仍然是0，线的面积当然为0.

5. Polygon

class Polygon(shell[, holes=None])
Polygon接受两个位置参数，第一个位置参数是和LinearRing一样，是一个有序的point元组。第二个位置参数是可选的序列，其用来指定内部的边界

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))  #外围坐标点

6. 集合

MultiPoint、MultiLineString、MultiPolygon分别指的是多个点、多个linestring和多个polygon形成的集合

7. 几何对象关系

关系判断	具体含义
object.contains(other)	如果object的外部没有其他点，或者至少有一个点在该object的内部，则返回True a.contains(b)与 b.within(a)的表达是等价的
object.crosses(other)	如果一个object与另一个object是内部相交的关系而不是包含的关系，则返回True
object.disjoint(other)	如果该对象与另一个对象的内部和边界都不相交则返回True
object.intersects(other)	如果该几何对象与另一个几何对象只要相交则返回True
object.convex_hull	返回包含对象中所有点的最小凸多边形（凸包）

coords = [(0, 0), (1, 1)]
print(geo.LineString(coords).contains(geo.Point(0.5, 0.5)))#线与点的关系
print(geo.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)])

print( geo.LineString(coords).crosses(geo.LineString([(0, 1), (1, 0)])))
print(geo.Point(0, 0).disjoint(geo.Point(1, 1)))
print( geo.LineString(coords).intersects(geo.LineString([(0, 1), (1, 0)])))

# 在下图中即为在给定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])

# 返回对象与对象之间的交集
polygon1 = Polygon( [(0, 0), (0, 2), (2, 2), (2, 0), (0, 0)])
hull1.intersection(polygon1)
#返回对象与对象之间的并集
hull1.union(polygon1)
#返回对象与对象之间的补集
hull1.difference(polygon1) 

8. 与numpy和python数组之间的关系

point、LineRing和LineString提供numpy数组接口，可以进行转换numpy数组

from shapely.geometry import asPoint,asLineString,asMultiPoint,asPolygon
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格式

另外还有一些非常有用但是不属于某个类方法的函数，如有需要可以在官网查阅

• ops.nearest_points 求最近点
• ops.split 分割线
• ops.substring 求子串
• affinity.rotate 旋转几何体
• affinity.scale 缩放几何体
• affinity.translate 平移几何体

1.2 geopandas库

库简介

GeoPandas提供了地理空间数据的高级接口，它让使用python处理地理空间数据变得更容易。GeoPandas扩展了pandas使用的数据类型，允许对几何类型进行空间操作。几何运算由shapely执行。Geopandas进一步依赖fiona进行文件访问，依赖matplotlib进行绘图。

geopandas和pandas一样，一共有两种数据类型：

• GeoSeries
• GeoDataFrame

它们继承了pandas数据结构的大部分方法，这两个数据结构可以当做地理空间数据的存储器，shapefile文件的pandas呈现。

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()

#根据每一个polygon的pop_est不同，便可以用python绘制图表显示不同国家的人数
fig, ax = plt.subplots(figsize=(9,6),dpi = 100)
world.plot('pop_est',ax = ax,legend = True)
plt.show()

相关学习资料：

https://www.cnblogs.com/giserliu/p/4988615.html

集锦：
https://mp.weixin.qq.com/mp/appmsgalbum?action=getalbum&__biz=MzA3ODYwNDkzOQ==&scene=1&album_id=1342860388945444864&count=3#wechat_redirect

1.3 Folium库

folium可以满足我们平时常用的热力图、填充地图、路径图、散点标记等高频可视化场景.

folium也可以通过flask让地图和我们的数据在网页上显示，极其便利。

# 没有安装folium的要先安装
import folium
import os
#首先，创建一张指定中心坐标的地图，这里将其中心坐标设置为北京。zoom_start表示初始地图的缩放尺寸，数值越大放大程度越大
m=folium.Map(location=[39.9,116.4],zoom_start=10)
m

1.3.1 相关参数

参数	输入	用途
location	tuple/list	用于控制初始地图中心点的坐标，格式为【纬度，经度】/(纬度，经度)，默认为None
width	int/str	- int型，传入的是地图宽度的像素值 - str型，传入的是地图宽度的百分比，形式为‘xx%’,默认100%
height	/	控制地图的高度，格式同width
tiles	str	用于控制绘图调用的地图样式，默认为‘OpenStreetMap’
max_zoom	int	控制地图可以放大程度的上限，默认为18
attr	str	当在tiles中使用自选URL内的osm时使用，用于给自选osm命名
control_scale	bool	控制是否在地图上添加比例尺，默认为False
no_touch	bool	控制地图是否禁止接受来自设备的触控时间，默认为False

1.3.2 绘制热力图

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()

m=folium.Map([48,5],tiles='stamentoner',zoom_start=6)
HeatMap(data).add_to(m)
m

folium的其他使用可以参考知乎的这篇文章，较为全面。
https://www.zhihu.com/question/33783546

1.4 Kepler.gl库

Kepler.gl与folium类似，也是是一个图形化的数据可视化工具，基于Uber的大数据可视化开源项目deck.gl创建的demo app。目前支持3种数据格式：CSV、JSON、GeoJSON。

基础教程：https://sspai.com/post/55655

1.4.1 获取数据

#获取文件夹中的数据
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'] = pd.to_numeric(tmp_df['direction'])
    return tmp_df
# 读入数据
df=get_data(r'hy_round1_train_20200102','train')

1.4.2 修改经纬度&时间

# 平面坐标转经纬度，供初赛数据使用
# 选择标准为NAD83 / California zone 6 (ftUS) (EPSG:2230)
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]
    return time_str_reformat

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'))

1.4.3 可视化展示

学习教程：https://www.yht7.com/news/127888

1.5 GeoHash库

在对于经纬度进行数据分析和特征提取时常用到的是GeoHash编码，该编码方式可以将地理经纬度坐标编码为由字母和数字所构成的短字符串，它具有如下特性：

1. 层级空间数据结构，将地理位置用矩形网格划分，同一网格内地理编码相同
2. 只要编码长度足够长，可以表示任意精度的地理位置坐标
3. 编码前缀匹配的越长，地理位置越邻近。

参考文档：https://blog.csdn.net/zhufenghao/article/details/85568340

1.5.1 Python操作

def geohash_encode(latitude, longitude, precision=12):
    """
    Encode a position given in float arguments latitude, longitude to
    a geohash which will have the character count precision.
    """
    lat_interval, lon_interval = (-90.0, 90.0), (-180.0, 180.0)
    base32 = '0123456789bcdefghjkmnpqrstuvwxyz'
    geohash = []
    bits = [16, 8, 4, 2, 1]
    bit = 0
    ch = 0
    even = True
    while len(geohash) < precision:
        if even:
            mid = (lon_interval[0] + lon_interval[1]) / 2
            if longitude > mid:
                ch |= bits[bit]
                lon_interval = (mid, lon_interval[1])
            else:
                lon_interval = (lon_interval[0], mid)
        else:
            mid = (lat_interval[0] + lat_interval[1]) / 2
            if latitude > mid:
                ch |= bits[bit]
                lat_interval = (mid, lat_interval[1])
            else:
                lat_interval = (lat_interval[0], mid)
        even = not even
        if bit < 4:
            bit += 1
        else:
            geohash += base32[ch]
            bit = 0
            ch = 0
    return ''.join(geohash)

#调用Geohash函数
DF[DF['ID']==1].apply(lambda x: geohash_encode(x['lat'], x['lon'], 7), axis=1)

1.5.2 注意事项

GeoHash的主要价值在于将二维的经纬度坐标信息编码到了一维的字符串中，在做地理位置索引时只需要匹配字

符串即可，便于缓存、信息压缩。在使用大数据工具（例如Spark）进行数据挖掘聚类时，GeoHash显得更加便捷

和高效。

但是使用GeoHash还有一些注意事项：

1. 由于GeoHash使用Z形曲线来顺序填充空间的，而Z形曲线在拐角处会有突变，这表现在有些相邻的网格的编码前缀比其他网格相差较多，因此利用前缀匹配可以找到一部分邻近的区域，但同时也会漏掉一些。
2. 一个网格内部所有点会共用一个GeoHash值，在网格的边缘点会匹配到可能较远但是GeoHash值相同的点，而本来距离较近的点却没有匹配到。这种问题可以这样解决：适当增加GeoHash编码长度，并使用周围的8个近邻编码来参与，因为往往只使用一个GeoHash编码可能会有严重风险！

2. 作业

2.1 背景介绍

实验数据共11000条(其中7000条训练数据、2000条testA、2000条testB)渔船轨迹北斗数据；

比赛数据在本次组队学习中只用到了hy_round1_testA_20200102与hy_round1_train_20200102文件。

其中DF.csv和df_gpd_change.pkl 分别是Task1中所需要的数据。DF.csv是将轨迹数据进行异常处理之后的数据，而df_gpd_change.pkl是将异常处理之后的数据进行douglas-peucker算法进行压缩之后的数据。

2.1 基础作业

1.尝试去使用kepler.gl可视化来分析不同类型船舶AIS数据的分布情况，并为接下来的特征工程的提取建立基础

 """8G内存的电脑请慎重执行下面代码"""
#创建一个KeplerGl对象
map1 = KeplerGl(height=800)
#激活KeplerGl对象到jupyter的窗口中
map1

map1.add_data(df,name='layer1')

在基本设置美化一下

2.2 进阶作业

尝试在原有剔除异常点的数据（DF）中保留douglas-peucker算法所识别的关键点的数据，删除douglas-peucker未保存的数据，并尝试对这些坐标点进行geohash编码。
感兴趣的小伙伴一起来做一下~

3. 参考内容

https://tianchi.aliyun.com/forum/postDetail?spm=5176.12586969.1002.3.163c24d1HiGiFo&postId=110644