基于 POI 数据的土地利用分类：以深圳市为例
0 前言
本文为【和鲸社区】Workshop | Python 气象海洋地理数据分析 Ⅲ课程第三期基于 POI 数据的土地利用分类：以深圳市为例的教案，作者为中国地质大学地理与信息工程学院教授、博导，本文转载已获得【和鲸社区】授权。
原地址如下，大家可通过【和鲸社区】ModelWhale在线运行代码，进行项目复现
https://www.heywhale.com/home/activity/detail/62de53410bbb187a397cbe76/content/3
0.1 项目背景
近年来，随着中国经济的不断发展和城市化步伐的加快，大量的土地被不断开发和征用，导致原有的土地覆盖发生了很大的变化，城市土地利用的高动态变化对决策者和管理者管理新城市社区提出了挑战。因此，如何准确识别城市土地利用覆盖，以便为土地利用规划、管理和生态修复提供科学有效的参考，有着十分重要的意义。科学规划国土空间，需要对土地利用现状和城市功能区空间分布情况进行深入分析，如何准确有效地识别城市土地利用类型也成为关键问题。
遥感作为一种可以快速获取地表覆盖信息的技术，具有快速、宏观、综合、准确、周期性、低成本等优势，已经成为城市用地信息提取与监测的重要手段 。但遥感影像中提取的特征仅仅描述了地物的自然属性，无法完全对接社会经济属性明显的城市土地利用类型。
随着互联网的快速发展，大量基于位置服务的社交媒体数据为城市土地利用分类和城市空间结构分析研究提供了丰富的数据源。社交媒体数据能够反映人类经济社会活动的内在特征，补充了高分辨率遥感影像数据无法描述的地物内部经济社会属性，有助于城市土地利用分类。
例如，兴趣点（Points-of-interest，简称POIs）数据为城市土地利用分类研究提供了大量的语义信息，包括名称、地址、功能、经纬度等，可以描述地理空间中各类商业性设施和社会服务性设施，蕴含着丰富的人文经济特征，补充了遥感影像所缺乏的语义信息，是城市空间分析中的重要基础地理数据之一。已有研究表明，POI 的分布特征可以有效地说明地块的功能。因此，研究 POI 数据支持下的城市土地利用分类，可以有效挖掘社交媒体数据中的社会经济信息，显著提高识别城市土地利用类型的精度。
0.2 任务简介
本期，我们将基于和鲸 ModelWhale 平台，手把手教大家动手学习如何利用兴趣点数据进行土地利用分类，在这里我们将向大家演示POI数据的爬取、处理、训练以及聚类可视化分析等一套完整的基本流程。我们将学习内容分为了三个notebook：
notebook1：基于百度地图开发者平台的API，利用ak秘钥爬取指定的兴趣点数据； notebook2：利用Word2Vec词向量训练模型对poi数据进行训练，进而得到poi向量； notebook3：对poi向量进行聚类以及可视化分析对其进行定性评价，并利用poi向量得到深圳市地块向量，聚类后进一步映射到地图上得到深圳市土地利用分类分布图。
1 兴趣点爬虫和实现（API）
本节使用镜像为Python 3.7 TensorFlow 2.6 PyTorch 1.8，使用的计算资源是2核8G CPU资源，Kernel类型为Python3。由于使用的镜像很基础，因此大家使用不包含torch库的镜像都可以。但爬取poi数据不涉及GPU的使用，因此大家使用CPU资源就可以了。
POI 是“ Point of Interest ”的缩写，中文可以翻译为“兴趣点”，其主要包含名称、地址、功能、经纬度等属性。在地理信息系统中，一个 POI 可以是一栋房子、一个商铺、一个邮筒、一个公交站等。POI数据分类众多，包括美食、购物、旅游景点、政府机构、交通设施等地理信息数据。利用兴趣点（POI）进行土地利用分类的前提是要先获取兴趣点数据。那么如何自动便捷地获取兴趣点数据呢？
传统的地理信息采集方法需要地图测绘人员采用精密的测绘仪器去获取一个兴趣点的经纬度，然后再标记下来，是一个非常费时费事的工作，而利用 Python 可以方便的批量调用 API 即 Application Programming Interface 应用程序接口，抓取数据、返回兴趣点的名称、经纬度坐标等数据。
本节利用百度地图 API ，实现批量抓取美食类 POI 数据。其他类型 POI 数据以此类推，以深圳市为例。那如何简单地爬取兴趣点数据呢？
1.1 任务准备
登录百度账号，在百度地图开发者平台的API控制台申请一个服务端的ak，主要用到的是Place API。
去百度地图开放平台申请密匙：http://lbsyun.baidu.com/
步骤：登录----->控制台----->注册登录成为开发者（已成为则跳过该步骤）----->创建应用------>填写应用名称----->应用类型选择浏览器端---->refer白名单输入“*”—>创建应用成功----->得到密匙（AK）
创建成功后就可以在”我的应用“中看到自己创建的ak了，把它复制下来。那么ak该如何使用呢？
1.2 爬取POI数据
1.2.1 确定url解析地址
在使用ak之前，首先我们要导入程序所需要的库，构建所要抓取的根 url ，根据百度开放平台中的服务文档的说明确定 url 的内容，可点击[百度服务平台地点检索服务介绍]查询。接下来我们要确定存放解析地址出来的坐标文件的工作目录。
[百度服务平台地点检索服务介绍](https://lbsyun.baidu.com/index.php?title=webapi/guide/webservice-placeapi)
在本例中，我们使用圆形区域进行检索，这样我们就需要先确定检索区域范围中一些圆心点的坐标。选取圆心点和半径大小时可根据获取数据的需求，覆盖到感兴趣的所有区域。
"""首先导入相关的库及一些初始的文件路径"""
from urllib.request import urlopen #调用urllib.request库里面的urlopen命令，向网页发送请求并打开网页
import re, os
import csv
import requests
import json #调用库，json是一种便于解析的格式
import time #调用库，time是时间，防止访问次数过于频繁被系统后台禁止
from urllib.parse import quote #调用urllib.parse库里面的查询命令
# 构建抓取的根URL
baseURL = 'http://api.map.baidu.com/place/v2/search?output=json&'#抓取网页的地址
ak = 'Tny0Gj8oz6LYvwguFDwoOEPbMFwrEU5v' # 上一步在百度地图中申请的ak，这里用自己申请的ak，每个ak每天都有访问限制
query = quote('美食') #爬取的POI类型是美食，quote是把中文转化成url的格式
scope = '2'#返回内容的详细程度，2是比较详尽的意思
# 确定工作目录
input_path = u'/home/mw/input/Data5799/数据集' # 爬取坐标文件所在的路径
output_path = u'/home/mw/project/' # 输出文件路径
outputFile = output_path + 'BaiduPOI_shenzhen.txt'#定义输出的POI文件，BaiduPOI_shenzhen自行命名
poiList = [] #创建一个poi列表存储爬出来的数据
# poiList.append(['name','lng','lat','type','province','city','area']) # 添加表头，包含poi点的名称、经纬度、所属类别、所属省、市、区县等信息
poiList.append('name'+','+'lng'+','+'lat'+','+'type'+','+'province'+','+'city'+','+'area') # 添加表头，包含poi点的名称、经纬度、所属类别、所属省、市、区县等信息
1.2.2 利用url链接抓取poi数据并解析
之后我们定义一个fetch小程序用于抓取 url 链接，接下来读取检索点(圆心点)的坐标文件，生成抓取列表。循环所生成的 coordinates 列表，对列表中的每一个coordinate 进行 fetch ，并存储所得数据 。
注：爬取poi时坐标默认为百度地图的坐标类型 BD09II ，否则使用其他坐标系统时(例如 WGS84 )，检索点会有偏移，我们的数据坐标系是 WGS84，因此需要先进行由 WGS84坐标系到百度坐标系的转换。调用百度地图坐标转换的 API 即可实现由非百度坐标转为百度坐标，设置转换前的坐标类型请参考百度开放平台中[坐标类型说明]。
https://lbsyun.baidu.com/index.php?title=webapi/guide/changeposition
下一步进行实际的抓取动作，根据 BaseURL ，构建抓取所用 URL，生成抓取结果并将结果赋给 response ， 同时确认此次调用返回数据的页数。
我们在编写程序之前，需要手动单次进行 POI 信息抓取，访问相应 url 后，我们发现 results 的格式是列表，因此我们在编写程序时，要先对列表里的字典进行提取，再提取字典中的信息。具体来说，对应每个坐标点，循环提取每页抓取到的 poi 信息后，将提取到的结果赋给 contents 。嵌套在上一个循环内，开始循环 contents 列表中的所需信息，对 contents 列表中的信息进行进一步提取，即 poi 的名字、经纬度和 uid。
# 定义抓取动作
# 分四步：访问-读取-解析-休眠
def fetch(url):#定义抓取的参数是url
feedback = requests.get(url)
feedback = feedback.content
data = feedback.decode("utf8") # 读取url
response = json.loads(data) # 对返回结果进行解析,将已编码的 JSON 字符串解码为 Python 对象，就是python解码json对象
time.sleep(2) # 暂停2秒，防止访问次数过于频繁被系统后台禁止return response #返回抓取结果
return response
# 官方转换函数
# 在调用百度 API 进行 POI 爬取时，由于默认输入坐标类型属于百度坐标系（BD09ll），而我们第一手数据往往并不是百度坐标。
# 所以，我们需要将WGS84坐标系转换成百度地图中使用的坐标，才可将转化后的坐标在百度地图JavaScriptAPI、静态图 API、Web 服务 API 等产品中使用。
def wgs2bd0911(wgs_x, wgs_y):
# from:1是从WGS84坐标系转换 to:5是转为bd0911
url = 'http://api.map.baidu.com/geoconv/v1/?coords={}+&from=1&to=5&output=json&ak={}'.format(str(wgs_x) + ',' + str(wgs_y), ak)
response = requests.get(url)
response = response.content
res = response.decode("utf8")
temp = json.loads(res) # 对返回结果进行解析
bd09mc_x = 0
bd09mc_y = 0
if temp['status'] == 0:
bd09mc_x = temp['result'][0]['x'] # 转换后的经度
bd09mc_y = temp['result'][0]['y'] # 转换后的纬度
return bd09mc_x, bd09mc_y
# 定义读csv文件的函数
def read_csv(filepath):
data = []
if os.path.exists(filepath):
with open(filepath, mode='r', encoding='gb18030') as f:
lines = csv.reader(f) # 此处读取到的数据是将每行数据当做列表返回的
for line in lines:
data.append(line)
return data
else:
print('filepath is wrong：{}'.format(filepath))
return []
# 读取坐标文件，生成抓取队列
data = read_csv(input_path + '/'+ 'poi_integrate_id.csv') #打开坐标文件即经纬度坐标文件，r是读取的意思
header = data[0]# 记录 header
data = data[1:]# 去掉 header
# 循环coordinates，对每一个coordinates进行fetch，并存储所得数据
n = 10 # 定义要中心坐标点的个数，此处作为示例，仅挑选前10个中心坐标点爬取
for i in range(n):#循环coordinates中的每一个坐标
print('当前处理到了第{0}/{1}个点'.format(i + 1,n))
wgs_x, wgs_y = data[i][1], data[i][2] # 提取坐标，0对应坐标文件中的第1列，1对应坐标文件中的第2列，2对应坐标文件中的第3列
print('wgs_x:{},wgs_y:{}'.format(wgs_x, wgs_y)) # 原始的wgs84坐标系坐标
lng, lat = wgs2bd0911(wgs_x, wgs_y) # lng为BD09ll坐标系下的经度，lat为BD09ll坐标系下的纬度
print('lng:{},lat:{}'.format(lng, lat)) # 转换后的百度坐标系坐标
# 抓取动作
# 根据BaseURL，生成抓取所用URL
initialURL = baseURL + 'ak=' + ak + '&query=' + query + '&scope=' + scope + \
'&location=' + str(lat) + ',' + str(lng) + '&radius=800' + '&page_size=20&page_num=0'#ak,query,scope已经介绍过，不再赘述，location是坐标即经纬度，radius=800即抓取半径是800米,&page_size=20&page_num=0即每页返回20条数据，超过20条数据进行翻页
# 确认此次调用返回数据的页数
response = fetch(initialURL) # 返回抓取结果
# print(response)
totalNum = response['total'] # 返回抓取结果数据的数量
numPages = int(totalNum/20)+1 # 计算页数取整后+1
print (str(numPages) + ' in Total') # 在屏幕上打印页数
# 开始翻页
for i in range(0, numPages+1): #循环每页
print ('Start to fetch page ' + str(i)) #在屏幕上打印抓取到第几页
URL = baseURL + 'ak=' + ak + '&query=' + query + '&scope=' + scope + '&location=' + str(lat) + ',' + str(lng) + '&radius=800' + '&page_size=20&page_num=' + str(i)
response = fetch(URL) #返回抓取url的结果
contents = response['results'] #将返回的结果放在列表contents中
# print(contents) # contents是列表的格式
# 开始循环抓取列表中的所需信息
for content in contents: #content是列表中的内容，是字典格式
# print(content)
name = content['name'] #由于抓取到的数据很多，我们进行进一步的提取；提取poi名字
lat = content['location']['lat'] # 提取纬度
lng = content['location']['lng'] # 提取经度
catagory = content['detail_info']['type'] # 所属分类，如’hotel’、’cater’
province = content['province'] # 所属省份
city = content['city'] # 所属城市
area = content['area'] # 所属区县
# 定义输出的结果为名字+经度+纬度+类别+省+市+区县
poiInfo = name+','+str(lng)+','+str(lat)+','+catagory+','+province+','+city+','+area
poiList.append(poiInfo) # 提取一个poiInfo加到poiList列表中
# poiList.append([name,str(lng),str(lat),catagory,province,city,area]) # 提取一个poiInfo加到poiList列表中
# 最后生成一个csv文件，输出所有抓去结果
# with open(outputFile, mode='w', encoding='gb18030') as f:
# writer = csv.writer(f)
# for i in poiList:
# writer.writerow(i)
# 最后生成一个 txt 文件，输出所有抓取结果。
with open(outputFile, 'w') as f:
for poiInfo in poiList:
f.write(poiInfo + '\n')
f.close() # 出于效率的考虑，只有使用close()函数关闭文件时，才会将缓冲区中的数据真正写入文件中
f = open(outputFile, 'a')
print('poi爬取完毕！')
1.3 小结
本节主要通过构建一个fetch小程序抓取网页中的 url 链接，并对输入的坐标文件中的每一个坐标进行循环、提取目标参数，最终生成 poi 数据文件，完成整个调用 API 爬取百度 POI 数据的过程。
本节我们主要学习到了这些知识：
百度地图开发者平台API的调用； 基于百度地图爬取poi数据的整个基本流程以及poi数据的坐标转换。
在百度地图开放平台中，我们可以获得美食、购物、旅游景点、交通设施等各类 POI 数据。每个类型的 POI 数据对应着一个关键词。例如，我们想爬取公司企业类的 poi 数据，则只需要在上述程序中，将 query = quote ('美食') 中‘美食’改成‘公司企业’即可。其他类型的 POI 数据以此类推。每个类型的关键词在百度开放平台网页中有具体说明，可点击 [百度地图开放平台 poi 分类]进行查看。
[百度地图开放平台 poi 分类]https://lbsyun.baidu.com/index.php?title=lbscloud/poitags
2 兴趣点的处理与向量训练
上一节说明了兴趣点的自动获取方式，接下来介绍利用所获取的兴趣点进行城市土地利用分类的方法。在此之前，需要先将兴趣点数据训练成向量，再利用各个地块内的平均兴趣点向量来对地块进行表征，最后利用各个地块的表征向量完成土地利用分类。poi_integrate_id.csv文件中“ZoneID”字段存储了每个兴趣点所属地块编号。
如前所述，本节大家需要注意：
如何定义兴趣点向量的训练模型和基本参数？ 如何对poi数据和地块数据进行预处理以及进行poi数据的训练？
本节使用镜像为Python 3.7 TensorFlow 2.6 PyTorch 1.8，使用的计算资源是2核8G CPU资源，Kernel类型为Python3。特别注意的是，由于后续需要进行poi的投影坐标转换和距离计算，因此需要额外安装pyproj库，版本为2.4.0（本节后面直接在notebook中用pip导入），如果想一劳永逸，可以先自行在自己的镜像导入pyproj库（自己的镜像务必包含torch库）。此外，训练poi数据可以用CPU也可以用GPU，由于本案例训练的模型不复杂，训练的epoch也不多，因此用CPU也不慢。如果大家设置较大的epoch值，推荐使用GPU训练。
特别注意，由于后面要在程序执行过程中导入pyproj库，因此建议大家逐cell运行，安装好pyproj库后重启核再一键重新运行即可。
2.1 Word2Vec词向量训练模型简介
将兴趣点转化成向量的方法有很多，例如经典的TF-IDF算法、Woed2Vec算法等，本项目以Word2Vec作为基础演示模型。
Word2Vec的方法是在2013年的论文《Efficient Estimation of Word Representations inVector Space》中提出的，作者来自google，
文章下载链接：https://arxiv.org/pdf/1301.3781.pdf
Word2Vec是一种将单词转为向量的方法，其包含两种算法，分别是Skip-gram和CBOW，它们的最大区别是Skip-gram是通过中心词去预测中心词周围的词，而CBOW是通过周围的词去预测中心词。
2.1.1 独热编码简介
词向量顾名思义就是每个单词可以用唯一的一个向量表征，词向量维度大小为整个词汇表的大小。对于每个具体的词汇表中的词，将它们随机排序后，对应的位置置为1，这种词向量的编码方式我们一般叫做独热编码。
比如下面的这个例子，在语料库中，杭州、上海、宁波、北京各对应一个向量，向量中只有一个值为1，其余都为0。杭州 [0,0,0,0,0,0,0,1,0,……0,0,0,0,0,0,0] 上海 [0,0,0,0,1,0,0,0,0,……0,0,0,0,0,0,0] 宁波 [0,0,0,1,0,0,0,0,0,……0,0,0,0,0,0,0] 北京 [0,0,0,0,0,0,0,0,0,……1,0,0,0,0,0,0]
可以看到，独热编码用来表示词向量非常简单，但是却有很多问题。最大的问题是我们的词汇表一般都非常大，比如达到百万级别，这样每个词都用百万维的向量来表示简直是内存的灾难。这样的向量其实除了一个位置是1，其余的位置全部都是0，表达的效率不高，能不能把词向量的维度变小呢？
分布式表征方法可以解决独热编码存在的问题，它的思路是通过训练，将每个词都映射到一个较短（维度较低）的词向量上来。所有的这些词向量就构成了向量空间，进而可以用普通的统计学的方法来研究词与词之间的关系。这个较短的词向量维度是多大呢？这个一般需要我们在训练时自己来指定。一般我们取8的倍数，例如64/128/256维。
2.1.2 Word2Vec模型概述
Word2Vec模型是一种分布式表征方法，它可以将独热编码的向量转化为低维度的连续值，也就是稠密向量，稠密向量中的数字不仅只有0和1，而是有一些其他的浮点数，例如0.9、0.8等。并且其中意思相近的词将被映射到向量空间中相近的位置。如果将embed后的城市向量通过PCA降维后可视化展示出来，那就是这个样子：
可以明显看到，北京、上海和东京的向量的距离较近，因为它们都位于亚洲。而纽约和华盛顿的词向量距离较近，因为它们都位于美国。而北京、上海和东京与纽约和华盛顿相隔较远，是因为它们分别属于不同的洲际。因此，经过Word2Vec表征后的单词，在空间中的距离远近就体现了它们的语义相似性。距离越近，它们的语义就越相似。那Word2Vec的结构是什么样的呢？它复杂吗？
Word2Vec模型其实就是简单化的神经网络。一共有三层：输入层、隐藏层（也叫投影层）和输出层。输入是独热编码后的向量，隐藏层没有激活函数，也就是线性的单元。输出层维度跟输入层的维度一样，用的是Softmax回归。我们要获取的训练好的词向量其实就是隐藏层的输出单元，有的地方定为输入层和隐藏层之间的权重，其实说的是一回事。
Word2Vec模型分为CBOW(Continuous Bag-of-Words)模型和Skip-gram模型。CBOW模型根据某个中心词前后n个连续的词，来计算该中心词出现的概率，即用上下文预测目标词。Skip-gram只是逆转了CBOW的因果关系而已，即已知当前词语，预测上下文。两种模型结构简易示意图如下：
本例选择Skip-gram模型训练兴趣点向量，因此此处重点介绍Skip-gram模型。我们将训练神经网络做以下工作：给定一个句子中间的一个指定单词（输入单词），查看附近的单词并随机选择一个。网络将告诉我们，词汇中每个单词是我们选择的"附近单词"的概率。“附近”指的是算法中的一个“窗口大小”参数，典型的窗口大小可能是5，意思是后面5个字，前面5个字（总共10个）。
输出概率将与在输入单词附近找到词汇表中每个单词的可能性有关。例如，如果我们把单词“Soviet”作为训练后模型的输入，那么单词，如“Union”和“Russia”，的概率将远远高于不相关单词，如“watermelon”和“kangaroo”。
我们将通过给神经网络提供从训练文档中找到的单词对来训练神经网络。下面的例子展示了一些来自句子“The quick brown fox jumps over the lazy dog”的训练样本（单词对）。在这个例子中，我使用了一个较小的窗口尺寸（2），蓝色高亮的单词是输入单词，右边的Training Samples表示输入网络的训练词对（pairs）。
这个网络将要从每个对出现的次数中学习统计量。因此，例如，这个网络可能会得到比(“Soviet”, “Sasquatch”)更多的训练样本(“Soviet”, “Union”)。当训练结束后，如果我们将单词“Soviet”作为输入，随后它将给“Union”和“Russia”输出比“Sasquatch”更高的概率。
以上是Word2Vec训练单词的模型和方法，而兴趣点和单词有所不同，兴趣点是空间实体，具有位置信息（经纬度），而单词则简单一些。那如何将单词层面的Word2Vec方法，迁移到空间层面的兴趣点方面呢？
2.1.3 由Word2Vec模型到Place2Vec模型
为了演示方便，本例最终以Place2Vec模型作为兴趣点向量训练模型。更多关于Place2Vec模型，可以参照其原始论文：[From ITDL to Place2Vec – Reasoning About Place Type Similarity and Relatedness by Learning Embeddings From Augmented Spatial Contexts]。
https://www.researchgate.net/publication/320778966_From_ITDL_to_Place2Vec_--_Reasoning_About_Place_Type_Similarity_and_Relatedness_by_Learning_Embeddings_From_Augmented_Spatial_Contexts
Place2Vec模型与Skip-gram模型区别不大，简而言之，Place2Vec模型是基于Skip-gram模型的。他们最主要的区别是，Skip-gram模型的输入是根据句子中词语的顺序来构造pairs的，而Place2Vec模型则是根据兴趣点在空间中的距离来构造pairs的。例如，假设某中心词为某餐饮店，取与它在空间中距离最近的10个兴趣点（数量可以随意定）作为周围词来构造兴趣点pairs，那么一共就可以构造10对pairs作为输入。那么这样利用Place2Vec模型便可以将Skip-gram模型由单词层面迁移到兴趣点层面。
为了便于说明，后文"利用Word2Vec模型训练的向量"和"利用Place2Vec模型训练的向量"指代的是同一个意思，因为Place2Vec模型本来就是在Word2Vec模型的基础上改进得到的，本质还是Word2Vec模型模型。
下面是具体的利用Place2Vec模型训练兴趣点向量的方法：
2.2 定义Skip-gram模型函数
在定义函数之前，由于此镜像不存在pyproj包，因此需要先下载一个pyproj包，用于后续的距离计算及投影坐标转换。执行pip安装后需要重启核并重新运行才能生效。
In [1]:
!pip install pyproj==2.4.0 -i https://pypi.douban.com/simple/ # pip安装，从指定镜像下载安装工具包，镜像URL可自行修改
Looking in indexes: https://pypi.douban.com/simple/
Requirement already satisfied: pyproj==2.4.0 in /opt/conda/lib/python3.7/site-packages (2.4.0)
In [2]:
"""导入相关的包"""
import numpy as np
import torch.nn as nn # torch核心模块，包括常用神经网络及损失函数等，方便搭建模型
import torch.optim as optim # 该模块定义了不同梯度下降优化器方法
import torch.utils.data as Data # torch.utils是辅助模型训练、测试和结构优化的模块
from scipy.spatial.distance import cdist # 该函数用于计算两个输入集合的距离
import os
import torch
import torch.nn.functional as F
import torch.utils.data as tud # torch.utils.data引入数据集（Dataset）和数据载入器（DataLoader）方便对训练中数据进行操作
from collections import Counter # Counter 是 dictionary 对象的子类。collections 模块中的 Counter() 函数会接收一个诸如 list 或 tuple 的迭代器，然后返回一个 Counter dictionary
import random
import pandas as pd
from pyproj import CRS, Transformer # CRS初始坐标系，Transformer进行坐标转换
In [3]:
# 定义Word2Vec模型的Skip-grim方法的函数
# pairs为点对
def skip_gram(vocabulary, word_pairs, word_to_idx, id_to_word, k):
"""
学习词向量的概念
用Skip-thought模型训练词向量
学习使用PyTorch dataset和dataloader
学习定义PyTorch模型
学习torch.nn中常见的Module
Embedding
学习常见的PyTorch operations
bmm
logsigmoid
保存和读取PyTorch模型
在这一份notebook中，我们会（尽可能）尝试复现论文Distributed Representations of Words and
Phrases and their Compositionality中训练词向量的方法. 我们会实现Skip-gram模型，并且使用论文中noice contrastive sampling的目标函数。
这篇论文有很多模型实现的细节，这些细节对于词向量的好坏至关重要。
我们虽然无法完全复现论文中的实验结果，主要是由于计算资源等各种细节原因，但是我们还是可以大致展示如何训练词向量。
以下是一些我们没有实现的细节
subsampling：参考论文section 2.3
"""
# 设备选择
device = 'cuda' if torch.cuda.is_available() else 'cpu'
USE_CUDA = torch.cuda.is_available()
# 为了保证实验结果的可以复现，我们经常会把各种random的seed固定在某一个值
random.seed(53113)
np.random.seed(53113)
torch.manual_seed(53113)
if USE_CUDA:
torch.cuda.manual_seed(53113)
# 设定一些超参数
K = 100 # number of negative samples 负采样单词数量
NUM_EPOCHS = 3 # 训练周期。模型将使用的训练集全部样本训练一遍即为一个epoch，设置3即模型将训练集重复训练3次
BATCH_SIZE = 128 # 批大小。小批量梯度下降算法中，设置的一次载入内存进行批梯度计算的样本数量
LEARNING_RATE = 0.2 # 学习率。梯度下降算法的步长或者收敛速度。过小收敛缓慢，过大则可能无法收敛，或者在最值附近震荡
EMBEDDING_SIZE = 256 # 训练的向量的维度
# 日志文件
LOG_FILE = 'word-embedding.log'
word_counts = np.array([count for count in vocabulary.values()], dtype=np.float32) # 每个单词出现次数的列表
"""
单词出现频率 单词出现次数/单词总的个数
把频次乘以3/4效果会好一些 这是论文里面的trick
word_freqs 用来做负采样
"""
word_freqs = word_counts / np.sum(word_counts)
word_freqs = word_freqs ** (3. / 4.)
word_freqs = word_freqs / np.sum(word_freqs)
# print(word_freqs)
VOCAB_SIZE = len(id_to_word)
"""
实现Dataloader
一个dataloader需要以下内容：
把所有text编码成数字，然后用subsampling预处理这些文字。
保存vocabulary，单词count，normalized word frequency
每个iteration sample一个中心词
根据当前的中心词返回context单词
根据中心词sample一些negative单词
返回单词的counts
这里有一个好的tutorial介绍如何使用PyTorch dataloader. 为了使用dataloader，我们需要定义以下两个function:
__len__ function需要返回整个数据集中有多少个item
__get__ 根据给定的index返回一个item
有了dataloader之后，我们可以轻松随机打乱整个数据集，拿到一个batch的数据等等。
"""
class WordEmbeddingDataset(tud.Dataset):
"""
text 单词文本 来自与train 训练数据
word_to_idx 单词到索引的映射 字典
idx_to_word 索引到单词的列表 下标即是索引
word_freqs 每个单词的频率 出现次数/总次数 ××（3/4）列表
word_counts 每个单词的出现次数
"""
def __init__(self, pair, word2idx, idx2word, freqs, counts):
super(WordEmbeddingDataset, self).__init__()
"""
将text转为数字索引编码的列表
并转为tensor 张量
"""
self.text_encoded = pair # 已经得到的（中心词，周围词）点对
self.text_encoded = torch.Tensor(self.text_encoded).long()
self.word_to_idx = word2idx
self.idx_to_word = idx2word
self.word_freqs = torch.Tensor(freqs)
self.word_counts = torch.Tensor(counts)
def __len__(self):
return len(self.text_encoded)
def __getitem__(self, idx):
"""
首先 skgram 这个模型的任务就是 根据中心词 预测周围的词
这个function 返回一下数据item 用于训练
---中心词 当前的词
---这个单词附近的（上下文正确的单词)的单词
---随机采样的K个单词作为负样本 就是负样本采样
"""
# ----中心词
center_word = self.text_encoded[idx][0]
# ---中心词窗口的附近正确的单词,共有k个
pos_words = self.text_encoded[idx][1:]
"""
负采样的单词 以及负采样的规则
word_freqs 每个单词的频率 出现次数/总次数 ××（3/4）列表
word_freqs 用来做负采样
torch.multinomial(input, num_samples,replacement=False, out=None) → LongTensor
作用是对input的每一行做n_samples次取值，输出的张量是每一次取值时input张量对应行的下标。
输入是一个input张量，一个取样数量，和一个布尔值replacement。
input张量可以看成一个权重张量，每一个元素代表其在该行中的权重。如果有元素为0，那么在其他不为0的元素
被取干净之前，这个元素是不会被取到的。
n_samples是每一行的取值次数，该值不能大于每一样的元素数，否则会报错。
replacement指的是取样时是否是有放回的取样，True是有放回，False无放回。
按权重 取得是对应的下标 在这里就是按着每个单词的频率作为权重 每次取出K=100*6个单词 的索引 也就是单词数字
编码 这就是负采样的规则 按照单词出现频率进行采样 出现频率越高 被采样的次数概率越高
并且是有放回的
"""
negative_places = list(
set(range(len(self.word_to_idx))) - set(pos_words))
neg_words = torch.multinomial(self.word_freqs[negative_places], K * pos_words.shape[0], True)
# random.sample从指定序列中随机获取指定长度的片断
# negative_context = random.sample(negative_places, int(len(self.second_class_walks[index]) * k))
return center_word, pos_words, neg_words
"""
创建dataset 和dataloader并且测试 返回一个样本 中心词 中心词附近的正确单词 负采样出来的负样本
"""
dataset = WordEmbeddingDataset(word_pairs, word_to_idx, id_to_word, word_freqs, word_counts)
dataloader = tud.DataLoader(dataset, batch_size=BATCH_SIZE, shuffle=True, num_workers=0)
"""
定义pytorch模型
"""
class EmbeddingModel(nn.Module):
def __init__(self, vocab_size, embed_size):
# 初始化输入和输出embedding
super(EmbeddingModel, self).__init__()
# 词典大小
self.vocab_size = vocab_size
# 词嵌入维度
self.embed_size = embed_size
# 初始化的权重范围 可以更好的拟合模型
initrange = 0.5 / self.embed_size
# 嵌入层
self.out_embed = nn.Embedding(self.vocab_size, self.embed_size, sparse=False)
# 初始化输出层嵌入层权重
self.out_embed.weight.data.uniform_(-initrange, initrange)
# 输入嵌入层
self.in_embed = nn.Embedding(self.vocab_size, self.embed_size, sparse=False)
# 初始化输入嵌入层权重
self.in_embed.weight.data.uniform_(-initrange, initrange)
def forward(self, input_labels, pos_labels, neg_labels):
"""
input_labels:中心词 [batch_size] [1,2,..........128..............] 输入是数字表示
pos_labels:中心词周围context window 出现过的单词[batch_size*window_size*2] [128*6,]
neg_labels: 中心词周围没有出现的单词 负采样过来的单词 [batch_size*window_size*2*K] >>[128*6*100,]
return: loss, [batch_size]
"""
batch_size = input_labels.size(0)
"""
#输入嵌入层 把输入单词词典维度嵌入到embed_size维的空间
#数据维度变化[batch_size,onehot_vocab_size]>>>>[batch_size,embed_size]
# 参数权重的维度 [vocab_size,embed_size] 这就是我们输入的词嵌入 词向量
#每一行就可以表示为这个单词的 embed_size维度的向量
input_embedding= [batch_size,embed_size]
"""
input_embedding = self.in_embed(input_labels)
"""
输出层词嵌入层 和上面词嵌入网络层原理一样
#数据维度变化[batch_size,2*c,onehot_vocab_size]>>>>[batch_size,embed_size]
# 参数权重的维度 [vocab_size,embed_size] 这就是我们输入的词嵌入 词向量
#每一行就可以表示为这个单词的 embed_size维度的向量
pos_embedding=[batch_size,2*C,embed_size]
"""
pos_embedding = self.out_embed(pos_labels)
"""
负采样单词 和预测输出单词共享权重变量
#数据维度变化[batch_size,onehot_vocab_size]>>>>[batch_size,embed_size]
# 参数权重的维度 [vocab_size*2*C×K,embed_size] 这就是我们输入的词嵌入 词向量
#每一行就可以表示为这个单词的 embed_size维度的向量
neg_embedding= [batch_size,2*C*K,embed_size]
"""
neg_embedding = self.out_embed(neg_labels)
"""
计算损失函数
这里的目标函数是重点
简单粗暴的理解下目标函数
这个模型就做了一件事，就是利用input嵌入层得到输入词向量
利用输出词向量得到中心词周围的词的词向量和负采样得到的 词向量
》》》》》》》》》》》》》》》》》》》》》》》》》》》》》》》》》》》》》》》》
那么假设中心词跟周围词关系越近，跟负采样的词关系很远。这个假设很关键，是个隐藏的条件
那么目标函数为：中心词和周围词进行矩阵乘法，值越大越好。中心词与负采样的词矩阵乘法，值越小越好，取负-，也是越大越好
这个模型非常的简单，只有嵌入层，没有任何的隐藏层，就是一个表示模型
"""
"""
具体的目标函数 会在文章末尾给出
第一步： 正确的 输出词向量 点积 输入词向量
pos_embedding=[batch_size,2*C,embed_size] input_embedding= [batch_size,embed_size]
化成相同的形状
pos_embedding=[batch_size,2*C,embed_size] input_embedding= [batch_size,embed_size,1]
得到
log_pos=[batch_size,2*C,1]>>squeeze()>>>[batch_Size,2*C]
第二步 负采样的 输出词向量 点积 负的 -输入词向量
neg_embedding=[batch_size,2*C×K,embed_size] input_embedding= [batch_size,embed_size]
化成相同的形状
pos_embedding=[batch_size,2*C*K,embed_size] input_embedding= [batch_size,embed_size,1]
得到
log_neg=[batch_size,2*C*K,1]>>squeeze()>>>[batch_Size,2*C*K]
第三步： 将log_pos 进行logsigmoid 并在第二维度上进行求和
第四步： log_neg 进行logsigmoid 并在第二维度上进行求和
第五步 将2者相加 [batch_Size]
"""
# 计算的是与周围词的相似度 越大越好
log_pos = torch.bmm(pos_embedding, input_embedding.unsqueeze(2)).squeeze()
# 计算的是负采样的词的相似度 取- 越大越好
log_neg = torch.bmm(neg_embedding, -input_embedding.unsqueeze(2)).squeeze()
log_pos = F.logsigmoid(log_pos).sum(1)
log_neg = F.logsigmoid(log_neg).sum(1) # batch_size
loss = log_pos + log_neg
# 取- 那么需要这个值最小
return -loss
def input_embeddings(self):
return self.in_embed.weight.data.cpu().numpy()
# 定义一个模型，以及如果有gpu，就把模型移动到gpu进行计算
model = EmbeddingModel(VOCAB_SIZE, EMBEDDING_SIZE)
if USE_CUDA:
model = model.cuda()
"""
训练模型：
模型一般需要训练若干个epoch
每个epoch我们都把所有的数据分成若干个batch
把每个batch的输入和输出都包装成cuda tensor
forward pass，通过输入的句子预测每个单词的下一个单词
用模型的预测和正确的下一个单词计算cross entropy loss
清空模型当前gradient
backward pass
更新模型参数
每隔一定的iteration输出模型在当前iteration的loss，以及在验证数据集上做模型的评估
"""
optimizer = torch.optim.SGD(model.parameters(), lr=LEARNING_RATE)
for e in range(NUM_EPOCHS):
for i, (input_labels, pos_labels, neg_labels) in enumerate(dataloader):
# todo
input_labels = input_labels.long()
pos_labels = pos_labels.long()
neg_labels = neg_labels.long()
if USE_CUDA:
nput_labels = input_labels.cuda()
pos_labels = pos_labels.cuda()
neg_labels = neg_labels.cuda()
optimizer.zero_grad()
# loss越小 模型越准确
loss = model(input_labels, pos_labels, neg_labels).mean()
loss.backward()
optimizer.step()
if i % 100 == 0:
with open(LOG_FILE, "a") as fout:
fout.write("epoch: {}, iter: {}, loss: {}\n".format(e, i, loss.item()))
print("epoch: {}, iter: {}, loss: {}".format(e, i, loss.item()))
# 所有模型训练完保存一次
embedding_weights = model.input_embeddings()
state = {
'embedding': embedding_weights,
} # 建立一个字典，保存训练好的embedding
torch.save(state, '/home/mw/project/embedding_k_is_{}_epoch_{}.pth'.format(k, NUM_EPOCHS)) # 保存最后一轮的参数
return embedding_weights
2.3 训练poi向量
定义好Skip-gram函数后，接下来就可以正式开始训练兴趣点向量了。训练的思路是：
首先加载兴趣点数据和深圳市的地块数据，兴趣点数据中包含深圳市2013年和2020年的百度地图兴趣点，深圳市地块数据包含了深圳市每个地块所属的具体编号，深圳市一共有6912个地块，“ZoneID”值为-1代表该poi点不属于任何地块，即不在任何地块范围内。 找出有poi的地块（poi_integrate_id.csv文件中不为-1的所有其它“ZoneID”值），逐地块进行循环。对于每个地块，先找出其包含的poi，再针对每个poi，求与其距离最近的k个poi。若地块内包含的poi个数N=1，则与其最近的poi均用此poi代替。若1 < N <= k，则不足的k-N个poi用第一个最近的poi代替。
In [4]:
integrate_path = u'/home/mw/input/Data5799/数据集/poi_integrate_id.csv' # 2013与2020年集成poi的文件存储路径
pth_parcel = u'/home/mw/input/Data5799/数据集/shenzhen_parcel.csv' # 地块文件路径
Years = [2013, 2020] # poi数据的年份
# 用于加载数据，path为文件路径
def load_data(path):
original_poi = pd.read_csv(path, encoding='gb18030') # 读取原始poi数据
return original_poi
# 如果要进行距离计算，需要先进行投影，因为经纬度是球面坐标，不方便直接计算距离，需要先将其转换为平面坐标才方便计算距离。
# 该函数用于将wgs84地理坐标转换为对应投影的平面直角坐标
# wgs_x为经度坐标，wgs_y为纬度坐标
def transform(wgs_x, wgs_y):
p1 = CRS.from_epsg(4326) # 定义数据地理坐标系，wgs_84地理坐标系的EPSG Code
p2 = CRS.from_epsg(32649) # WGS_1984_UTM_Zone_49N投影坐标系的EPSG Code
transformer = Transformer.from_crs("epsg:4326", "epsg:32649") # 构造转换对象
x, y = transformer.transform(np.array(wgs_y), np.array(wgs_x)) # 参数为元组形式，不能为list。x，y为投影坐标
return x, y
# 加载数据
parcel = load_data(pth_parcel)
integrate_corpus = load_data(integrate_path) # 集成的语料库
FirstLevel = integrate_corpus['FirstLevel'].drop_duplicates().values.tolist() # 第一大类，用于可视化，drop_duplicates作用是去重
zone_all = integrate_corpus.loc[integrate_corpus['ZoneID'] > -1] # 有poi的地块编号都不为-1，该变量表示满足此条件的所有poi数据
k = 10 # 取k个最邻近的poi作为周围词的数量
pairs = [] # 初始化包含点对的列表
# 开始按年份循环地块
for temp_year in Years:
if temp_year == 2013: # 作业2中需修改此处的2013为2020
temp_year_poi = integrate_corpus.loc[integrate_corpus['Year'] == temp_year] # 所有当前年份的poi数据
# 构造字典
sequence = temp_year_poi['SecondLevel'].values.tolist() # 当前年份所有的poi二级类，有重复
vocab = dict(Counter(sequence)) # 统计每个类别出现的次数，字典形式
idx_to_word = [word for word in vocab.keys()]
word2id = {word: i for i, word in enumerate(idx_to_word)} # poi词典，每个poi名称对应一个id序号
print('当前年份的poi数据词典长度为：', len(word2id))
poi_with_zone = temp_year_poi.loc[temp_year_poi['ZoneID'] > -1] # 当前年份落在地块内的poi，即地块编号>-1的poi
temp_year_zone = poi_with_zone['ZoneID'].drop_duplicates().values.tolist() # 有poi的当前年份的地块列表
count = 0
# 外层循环为地块的顺序
for order in temp_year_zone: # 设置地块循环是为了方便记录执行到哪些点了，就是提示进行到哪一步了
# 统计并记录每个地块包含的poi
poi_class = poi_with_zone.loc[poi_with_zone['ZoneID'] == order, 'SecondLevel'].values.tolist() # poi类别
poi_class = np.array(poi_class) # 将列表转为数组，便于将编号转为类别文字
poi_Lon = poi_with_zone.loc[poi_with_zone['ZoneID'] == order, 'Lon'].values.tolist() # 经度
poi_Lat = poi_with_zone.loc[poi_with_zone['ZoneID'] == order, 'Lat'].values.tolist() # 纬度
# 坐标转换，同一年数据可以直接转换
poi_x, poi_y = transform(poi_Lon, poi_Lat) # wgs84地理坐标转换为平面坐标
poi_xy = np.array(list(zip(poi_x, poi_y))) # 当前地块所有poi的经纬度数组
N = poi_xy.shape[0] # 当前地块poi的总个数
'''构造（中心poi，周围poi）点对，每个poi构造10个对'''
# 如果地块只有一个poi，则复制后直接添加
if N == 1:
pairs.append([]) # 创建二维列表，第一个存储中心词，后面的是k个周围词
center = word2id[poi_class[0]] # center word
context = [center] * k
pairs[-1].append(center) # 添加中心词
for w in context:
pairs[-1].append(w)
else:
dist_all = cdist(poi_xy, poi_xy) # 求每对点间的距离
row, col = np.diag_indices_from(dist_all)
dist_all[row, col] = np.zeros((dist_all.shape[0],)) # 替换对角值为0，自己和自己距离为0
sorted_index = np.argsort(dist_all) # 对距离矩阵每一行按从小到大排序，返回array
sorted_index = sorted_index[:, 1:] # 排除第一列0
# 对每一行求点对
for idx in range(sorted_index.shape[0]): # 每一行索引代表一个中心词id
pairs.append([])
center = word2id[poi_class[idx]] # center word
pairs[-1].append(center) # 末尾添加中心词
if 1 < N <= k: # 有几个取前几个，并补齐至k个
for i in range(k - N + 1):
context = word2id[poi_class[sorted_index[idx, 0]]] # 以第一个最近的poi补齐
pairs[-1].append(context) # 周围词列表，第0个不取
n = N - 1
else: # 取前k个
n = k
for i in range(n):
context = word2id[poi_class[sorted_index[idx, i]]]
pairs[-1].append(context) # 周围词列表，第0个不取
count += 1
# \r 表示将光标的位置回退到本行的开头位置,配合end=''不换行，实现打印新内容时删除旧内容
print('\r已构造完{0}年的第{1}个地块的点对'.format(temp_year, count), end="")
else:
continue
state = {
'pairs': pairs,
} # 建立一个字典，保存训练好的embedding
torch.save(state, '/home/mw/project/pairs_k_is_{}.pth'.format(k)) # 保存训练的点对
# pairs = state['pairs']
# count = len(pairs) # 查看兴趣点pairs的数量
embedding = skip_gram(vocab, pairs, word2id, idx_to_word, k) # 返回词向量
print('词向量训练完成')
2.4 小结
本节介绍了关于Word2Vec、Place2Vec等模型的基本知识，并用实例代码演示了深圳市2013年兴趣点的训练过程和方法，最终得到了深圳市2013年的兴趣点向量。
本节我们主要学习到了这些知识：
词向量、独热编码的基本概念，以及Word2Vec、Place2Vec模型的基本结构； poi数据的处理以及利用Word2Vec/Place2Vec模型训练poi向量的基本流程和方法。
下一节，我们将学习聚类和可视化的有关知识，以及如何将这些刚训练的兴趣点向量可视化地展示出来。此外，我们还将介绍如何得到地块向量，以及利用地块向量进行聚类得到土地利用分类数据的方法。
3 兴趣点向量的定性精度评价与土地利用分类
本节使用镜像为Python 3.7 TensorFlow 2.6 PyTorch 1.8，使用的计算资源是2核8G CPU资源，Kernel类型为Python3。由于需要导入训练的向量文件（.pth）以及做聚类分析，因此推荐大家使用包含torch库和Scikit-learn库的镜像。但向量的聚类及可视化不涉及gpu的使用，因此大家使用CPU资源就可以了。
上一节只得到了兴趣点向量，但并不知道所训练的兴趣点向量的效果如何。那如何评价所训练的兴趣点向量的质量效果呢？ 此外，得到了poi向量后，如何进一步获得地块向量及利用它进行土地利用分类呢？
首先，可以仿照上一节中Word2Vec模型中关于“北京、上海和东京”等词的词向量可视化来定性评价所训练的兴趣点向量。而由于我们的向量数量较多，不方便直接在图中将文字和点同时标注在一起。因此在可视化之前，我们首先对向量进行聚类，这样可视化时在同一个聚类簇里面的向量距离相近，具有较高的相似性。再对不同的聚类簇予以不同的颜色显示，就能很直观地看到聚类的效果了。那什么是聚类分析呢？
特别注意：
由于聚类前需要先计算Silhouette值以确定最佳的聚类数，因此建议大家在做作业时逐个cell依次运行，根据silhouette值调整聚类数。而不要一键运行，一键运行所有的cell在原参数没有改变的情况下可能无法得到最好的结果。 程序运行到visualize函数时会出现FutureWarning警告，这是包的版本问题，不影响程序运行，因此可以忽略这里的警告。
3.1 K-Means聚类算法简介
聚类分析就是将大量数据中具有"相似"特征的数据点或样本划分为一个类别，聚类分析提供了样本集在非监督模式下的类别划分。聚类的基本思想是"物以类聚、人以群分"，将大量数据集中相似的数据样本区分出来，并发现不同类的特征。聚类模型可以建立在无类标记的数据上，是一种非监督的学习算法。聚类根据数据自身的距离或相似度将他们划分为若干组，就是给样本打上不同类的标签，划分原则是组内样本最小化而组间距离最大化。常用的聚类算法有K-Means聚类算法、DBSCAN（具有噪声的密度聚类算法）等，常见的聚类算法及聚类效果可见下图：
本例中我们选择K-Means聚类算法进行聚类。那什么是K-Means聚类算法呢？
简而言之，K-Means聚类的基本思想就是：对于给定的样本集，按照样本之间的距离(也就是相似程度)大小，将样本集划分为K个簇(即类别)。让簇内的点尽量紧密的连在一起，而让簇间的距离尽量的大。
算法步骤：步骤1：随机取k个初始中心点 步骤2：对于每个样本点计算到这k个中心点的距离，将样本点归到与之距离最小的那个中心点的簇。这样每个样本都有自己的簇了 步骤3：对于每个簇，根据里面的所有样本点重新计算得到一个新的中心点，如果中心点发生变化回到步骤2，未发生变化转到步骤4 步骤4：得出结果
就像这样：
知道了K-Means聚类的基本原理，那该如何实现它呢？
In [29]:
from sklearn.cluster import KMeans # 导入聚类所需的包
In [30]:
# K-means聚类,k为聚类数，X为向量数据
# 使用欧氏距离作为距离度量标准
# 返回聚类结果
def kmeans(k, X):
# 正式定义模型
# model1 = KMeans(n_clusters=k, n_init=10) # n_init：用不同的初始化质心运行算法的次数。
model1 = KMeans(n_clusters=k, random_state=500) # random_state：用于初始化质心的生成器（generator）,设为固定值保证结果可复现。
# 跑模型
model1.fit(X)
# 需要知道每个类别有哪些参数
cluster_result = model1.predict(X) # 每个向量所属的类别,[1,1,2,0,3...]这种形式，类别从0开始
return cluster_result
3.2 确定聚类数k
通过上述程序可以看到，K-Means聚类的聚类数k是一个超参数，是需要事先给定的。那么该如何确定聚类数k的值呢？
这里我们选择利用Silhouette值（轮廓系数）来决定k值。Silhouette值是一个聚类的评价指标，用来描述一个目标对于目标所在簇与其他簇之间的相似性。其范围是从-1~+1，这个值越大表明目标与自己所在簇之间的匹配关系度越高，与其他簇的匹配关系度越低。如果这个值越高，那么聚类结果越好，如果是很小或是负值，那么可能是分簇太多或是太少造成的。
因此，为了确定具体的k值，我们选择将k值取[2,26]之间的整数，分别比较取各个k时的轮廓系数，选择其中轮廓系数相对较高的k值作为我们最终的聚类数。具体的计算程序及轮廓系数可视化的程序如下：
In [31]:
# 先导入计算轮廓系数相关的包和绘图可视化包
from sklearn.metrics import silhouette_score # score计算所有样本的平均轮廓系数
import matplotlib.pyplot as plt
In [32]:
# 计算各个分类的平均轮廓系数(silhouette值),X为向量,m为不同的聚类数2~m构成一个聚类区间
# m的作用相当于对k分别取值2~m之间的整数进行聚类，找出其中Silhouette值较大的一个k值
def compute_silhouette(X, m):
# 依次对不同的k值进行聚类,k取值为2~m之间的整数（不包括m）
# n_int指定了k均值算法运行的次数。每一次都会选择一组不同的初始化均值向量，最终算法会选择最佳的分类簇来作为最终的结果。
kmeans_per_k = [KMeans(n_clusters=k, random_state=500).fit(X) for k in range(1, m)]
# 必须至少有两类才能计算silhouette值，所以从1开始
silhouette_scores = [silhouette_score(X, model.labels_) for model in kmeans_per_k[1:]]
# 画图
# plt.figure(figsize=(8, 3)) # figure的宽和高，单位为英寸
plt.plot(list(range(2, m)), silhouette_scores, color='red', marker='o')
plt.xlabel("k", fontdict={'size': 16})
plt.ylabel("Silhouette score", fontdict={'size': 16})
plt.xticks(range(2, m, 1)) # 设置x轴坐标间隔
plt.show()
3.3 向量降维及可视化
确定了聚类数k以及利用k值进行聚类后，由于向量的维度较高，达到256维，因此不能直接在二维平面进行显示，需要先将其降至2维才方便显示。那么该如何对向量进行降维呢？
常用的降维方法有PCA（主成分分析）、T-SNE（t-分布领域嵌入算法）方法等。这里我们选择T-SNE方法进行降维。它的主要想法就是：将高维分布点的距离，用条件概率来表示相似性，同时低维分布的点也这样表示。只要二者的条件概率非常接近（用相对熵来训练，所以需要label），那就说明高维分布的点已经映射到低维分布上了。
具体的数学原理这里不做过多介绍，感兴趣的同学可以下去具体了解。这里主要介绍利用T-SNE方法对向量进行降维及向量降维后的可视化方法，具体程序如下：
In [33]:
import numpy as np
from sklearn.manifold import TSNE # 导入T-SNE方法的包
In [34]:
# 定义向量降维可视化的函数
# embedding为需要可视化的向量，category为聚类类别
def visualize(embedding, category):
"""T-SNE降维"""
# 对原始向量降维
X_tsne = TSNE(n_components=2, init='pca', random_state=33).fit_transform(embedding) # 降至二维
x_min, x_max = X_tsne.min(0), X_tsne.max(0) # min(0)返回该矩阵中每一列的最小值，max(0)返回该矩阵中每一列的最大值
X_norm = (X_tsne - x_min) / (x_max - x_min)
# 计算质心坐标
centroids_xy = [] # 初始化每一类簇的质心坐标
for temp_class in range(np.max(category) + 1):
temp_poi_id = np.argwhere(category == temp_class) # 为二维数组，shape只有一列
temp_poi_id = temp_poi_id.flatten() # 展平为一维数组
temp_class_embedding = X_norm[temp_poi_id.tolist()] # 找出当前类的poi嵌入向量，二维数组
temp_centroid_xy = np.mean(temp_class_embedding, axis=0) # 找出当前类点的质心坐标，一维数组
centroids_xy.extend([temp_centroid_xy.tolist()])
# 画图
plt.figure(figsize=(12, 12)) # 设置画布大小,可以根据需要调整
plt.scatter(X_norm[:, 0], X_norm[:, 1], c=category, cmap=plt.cm.rainbow)
# 添加注释
cluster_id = ['C' + str(i) for i in range(np.max(category) + 1)] # 聚类名列表，用于在每个类簇的质心指示具体的类簇是哪个
for i in range(len(cluster_id)):
plt.text(centroids_xy[i][0], centroids_xy[i][1], cluster_id[i], fontsize=16,
color="black",
style="normal",
weight="bold", verticalalignment='center',
horizontalalignment='right')
# 设置坐标刻度值的大小
plt.yticks(size=16)
plt.xticks(size=16)
plt.show()
3.4 poi向量聚类可视化
In [35]:
# 定义好聚类及可视化的函数以后，接下来就可以开始主程序的撰写了
# 先导入相关包
import torch
import pandas as pd
from collections import Counter
In [36]:
place2vec_path='/home/mw/project/embedding_k_is_10_epoch_3.pth' # 上一节训练的向量的保存路径
# 最大最小值归一化处理，便于绘图显示，规范坐标
def normalization(X):
# 特征缩放
X_min, X_max = np.min(X), np.max(X)
X = (X - X_min) / (X_max - X_min)
# 根据数据需要进行不同的处理
return X
# 获取向量
state = torch.load(place2vec_path)
poi_embedding = state['embedding'] # place2vec方法得到的原始向量
poi_embedding = normalization(poi_embedding) # 最大最小值归一化
k1 = 26 # 初始向量聚类数k的试验区间为[2,26],k为整数
compute_silhouette(poi_embedding, k1) # 计算每个k值下的平均silhouette值，确定最优k值
如上图可见，选取平均Silhouette值曲线图中较大Silhouette值对应的k值作为下面聚类的聚类数，即18。因此，如下方程序所示，kmeans函数的第一个参数设置为18。
特别注意，在做作业时这里需要大家根据具体的数据以及具体的平均Silhouette值曲线图选择合适的聚类数k，因为2013年和2020年的poi向量存在差异，它们的最佳聚类数也不一定相同。
确定了聚类数以后便可以继续下面对于poi向量的kmeans聚类以及可视化了。
In [37]:
cluster_category = kmeans(18, poi_embedding) # 使用K-Means聚类
print(cluster_category) # 查看聚类结果
print(cluster_category.shape) # 查看向量维度
visualize(poi_embedding, cluster_category) # 可视化poi聚类向量
上图可以根据颜色明显看出类簇呈现相对聚集的分布，向量化效果还不错。定性评价向量化效果没有特别固定的标准，若相同颜色的点聚在一起地越紧凑，并与不同颜色的类簇之间有明显的区别，则聚类效果越好，向量化的效果也越好。
3.5 地块向量的获取及聚类可视化输出
In [38]:
# 获取d年份地块向量，最终得到每个地块的变化向量
# embedding为通过Word2Vec模型获取的poi嵌入向量
def parcel_embedding(embedding, integrated_poi, year):
temp_year_poi = integrated_poi.loc[integrated_poi['Year'] == year]
zone = temp_year_poi.loc[
integrated_poi['ZoneID'] > -1, 'ZoneID'].drop_duplicates().values.tolist() # 有poi的地块编号都不为-1
# 构造字典
sequence = temp_year_poi['SecondLevel'].values.tolist() # 当前年份所有的poi二级类，有重复
vocab = dict(Counter(sequence)) # 统计每个类别出现的次数，字典形式
idx_to_word = [word for word in vocab.keys()]
word2id = {word: i for i, word in enumerate(idx_to_word)} # poi词典，每个poi名称对应一个id序号
count = 0
parcel_mean = np.zeros(shape=(0, embedding.shape[1])) # 初始化地块向量的空数组，用于拼接
parcel_contain_all_poi = [] # 初始化包含poi的地块列表
# 外层循环为地块的顺序
for order in zone:
parcel_contain_all_poi.append(order) # 添加当前地块编号
# 统计并记录每个地块包含的poi
class_temp_year = temp_year_poi.loc[temp_year_poi['ZoneID'] == order, 'SecondLevel'].values.tolist() # 当前地块2013年poi类别
id_temp_year = [word2id[level] for level in class_temp_year] # 根据poi类别获取其在词典中的id
'''求地块平均向量'''
# 首先获取向量
embedding_temp_year = embedding[id_temp_year]
# axis=0表示对每行求平均
mean_temp_year = embedding_temp_year.mean(axis=0)
# 向量拼接,需要对embedding_change扩充一个维度才能拼接
parcel_mean = np.concatenate((parcel_mean, mean_temp_year[np.newaxis, :]), axis=0) # 纵向拼接
count += 1
# \r 表示将光标的位置回退到本行的开头位置,配合end=''不换行，实现打印新内容时删除旧内容
print('\r已处理完第{}个地块'.format(count), end="")
return parcel_mean, parcel_contain_all_poi
In [39]:
integrate_path = u'/home/mw/input/Data5799/数据集/poi_integrate_id.csv' # 2013与2020年集成poi的文件存储路径
original_poi = pd.read_csv(integrate_path, encoding='gb18030') # 读取原始poi数据
k2 = 16 # 地块向量聚类数k的试验区间为[2,26],k为整数
current_year = 2013 # 2013年poi, 作业2中需修改此处2013为2020
# parcel_embed为地块向量，parcel_contain_pois为包含当前年份poi的地块
parcel_embed, parcel_contain_pois = parcel_embedding(poi_embedding, original_poi, current_year)
compute_silhouette(parcel_embed, k2) # 计算每个k值下的平均silhouette值，确定最优k值
这里也是一样，选取平均Silhouette值曲线图中较大Silhouette值对应的k值作为下面聚类的聚类数，即令parcel_k=5。
确定了聚类数以后便可以继续下面的地块向量的kmeans聚类以及聚类结果分布映射可视化了。
In [40]:
parcel_k = 5
parcel_cluster_category = kmeans(parcel_k, parcel_embed) # 使用K-Means聚类
In [41]:
# 将地块向量差的聚类结果、以及两期地块向量的夹角余弦距离写入csv
# change_result为地块变化向量的聚类结果，cosine_dist为两期地块的余弦距离,parcel_contains_all_poi代表两期poi都包含的地块
# ZoneCluster_change_pth为变化向量代表的地块编号及该地块所属的聚类结果的存储路径
def save_result(parcel_path, save_pth, change_result, parcel_contains_all_poi):
parcel = pd.read_csv(parcel_path, encoding='gb18030') # 加载地块数据
'''新建列'''
cluster_num = np.max(change_result) + 1 # 聚类数量，因为从0开始，所以要+1
parcel['Cluster'] = [cluster_num] * parcel.shape[0] # 以最后一类作为不包含poi地块的分类，初始化这一列时数量等于parcel的行数量
count = 0
# 外层循环为地块的顺序,对有poi的地块重新赋值
for order in parcel_contains_all_poi:
parcel.loc[parcel['FID'] == order, 'Cluster'] = change_result[count] # 写入变化向量的聚类结果
count += 1
# \r 表示将光标的位置回退到本行的开头位置,配合end=''不换行，实现打印新内容时删除旧内容
print('\r已处理完第{}个地块'.format(count), end="")
parcel.drop(columns=['城市'], inplace=True) # 删除城市列, 便于arcgis属性连接
# 输出
parcel.to_csv(save_pth, index=False, encoding='gb18030') # 设置index=False是为了不输出unnamed列
In [42]:
parcel_pth = u'/home/mw/input/Data5799/数据集/shenzhen_parcel.csv' # 地块文件路径
save_pth = '/home/mw/project/parcel_cluster_result_{}.csv'.format(current_year)
save_result(parcel_pth, save_pth, parcel_cluster_category, parcel_contain_pois)
print('\n地块向量聚类结果输出完成！')
已处理完第4125个地块
地块向量聚类结果输出完成！
土地利用分类效果展示
利用输出的csv文件，结合ArcGIS或者QGIS软件，将地块向量的聚类结果映射到深圳市地图中。最终的效果图是这样：
其中不同的颜色代表不同的类别，灰色的地块代表的是没有poi存在的地块。
排版：点点GIS
文章授权转载
- END -
ArcGIS制图技巧：制图入门与点、线、面状符号制作
ArcGIS综合制图完整版.doc（文档及练习数据可下载）
如何采用POI进行道路选取？
全国0.8米影像、全球矢量、POI...各种数据高频更新！星图地球数据云正式发布
实景三维模型和地形三维模型有什么区别？