PySpark空间计算学习总结—（6）PySpark实战例一处理7500w人口大数据

一、前沿说明

需求：目前手上有 7500w 条全国人口数据（点数据，有四个字段，分别是时间、经度、维度、标记），csv 格式。现在要找到哪些数据在广东省内。

分析：可以将人口数据与广东省行政区(这里用的shp文件)进行求交集，即点数据和面数据求交集，则可以得出在广东省境内的人口数据。
方法：这里用 GDAL 中的 ogr 求交集方法。首先读取广东省行政区划数据（广东省行政区划数据为 shp 格式，是一个面数据），构造成一个可以用来进行求交集并集等计算的几何结构（这里用 wkt，因为 python 无法序列化 GDAL 的对象为一个列表）；然后读取人口数据，筛选出不符合要求的数据，将点数据转为一个 rdd；最后将点数据和面数据都转为空间元素之后，再进行求交集运算。
注意：这里麻烦的点有两个：一个是读取 shp 文件，将其构造为 wkt 结构，另一个点是将点数据和面数据转为空间元素。

二、具体实现说明

目前只是在较少数据中实现了（人口数据只用了 500 条）
1、首先要下载安装所需要的依赖

import ogr,os
import shapely.geometry as sg
import findspark
findspark.init()
import pyspark

2、代码实现
（1）读取 shp 文件，并将其构造为一个用来计算的几何结构

# 读取shp文件——广东省区域
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open("./data/gd/Guangdong.shp",0)
layer = dataSource.GetLayerByIndex(0)
# 构造一个用来计算的几何结构，注意这里用wkt，是因为Python无法序列化GDAL的对象为一个列表
featArr = []
layer.ResetReading()
for f in layer:
    featArr.append((f.GetField("NL_NAME_2"),
                    f.geometry().ExportToIsoWkt()))

（2）过滤点数据，将点数据转化为一个 rdd

# 过滤方法——点
def isPoint(line):
    pntline = line.split("\t")
    try:
        wkt = "POINT({0} {1})".format(float(pntline[1]),float(pntline[2]))
        geom = ogr.CreateGeometryFromWkt(wkt)
        return geom.IsValid()
    except:
        return False
rdd = sc.textFile("./data/20190510.csv").filter(lambda line:isPoint(line))

（3）将点数据和面数据都转为空间元素之后，再进行求交集运算。

# 快速查询，将查询出来的数据插入到 resultsAll 里面
def myMap(pnt,featArr):
    resultsAll = []
    for feat in featArr:
        key = feat[0]
        geo = ogr.CreateGeometryFromWkt(feat[1])
        pntline = pnt.split("\t")
        wkt = "POINT({0} {1})".format(float(pntline[1]),
                                      float(pntline[2]))
        pntGeom = ogr.CreateGeometryFromWkt(wkt)
        if geo.Contains(pntGeom):
            resultsAll.append(geo.Intersects(pntGeom))
    return resultsAll
# 计算
res = rdd.filter(lambda line : myMap(line,featArr)).collect()

三、在 PyCharm 里面运行

注意：这里没有把文件地址写死在程序里面，通过 PyCharm 里面传入参数，并且把结果保存到本地文件中了。
首先在运行代码里面在 Edit-Configurations 里面设置传入参数：

则完整的代码修改为：

import ogr,os,sys
import shapely.geometry as sg
import findspark
findspark.init()
from pyspark import SparkConf, SparkContext

# 将计算结果保存到本地
if __name__ == '__main__':

    if len(sys.argv) != 4:
        print("Usage: wordcount <input> <input> <output>",file=sys.stderr)
        sys.exit(-1)

    conf = SparkConf()
    sc = SparkContext(conf=conf)

    # 001—读取shp文件——广东省区域
    driver = ogr.GetDriverByName('ESRI Shapefile')
    dataSource = driver.Open(sys.argv[1], 0)
    layer = dataSource.GetLayerByIndex(0)

    # 构造一个用来计算的几何结构，注意这里用wkt，是因为Python无法序列化GDAL的对象为一个列表
    featArr = []
    layer.ResetReading()
    for f in layer:
        featArr.append((f.GetField("NL_NAME_2"),f.geometry().ExportToIsoWkt()))

    # 过滤方法———点
    def isPoint(line):
        pntline = line.split("\t")
        try:
            wkt = "POINT({0} {1})".format(float(pntline[1]), float(pntline[2]))
            geom = ogr.CreateGeometryFromWkt(wkt)
            return geom.IsValid()
        except:
            return False

    rdd = sc.textFile(sys.argv[2]).filter(lambda line: isPoint(line))

    #  快速查询，将查询出来的数据插入到 resultsAll 里面
    def myMap(pnt, featArr):
        resultsAll = []
        for feat in featArr:
            key = feat[0]
            geo = ogr.CreateGeometryFromWkt(feat[1])
            pntline = pnt.split("\t")
            wkt = "POINT({0} {1})".format(float(pntline[1]),
                                          float(pntline[2]))
            pntGeom = ogr.CreateGeometryFromWkt(wkt)
            if geo.Contains(pntGeom):
                resultsAll.append(geo.Intersects(pntGeom))
        return resultsAll

	# 将数据保存到本地
    def saveFile():
        sc.textFile(sys.argv[2])\
            .filter(lambda line: isPoint(line))\
            .filter(lambda line : myMap(line,featArr))\
            .saveAsTextFile(sys.argv[3])

    saveFile()

    sc.stop()

注：上述代码可以直接放到服务器上运行：在服务器中新建一个 spark0101.py 文件，直接在服务器终端输入（路径也可以是 hdfs 的路径）：

./spark-submit --master local[2] --name spark0101.py /home/spark0101.py file:///home/data/guangdong.shp file:///home/data/20190510.csv file:///home/data/result