一、前沿说明
需求:目前手上有 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