使用 GPT 3.5 和长期记忆的生成式问答
原文：
towardsdatascience.com/generative-question-answering-with-long-term-memory-c280e237b144
探索增强检索的机器学习新世界
·发布于 Towards Data Science ·6 分钟阅读·2023 年 2 月 17 日
–
图片由 Bret Kavanaugh 在 Unsplash 提供。最初发布于 pinecone.io，作者在该网站工作*。
生成式人工智能在 2022 年激发了几次*“哇”*的时刻。从像 OpenAI 的 DALL-E 2、Midjourney 和 Stable Diffusion 这样的生成艺术工具，到像 OpenAI 的 GPT-3.5 这一代模型、BLOOM 以及 LaMDA 和 ChatGPT 这样的聊天机器人。
生成式人工智能（Generative AI）正经历兴趣和创新的蓬勃发展，这并不令人惊讶[1]。然而，这仅仅是生成式人工智能广泛应用的第一年：这是一个新兴领域的早期阶段，准备颠覆我们与机器互动的方式。
最具启发性的用例之一属于Generative Question-Answering（GQA）。通过使用 GQA，我们可以为信息检索（IR）塑造类人互动。
我们每天都在使用信息检索系统。谷歌搜索对网络进行索引，并根据你的搜索词检索相关信息。Netflix 利用你在平台上的行为和历史推荐新的电视节目和电影，而亚马逊在产品方面也有类似的做法[2]。
这些信息检索（IR）的应用具有改变世界的潜力。然而，它们可能仅仅是我们在未来几个月和几年中看到的 IR 与生成式问答（GQA）结合的微弱回响。
想象一下，一个谷歌可以基于前 20 页的智能和深刻的摘要回答你的查询——突出关键点和信息来源。
现有技术已经使这变得可能且出乎意料的简单。本文将探讨增强检索的 GQA 及其如何与 Pinecone 和 OpenAI 一起实现。
最直接的 GQA 系统只需要用户文本查询和一个大型语言模型（LLM）。
最简单的 GQA 系统。
我们可以通过 OpenAI 访问世界上最先进的 LLM 之一。首先，我们需要注册一个 API 密钥。
注册账户后，可以通过点击你的账户（右上角）> 查看 API 密钥 > 创建新秘密密钥来创建 API 密钥。
然后我们切换到 Python 文件或笔记本，安装一些先决条件，并初始化与 OpenAI 的连接。
从这里，我们可以使用 OpenAI 完成端点来询问类似于 “谁是第十二位登月者，他们何时着陆？” 的问题：
我们立即得到一个准确的答案。然而，这个问题相对简单，如果我们问一个鲜为人知的话题会发生什么呢？
尽管这个答案在技术上是正确的，但它并不是一个真正的答案。它告诉我们使用监督学习方法来学习句子之间的关系。这两个事实都是正确的，但并没有回答原始问题。
有两种方法可以让我们的 LLM 更好地理解主题，并更准确地回答问题。
-
我们对包含微调句子转换器领域的文本数据进行 LLM 微调。
-
我们使用检索增强生成，即我们在 GQA 过程中添加一个信息检索组件。添加检索步骤可以让我们检索相关信息，并将其作为二级信息来源输入到 LLM 中。
在接下来的部分中，我们将概述如何实现第二种选择。
构建知识库
实现检索的第二种选择，我们需要一个外部的“知识库”。知识库既是存储信息的地方，又是有效检索这些信息的系统。
知识库是一个信息存储库，可以作为 GQA 模型的外部参考。我们可以将其视为 AI 系统的*“长期记忆”*。
我们将可以检索语义相关信息的知识库称为向量数据库。
向量数据库存储使用特定 ML 模型编码的信息的向量表示。这些模型对语言有一定的“理解”，可以将具有相似含义的段落编码到相似的向量空间，将不相似的段落编码到不同的向量空间。
我们可以通过 OpenAI 的嵌入端点实现这一点：
我们需要对许多记录重复执行嵌入过程，这些记录将作为我们管道的外部信息来源。这些记录仍需下载并准备好以进行嵌入。
数据准备
我们在知识库中使用的数据集是托管在 Hugging Face Datasets上的jamescalam/youtube-transcriptions数据集。它包含来自几个 ML 和技术 YouTube 频道的音频转录。我们用以下方式下载它：
数据集包含许多小的文本片段。我们需要合并几个片段，以创建包含更多有意义信息的更大文本块。
创建文本块后，我们可以开始初始化我们的知识库，并用数据填充它。
创建向量数据库
向量数据库是我们管道中的存储和检索组件。我们使用 Pinecone 作为我们的向量数据库。为此，我们需要注册一个免费的 API 密钥，并在下面输入它，以创建用于存储数据的索引。
然后我们这样嵌入和索引数据集：
我们准备将 OpenAI 的Completion和Embedding端点与我们的 Pinecone 向量数据库结合，创建一个检索增强 GQA 系统。
OP Stack
OpenAI Pinecone（OP）技术栈是构建高性能 AI 应用（包括检索增强 GQA）的越来越受欢迎的选择。
我们在查询时间的管道包括以下内容：
-
OpenAI
Embedding端点用于创建每个查询的向量表示。 -
Pinecone 向量数据库用于从以前索引的上下文数据库中搜索相关的段落。
-
OpenAI
Completion端点用于生成自然语言答案，考虑到检索到的上下文。
我们首先使用相同的编码器模型对查询进行编码，以创建查询向量xq。
查询向量xq用于通过index.query查询 Pinecone，并与之前索引的段落向量进行比较，以找到最相似的匹配项——在上面的res中返回。
使用这些返回的上下文，我们可以构建一个提示，指示生成 LLM 根据检索到的上下文回答问题。为了简单起见，我们将在一个名为retrieve的函数中完成所有这些工作。
请注意，为了可读性，生成的扩展查询（query_with_contexts）已被缩短。
从retrieve中，我们生成一个更长的提示（query_with_contexts），其中包含一些指令、上下文和原始问题。
然后，提示通过 OpenAI 的Completion端点输入生成 LLM。像以前一样，我们使用complete函数处理所有内容。
由于额外的*“源知识”*（直接输入模型的信息），我们消除了 LLM 的幻觉——生成了准确的回答。
除了提供更多的事实性回答外，我们还有 Pinecone 提供的信息来源用于生成我们的回答。将其添加到下游工具或应用中可以帮助提高用户对系统的信任，让用户确认所提供信息的可靠性。
这就是关于检索增强型生成式问答（GQA）系统的演示。
如前所示，单独使用大型语言模型（LLMs）效果非常好，但在面对更小众或具体的问题时却常常表现不佳。这通常会导致幻觉，这些幻觉不易被发现，且可能会被系统用户忽视。
通过为我们的 GQA 系统添加*“长期记忆”*组件，我们可以利用外部知识库来提高系统的准确性和用户对生成内容的信任。
自然，这类技术具有广泛的潜力。尽管它是一项新技术，但我们已经在YouChat、一些播客搜索应用中看到它的应用，并且有传闻称它将成为对抗谷歌的挑战者[3]。
任何需要信息的地方都有可能出现颠覆，而检索增强型 GQA 代表了利用当前过时信息检索系统的最佳机会之一。
[1] E. Griffith, C. Metz, 《人工智能的新领域蓬勃发展，即使在技术低迷时》（2023），纽约时报
[2] G. Linden, B. Smith, J. York, 《亚马逊推荐：物品间协作过滤》（2003），IEEE
[3] T. Warren, 《微软通过将 ChatGPT 与 Bing 搜索整合来挑战谷歌》（2023），The Verge
最初发表于 https://www.pinecone.io 。
所有图片均为作者提供，除非另有说明
遗传算法：寻找最佳电子邮件投递时间表以最大化互动
使用进化算法优化消费者银行的 D2C 活动
·发表于 Towards Data Science ·阅读时间 9 分钟·2023 年 11 月 1 日
–
作者使用 Bing Image Creator 生成的图像
某些电子邮件投递时间是否会导致更高的互动？
电子邮件营销人员常面临的一个最常见问题是何时发送电子邮件以最大化打开率、点击率和转化率。这个问题没有明确的答案，因为不同的受众可能有不同的偏好和行为。他们处于哪个时区？他们使用什么设备查看电子邮件？他们的日常生活和时间表是什么？他们多频繁检查电子邮件？这些因素可能会影响他们何时最有可能打开和互动你的电子邮件。
你可以使用像 A/B 测试或分割测试这样的工具来比较在不同时间发送的电子邮件活动的表现。你还可以使用像 Google Analytics 或 Mailchimp 这样的分析工具来跟踪电子邮件活动的指标，例如打开率、点击率、跳出率和转化率。通过分析数据，你可以确定适合你受众和目标的最佳投递时间。
当你对客户在不同时间的点击率和打开率有了充分了解后，下一步是创建一个最佳的投递时间表，以最大化这些指标，同时避免负面影响，例如用户退订——我们在电子邮件营销中称之为“疲劳”现象。
在这篇文章中，我将尝试使用遗传算法来解决这个优化问题——这种方法在营销领域并不常见。
我理解遗传算法可能会让人感到困惑。我们将看到如何使用简单的 Pandas 操作来实现该算法的核心概念。
定义问题陈述
我在下一步行动的文章中介绍了虚构但受欢迎的消费者银行 ULFC。在那次努力中，我们创建了一个强化模型，根据客户过去的响应，建议向其按揭客户推出下一个最佳优惠。现在，ULFC 希望数据科学团队推荐一个最优的时间表——日期、时间、数量——以作为季度电子邮件活动的一部分。目标是最大化客户对电子邮件的参与，而不流失他们。
我将使用这份合成数据，这是我创建的。它包含了针对 1000 名客户的历史点击率、打开率和疲劳率数据，这些客户将作为本次活动的目标。
让我们先来看一下数据。
-
打开率通常较高，可达 30%。
-
点击率更具意图性，范围在 0-3%之间
-
疲劳率（客户取消电子邮件营销的概率）更具上下文性，并且可能因细分/活动而异。我们将假设在 0-5%之间。
import numpy
import pandas
open_rate_df = pd.read_csv("open_rate.csv").set_index("Customer_ID")
click_rate_df = pd.read_csv("click_rate.csv").set_index("Customer_ID")
fatigue_rate_df = pd.read_csv("fatigue_rate.csv").set_index("Customer_ID")
作者提供的图片 — 1000 名客户的每日打开率
作者提供的图片 — 1000 名客户的每日点击率
作者提供的图片 — 1000 名客户的每日疲劳率
遗传算法简介
遗传算法（GA）是一种通过模拟自然进化过程来寻找问题解决方案的方法。遗传算法通过创建一个可能解决方案的种群（称为个体），并对其应用选择、交叉和变异等生物学操作符来工作。这些个体通过一个适应度函数进行评估，该函数衡量它们解决问题的效果。最适应的个体更可能生存和繁殖，而较不适应的个体则被淘汰。通过这种方式，遗传算法可以在搜索空间中进行探索，并随着时间的推移趋向于最优或接近最优的解决方案。以下是任何遗传算法用例中涉及的核心步骤。
作者提供的图片 — 任何遗传算法中的核心步骤
设置参数
客户的数量 K 是 1000。在这次电子邮件营销活动中，我们希望将打开率和点击率最大化。它们的权重由 w₀ 和 wc 分别设定。
假设我们不希望在时间范围内向同一客户发送超过 2 封邮件（由 E 设置）。另一方面，我们希望将每个客户的预期疲劳率限制在 10%以下（由 T 设置）。
K = 1000
w0 = 1
wc = 1
T = 0.1
E = 2
创建初始种群
在这里，我们基本上创建了初始候选者，这些候选者作为我们算法通过选择、交叉和变异来改进的起始种群。每个候选者是一组 0 和 1 的二进制决策变量，表示是否在该时间段发送邮件。
作者提供的图片 — 示例候选者
def initialize_population(num_candidates):
all_candidates = []
for i in range(num_candidates):
# Create random decision variables such that only 1 email is delivered to a person per day
base = np.random.choice([-1,0,1,2], size=(K,5), p=[0.95, 0.05/3, 0.05/3, 0.05/3]).tolist()
data = []
for k in range(K):
vals = []
for t in base[k]:
if t == -1:
vals.extend([0,0,0])
elif t == 0:
vals.extend([1,0,0])
elif t == 1:
vals.extend([0,1,0])
elif t == 2:
vals.extend([0,0,1])
data.append(vals)
all_candidates.append(
pd.DataFrame(
data = np.array(data),
index = customer_ids,
columns = timeslots
)
)
return all_candidates
定义适应度函数
目标是根据每个接收者的偏好和行为选择最佳的邮件发送日期。我们希望增加他们打开和点击邮件的可能性。我们假设每个接收者在每一天和时间段的打开率和点击率不同。我们为打开率和点击率添加权重，这可以根据是转换还是意识提升活动进行调整 [1]。因此，我们需要找到每个接收者得分最高的日期和时间段。
我们需要考虑两个限制：
-
疲劳限制： 每个接收者在一定的频率下会感到疲倦，这种频率因人而异。如果我们发送过多的邮件，他们可能会退订。我们需要将每个接收者的总疲劳率保持在一定的限制内。
-
总投递限制： 这控制了我们可以发送给一个接收者的最大邮件数量。
作者提供的图片 — 目标函数
在下面的代码中，我已经将这一目标转换为一个适应度函数，该函数返回目标分数，并根据约束 1 和 2 是否满足返回布尔值。
def fitness_function(candidate_df):
obj = ((w0 * open_rate_df + wc * click_rate_df) * candidate_df).values.sum()
const1 = ((fatigue_rate_df * candidate_df).sum(axis=1) > T).sum() == 0
const2 = (candidate_df.sum(axis=1) > E).sum() == 0
return obj, const1, const2
def evaluate_fitness(x):
score, const1, const2 = fitness_function(x)
if const1 and const2: return score
else: return -1
现在让我们创建一个简单的函数来测量一些关键指标，以帮助我们了解投递时间表的质量。
def print_candidate_metrics(candidate):
open_rates = (open_rate_df * candidate).values.flatten()
open_rate = round(open_rates[open_rates > 0].mean() * 100, 1)
expected_opens = round(open_rates[open_rates > 0].sum(), 1)
print(f"Open Rate = {open_rate}%")
print(f"Expected number of opens = {expected_opens}")
click_rates = (click_rate_df * candidate).values.flatten()
click_rate = round(click_rates[click_rates > 0].mean() * 100, 1)
expected_clicks = round(click_rates[click_rates > 0].sum(), 1)
print(f"Click Rate = {click_rate}%")
print(f"Expected number of clicks = {expected_clicks}")
fatigue_rates = (fatigue_rate_df * candidate).values.flatten()
fatigue_rate = round(fatigue_rates[fatigue_rates > 0].mean() * 100, 1)
print(f"Fatigue Rate = {fatigue_rate}%")
首先，我们将检查一些初始随机选择的候选者的指标作为基准。
initial_candidates = initialize_population(3)
for i in initial_candidates:
print_candidate_metrics(i)
print()
作者提供的图片 — 初始随机生成的候选者时间表的指标
我看到初始候选者在时间范围内平均生成 35–40 次打开。我们还看到大约 4 次预期点击。
让我们跟随这个过程，构建算法，看看它是否能改善这些指标。
自然选择中的过程
驱动代际过程的主要过程是 选择。在遗传算法的实践中，这有许多定义。最简单的一种是根据适应度分数对所有候选者进行排名，并选择前 N 名将其属性传递给下一代。这种方法被称为精英主义或选择精英候选者的过程。
另外两个过程是交叉和变异。这两个操作是遗传算法中的基本操作，有助于从现有的解决方案中创造新的解决方案。交叉将两个或多个父代候选者的特征结合起来，生成子代候选者，而变异则在候选者中引入随机变化，以创造多样性并避免局部最优解。交叉和变异的概率控制着种群中候选者的变化率。较高的交叉概率意味着将通过结合父代生成更多的子代，而较低的变异概率则意味着较少的子代会随机改变。这些概率可以是固定的，也可以是自适应的，具体取决于问题和适应度函数。
我打算实施一个简单的交叉操作，我将随机选择一个行号（即交叉点），并从第一个子代中提取该行号之前的所有行，从第二个子代中提取该行号之后的所有行。
def crossover(candidate1, candidate2):
crossover_point = random.randint(1,K-1)
child1 = pd.concat([candidate1.iloc[:crossover_point,:], candidate2.iloc[crossover_point:,:]])
child2 = pd.concat([candidate1.iloc[:crossover_point], candidate2.iloc[crossover_point:]])
return child1, child2
变异有点有趣。我的实现是单点变异。我随机选择一个时间段中的一天（即变异点），并改变那一天的所有决策（0 和 1）。
def mutate(candidate):
mutation_point = random.randint(0,14)
new_candidate = candidate.copy()
new_candidate.iloc[:,mutation_point] = np.random.choice([0, 1], size=K, p=[.99, .01])
return new_candidate
通常，变异对于探索搜索空间的新区域非常有用，而交叉则有助于利用局部最优解。要更好地理解探索和利用，查看一下我的文章关于下一步行动模型。利用这些功能，让我们构建 ULFC 的算法，将所有部分拼接在一起。
最终算法
首先，我们有几个参数控制代数（把它看作是迭代）、每代结束时允许的种群大小、交叉和变异的概率，以及精英数量（用于选择过程）。
generations = 100
population_size = 50
crossover_rate = 0.8
mutation_rate = 0.2
elitism_size = 2
如之前的流程图所示，我们的算法有 5 个步骤。
-
创建随机初始种群。
-
使用适应度函数评估候选者。
-
选择精英。
-
应用交叉以创建新的子代候选者。
-
对新创建的候选者应用变异。
重新评估并重复多代。
def optimal_schedule():
# Create initial candidates
population = initialize_population(population_size)
for generation in range(generations):
new_population = []
# Elitism
population.sort(key=lambda x: evaluate_fitness(x), reverse=True)
new_population.extend(population[:elitism_size])
while len(new_population) < population_size:
parent1, parent2 = random.choices(population, k=2)
if random.random() < crossover_rate:
child1, child2 = crossover(parent1, parent2)
else:
child1, child2 = parent1.copy(), parent2.copy()
if random.random() < mutation_rate:
child1 = mutate(child1)
if random.random() < mutation_rate:
child2 = mutate(child2)
new_population.extend([child1, child2])
population = new_population
return population[0]
最后，让我们运行刚刚构建的算法。
result = optimal_schedule()
print_candidate_metrics(result)
作者提供的图片 — 最优时间表的指标
最优时间表的预期打开次数和点击次数比初始候选者更高。你可以在更多代中运行此算法以进一步优化。你甚至可以进行时间段分析，以跟踪推荐的时间表变化。
作者提供的图片
总结
媒介如电子邮件和短信可以是非常具有成本效益和说服力的渠道来接触受众。但进一步优化电子邮件营销的时间安排可以确保企业从其电子邮件中获得最大互动。这些项目还解锁了许多后续分析。例如，流行的 CRM 软件 HubSpot 的团队每年发布关于“发送电子邮件的最佳时间”的报告 [2]。根据他们 2023 年的发现，电子邮件从晚上 9 点到 12 点获得的互动最多。然而，这些数据来自对他们客户的调查。如果你在公司内部进行时间安排优化，你可以以更高的准确性和细致度产生这样的见解，从而适应你的业务案例。
最后，作为一种警告，我建议谨慎使用遗传算法。只有当你知道如何在你的问题背景下解释和说明它时。它只是另一种优化算法，还有许多可能更容易解释和更快的算法。另一方面，遗传算法是稳健的，并且可以很好地达到全局最优。
感谢阅读 😁
参考文献：
-
Zhang, L., He, J., Yan, Z., Dai, W., Pani, A. (2020). 遗传算法在多目标电子邮件营销投递问题中的应用
-
发送电子邮件的最佳时间 [2023 研究] 由 Flori Needle 撰写 —
blog.hubspot.com/marketing/best-time-to-send-email
天才圈子：绘制诺贝尔网络
原文：
towardsdatascience.com/genius-cliques-mapping-out-the-nobel-network-e9350552084
诺贝尔奖获得者的网络。图像由作者提供。
·发表在 Towards Data Science ·7 min read·2023 年 9 月 27 日
–
结合网络科学、数据可视化和维基百科，揭示所有诺贝尔奖获得者之间隐藏的联系。
这篇文章最初发表在《Nightingale》第二期，这是数据可视化协会的印刷杂志。
启发性的科学家：诺贝尔圈内外
尽管我在网络与数据科学方面获得了博士学位，但每当寻找灵感时，我总是紧密关注我在物理学方面的根源。成长在匈牙利，我特别惊叹于“火星人”这一群体的成就，这些著名科学家在二战期间从匈牙利移居到美国。有趣的是，他们中的一些人甚至上过同一所高中。
火星人包括例如列奥·西拉德，他不仅发现了核链反应理论，还与阿尔伯特·爱因斯坦和尤金·维格纳——曼哈顿计划的关键科学家——共同申请了冰箱专利，推动了第一台核反应堆的发展。由于他的贡献，维格纳在 1963 年获得了诺贝尔物理学奖，这 18 个与匈牙利背景的思想家相关的诺贝尔奖之一。
那 18 个奖项仅占所有诺贝尔奖的约百分之三。事实上，自 1901 年以来，大约 600 个奖项授予了不到一千位获奖者，涉及领域包括物理学、化学、生理学或医学、文学、和平，以及从 1969 年起的经济学。网站 NobelPrize.org 突出了有关奖项及其获奖者的其他有趣统计数据：从最年轻（17 岁）和最年长（97 岁）的获奖者，到像约翰·巴尔丁（物理学，1956 年和 1972 年）、林纳斯·鲍林（化学 1954 年，和平 1962 年）以及玛丽·居里（物理学 1903 年，化学 1911 年）等多次获奖者。
居里家族主导了诺贝尔奖。玛丽·居里首次与丈夫皮埃尔分享了一个奖项，后来又获得了第二个奖项。此外，这对伟大的夫妇还培养了一位诺贝尔奖得主。他们的女儿伊雷娜·居里与丈夫弗雷德里克·朱利奥共同获得了 1935 年化学奖。玛丽·居里还曾是另一个精彩的诺贝尔得主小圈子的成员（遗憾的是，皮埃尔于 1906 年去世）：1911 年的索尔维会议。这可能是科学史上最令人印象深刻的阵容之一：29 名与会者中有 27 人已经获得或后来获得了诺贝尔奖。
建立诺贝尔网络
火星人、居里家族、曼哈顿计划和索尔维会议的故事都表明，在幕后，一些紧密交织的社会网络在诺贝尔奖得主之间运作。为了追踪这个网络，我访问了最广泛使用的在线百科全书维基百科，并收集了每位获奖者的维基页面文本。
然后，在每位获奖者的页面文本中，我统计了所有其他名字的提及情况，并注意是否有任何配对共享了在维基百科上值得注意的共同历史。通过这种方式，我建立了一个包含 682 个节点和 588 条链接的网络，其中节点对应于获奖者，两个节点之间的链接强度与他们的维基页面互相引用的次数成正比。此外，我下载了每位获奖者页面的总浏览量，并将其网络节点大小设置为该数字的对数。这个节点缩放最终突出了那些已经成为家喻户晓的名字。为了完成网络可视化，我应用了与科学学科相对应的颜色编码。您可以在图 1 中找到结果。
图 1. 诺贝尔网络。根据维基百科页面之间的交叉引用，展示了至少有一个连接的诺贝尔奖得主网络。每个节点对应一个获奖者，边的宽度表示交叉引用的数量，节点的大小与其维基页面的总浏览量成正比。颜色编码表示他们所获奖项的学科（在有多个不同奖项的情况下，从获奖学科中随机选择一种颜色）。视图数量最高的节点被标记出来。
诺贝尔网络的启示
对我来说，作为一个网络科学家，关于网络的第一个也是最引人注目的观察是它的核心-边缘分离：中心的一个大连接组件（所谓的巨型组件），包含超过 30%的节点，以及其周围的一个碎片化环圈，包含较小的网络组件，大小最多为十个节点。最常见的组件大小为两个和三个节点，这与诺贝尔奖可以在最多三位获奖者之间分享的事实非常一致，而共享奖项在大多数领域变得越来越普遍。
我还发现，巨型组件中的节点较大，这意味着这些获奖者的可见性显著更高，并且搜索命中数量也更多，这通过他们的 Wiki 浏览量的对数来衡量。通过查看数据，结果显示中央组件的 Wiki 浏览量中位数为 351,005，而外环组件仅为 170,510，并且中央组件的平均浏览量大约是外环组件的 2.5 倍。因此，看起来中央 clique 的人气远远更高！
他们是谁？绿色至黄色的色调与红色调的对比旨在区分科学与人文学科，这与巨型组件的左右两侧相对应。这些侧面由詹姆斯·查德威克爵士连接，他因发现中子而获得了 1935 年的诺贝尔物理学奖，并且还成为了联合国的科学顾问。科学一侧（图 2a），以阿尔伯特·爱因斯坦和马克斯·普朗克等研究人员为代表，似乎在普鲁士科学院（1700–1945）中有着深厚的根基，并且在现代物理学的奠基人中也很强大，从居里夫人到恩里科·费米和尤金·维格纳或乔尔吉·赫维希（都有匈牙利和火星的根源）。
图 2a. 图 1 的放大视图，重点关注科学领域的 clique。
在人文学科一侧（图 2b），我们可以看到一些相当受欢迎的人物。显然，科学并不是通向世界声誉的途径！有两个突出的中央获奖组织立即引人注目：欧盟和联合国，它们都获得了诺贝尔和平奖。值得注意的个人包括杰出的政治家，如巴拉克·奥巴马或亨利·基辛格，人权活动家纳尔逊·曼德拉，以及经济学家米尔顿·弗里德曼（有匈牙利，但非火星的根源）。
至于外环部分，有一些著名的社交人物，如欧内斯特·海明威、温斯顿·丘吉尔、富兰克林·D·罗斯福和理查德·费曼——个人来说，我最喜欢的诺贝尔奖获得者，因为他既有科学贡献，又有顽皮和古怪的个性。这些人尽管生活繁忙，但在网络中相对孤立，这可能是由于他们活跃时期的时间和地理位置与其他获奖者的差异。此外，这里的数据可能不完整，因为维基百科并不完美，也不百分之百准确地记录社交关系，而且遗憾的是（或者说并非如此？）当时并没有 Facebook。
最后，匈牙利人和火星人：从数据来看，许多人甚至与诺贝尔奖得主没有任何联系，而那些属于网络中的成员也只是分散在各处。这背后的原因不清楚——也许火星人的传说被夸大了，或者可能他们中没有足够多人获得诺贝尔奖以至于在网络中显现。虽然可以确定的是，曼哈顿计划在运作期间有七位诺贝尔奖得主，后来又有十几位，但其中只有维根来自火星人。
图 2b. 图 1 的放大，重点关注人文学科中的小圈子。
结论
尽管在《诺贝尔网络》中扫描这些名字和联系非常激励人心，并且让我真切地感觉自己是在“站在巨人的肩膀上”，但网络还是存在缺陷。除了边缘的东欧人，我们看到一个精英俱乐部在中心出现，主要的知名名字集中在巨大的组件中，排除了三分之二的网络。这表明，三分之二的获奖者仅仅是领奖后回到自己的工作中，只有剩下的三分之一参与了可见的联系，无论是友谊还是合作。正如“整体大于部分之和”，错过了超过 60%的那些才智横溢的头脑从中心的思想流动中，似乎是一种遗憾。
更多的错失机会出现了。中心组件本身被分为两个阵营：科学和人文学科。这种极化与今天的主要方向——跨学科研究非常不符，这使我们能够解决前所未有的重大社会问题。此外，网络揭示了女性获奖者的数量较少。尽管玛丽·居里的历史卓越，但女性获奖者仅约占 6%，其中大多数获得了和平奖（109 位中的 16.5%），而在物理学奖中获奖者最少（219 位中的 1.8%）。
但一切并未丧失。像这样的映射练习可以帮助揭示这些问题，否则，即使是最狂热的诺贝尔奖迷也几乎看不到。放大并利用网络科学可以突出那些隐藏的模式，并启发理解，这是识别未来解决方案的第一步，无论是关于性别差距还是精英小圈子。
群体上的几何深度学习
Continuous vs discrete approaches on the sphere
·
关注 发表在 Towards Data Science ·6 min read·Mar 6, 2023
–
理想情况下，群体上的几何深度学习技术应该能够编码对群体变换的等变性，以提供良好的表示空间和出色的性能，同时也要具备计算效率。然而，没有一种方法能够同时提供这两种理想的属性。连续的方法提供了优秀的等变性，但计算成本非常高。离散的方法通常计算效率相对较高，但牺牲了等变性。我们指出了未来能够兼具这两者优点的技术。
照片由 Serg Antonov 提供，发布在 Unsplash
群体上的深度学习是几何深度学习一个快速增长的领域（参见我们最近的 TDS 文章 几何深度学习简要介绍）。群体包括具有全局对称性的同质空间，其中最典型的例子是球面。
几何深度学习在群体上的实际应用很普遍，尤其是在球面上。例如，球形数据在许多应用中出现，不仅在数据直接获取自球面时（例如地球上的数据或通过 360° 相机捕捉全景照片和视频时），还包括考虑球面对称性时（例如分子化学或磁共振成像）。
我们需要在群体上既高效又可扩展到大规模高分辨率数据集的深度学习技术。一般来说，这个问题仍然没有解决。
一个球形数据的例子。[照片由 NASA 提供，来自 Unsplash]
目标
深度学习技术之所以如此有效，部分原因是现代架构中编码的归纳偏差。
一个特别强大的归纳偏差是编码数据已知满足的对称性（如我们 TDS 文章 爱因斯坦可以教我们什么关于机器学习*）所阐述的）。例如，卷积神经网络（CNNs）编码了平移对称性，或更确切地说，平移等变性，如下图所示。
平移等变性的说明。给定一张图像（左上），应用卷积核（𝒜）以获得特征图（右上），然后平移（𝒯）特征图（右下），等同于首先平移图像（左下），然后应用卷积核（右下）。[原始图形由作者创建。]
在深度学习架构中编码等变性会产生行为良好的特征空间，使得学习可以非常有效地进行。
因此，对于群体上的几何深度学习，我们希望编码对各种群体变换的等变性，这通常会带来非常好的性能。然而，在一般的群体设置中，这变得非常计算密集——在许多情况下是不可行的。
如何在深度学习架构中以计算上可扩展的方式编码等变性是一个活跃的研究领域。
群体卷积
卷积的概念，负责 CNN 架构在平面图像上的巨大成功，自然地编码了等变性，并且可以推广到群体设置。
群上信号（即数据、特征图）f 与滤波器 𝝭 的群卷积表示为
其中 g 是群 G 的一个元素，而 dµ(u) 是积分的（Haar）测度。上述表达式与更常见的平面设置中的卷积完全类似。我们对滤波器进行变换（对于平面 CNN 来说是平移），与感兴趣的信号相乘，然后求和，即积分。
在球面上，我们考虑由 3D 旋转给出的变换，因此球面上的信号卷积表示为
其中 R 表示旋转，ω 为球面坐标。
一旦定义了群上的卷积，我们就可以以类似于标准平面 CNN 的方式在群上构建 CNN。即，通过组合卷积和逐点非线性激活（还包括适当地构建在群上的池化和归一化层）。
问题仍然是：我们如何在实践中计算群卷积？
一方面，我们希望实现能够准确捕捉卷积上的等变性特性。另一方面，我们希望实现具有高计算效率。如我们将看到的，现有方法通常只能捕捉这些要求中的一个，而不能同时满足两个要求。
离散球面 CNN 方法
现有的方法可以大致分为离散和连续方法。
离散方法处理数据的离散版本，通常是像素或图形表示，这通常具有很高的计算效率。然而，通常不存在规则的离散化。
以球面为例，众所周知，球面的规则离散化并不存在。因此，没有方法可以以不变于旋转的方式对球面进行离散化，如下图所示。
在球面上旋转一组像素会得到一组不能与现有集合重叠的像素。这在球面的所有采样中都是正确的。[原始图形由作者创建。]
在离散空间上直接定义的操作无法严格捕捉等变性。
因此，离散方法在计算性能上具有优势，但代价是等变性。
连续球面 CNN 方法
作为上述离散方法的替代方案，也可以考虑信号的连续表示。
球面上的函数可以通过球面谐波的展开表示（如下图所示）。对于带限信号，可以在有限的样本集合中捕捉到信号的所有信息内容，从中可以准确地计算球面谐波系数[1]。这类似于著名的奈奎斯特-香农采样定理扩展到球面。
球面谐波函数。[图像来源于 Wikimedia Commons。]
由于球面是一个紧致流形，其谐波空间是离散的。因此，通过使用有限的球面谐波空间表示，可以访问底层的连续信号。
各种球面 CNN 架构已经被构建，其中卷积通过其谐波表示来计算[2–6]。通过访问底层的连续信号，这些方法实现了卓越的等变性。然而，它们涉及反复执行球面谐波变换，这在计算上是昂贵的。
连续方法准确捕捉旋转等变性，但计算要求高。
二分法：离散 vs 连续方法
如上所示，离散和连续方法之间存在一种二分法，如下图所示。理想情况下，我们希望技术既具有等变性又具有计算可扩展性。
然而，连续方法提供等变性，但计算成本很高。另一方面，离散方法通常计算效率较高，但牺牲了等变性。
连续和离散几何深度学习技术在群体上的二分法。[原始图由作者创建。]
破除二分法
我们期望在群体上进行几何深度学习的技术能够提供等变性（这通常转化为良好的表示空间和卓越的性能），同时也具有计算可扩展性。
在我们的下一篇文章中，我们将描述一种新的混合离散-连续（DISCO）方法，最近被 ICLR 接受，这种方法正好实现了这些目标[7]。通过保持表示的一些组件为连续的，我们实现了卓越的等变性，而其他组件则被离散化以提供高效的可扩展计算。
参考文献
[1] McEwen & Wiaux, A novel sampling theorem on the sphere, IEEE TSP (2012), arXiv:1110.6298
[2] Cohen, Geiger, Koehler, Welling, Spherical CNNs, ICLR (2018), arxiv:1801.10130。
[3] Esteves, Allen-Blanchette, Makadia, Daniilidis, 学习 SO(3)等变表示的球面卷积神经网络，ECCV（2018），arXiv:1711.06721。
[4] Kondor, Lin, Trivedi, Clebsch-Gordan 网络：一种完全傅里叶空间的球面卷积神经网络，NeurIPS（2018），arXiv:1806.09231
[5] Cobb, Wallis, Mavor-Parker, Marignier, Price, d’Avezac, McEwen, 高效的广义球面卷积神经网络，ICLR（2021），arXiv:2010.11661
[6] McEwen, Wallis, Mavor-Parker*，球面上的散射网络：用于可扩展和旋转等变的球面卷积神经网络，ICLR（2022），arXiv:2102.02828
[7] Ocampo, Price, McEwen, 通过离散-连续（DISCO）卷积实现可扩展和等变的球面卷积神经网络，ICLR（2023），arXiv:2209.13603
机器学习与经典统计中的线性回归几何解释
通过视觉和分析来揭示线性回归的困惑
·
关注 发表在 Towards Data Science ·9 分钟阅读·2023 年 12 月 12 日
–
图片：线性回归插图，作者 Stpasha，通过维基共享资源（公共领域）。原始图片链接：upload.wikimedia.org/wikipedia/commons/8/87/OLS_geometric_interpretation.svg
上述图像表示了普通最小二乘法（OLS）或线性回归的几何解释（在经典统计中这两个词可以互换使用）。让我们以直观的方式解析我们所看到的内容：
-
变量（X1 和 X2）： 想象你有两个变量 X1 和 X2。这些变量可以代表任何东西，比如你学习的小时数和你参加的练习考试次数。
-
数据点（y）： 现在，你有了你的结果或你试图预测的内容，我们称之为‘y’。在我们的例子中，这可能是你在实际考试中的成绩。
-
平面（colX）： 平面代表通过组合不同数量的变量 X1 和 X2 可以得到的所有可能预测值。在我们的例子中，它可能代表基于不同的学习时间和练习考试次数你可能预测的所有可能考试成绩。
-
估计系数（Beta1 和 Beta2）： 这些是 OLS 方法对每个变量对你的成绩影响程度的最佳猜测。因此，Beta1 可能告诉你每多学一个小时成绩预计会增加多少，而 Beta2 可能告诉你每增加一个练习考试成绩会增加多少。
-
预测点（XB^）： 这是基于估计系数你得到的预测成绩。它位于平面上，因为它是使用 OLS 估计的变量 X1 和 X2 的组合。
-
实际点（y）： 这是你的实际考试成绩。
-
误差（ε）： 这是你实际成绩和预测成绩之间的差异。换句话说，就是预测与现实的偏差。
那么，OLS 如何运作呢？
OLS 尝试找到 Beta1 和 Beta2 的值，使得当你使用 X1 和 X2（学习小时数和练习考试次数）来预测‘y’（考试成绩）时，所有数据点的误差（ε）尽可能小。在图形中，这就像是调整平面，直到垂直虚线（表示误差）总长度尽可能短。实际数据点（y）到平面（colX）的最短距离总是垂直于平面的直线。OLS 找到一个特定的平面，使这些垂直距离对于所有点都达到最小。
换句话说，OLS 试图“拟合”平面，使其尽可能接近你的实际成绩，同时认识到它通常不会通过所有实际点，因为现实生活中很少有那么完美的情况。
这就像把最好的纸张贴在一群铅笔点下面，使纸张尽可能接近所有点。
让我们回顾 OLS 的主要假设，并将其与上述视觉图像联系起来：
1. 线性性
假设： 自变量（X1，X2）与因变量（y）之间的关系是线性的。
视觉解释： 在图像中，这就是我们使用平面（colX）来表示 X1 和 X2 组合的原因。如果关系不是线性的，我们无法用平面来表示，它会是曲线或其他形状。
2. 独立性
假设： 观察值彼此独立。
视觉解释： 每个数据点（表示一个观察值）是独立绘制的。如果存在依赖性，我们会在误差（ε）中看到系统的模式，比如它们都集中在平面的一侧，这表明一个数据点的位置可能预测另一个，违反了这一假设。
3. 同方差性
假设： 误差项（ε）的方差在所有自变量水平上保持不变。
视觉解释： 理想情况下，实际数据点（y）到预测平面（colX）的垂直距离应该均匀分布。这些距离不应有漏斗形状或模式；它们应显得随机。如果误差随着 X1 或 X2 的增加而增大或减小，这将违反同方差性。
4. 无完美多重共线性
假设： 自变量之间没有完全相关性。
视觉解释： 在图示中，X1 和 X2 由两个方向不同的箭头表示。如果它们完全相关，它们会指向完全相同的方向，我们将不会有平面而是直线。这将使我们无法估计 X1 和 X2 对 y 的独特影响。
5. 无自相关性
假设： 误差项之间不相关。
视觉解释： 这个假设涉及误差项，虽然在图像中没有明确显示，但我们推断每个误差项（ε）是随机的，不受前一个或下一个误差项的影响。如果存在模式（例如，一个误差总是跟随另一个相似大小的误差），我们会怀疑自相关性。
6. 外生性
假设： 误差项的期望值为零。
视觉解释： 这意味着平面应当被定位，使得误差在平均值上互相抵消。一些数据点会在平面之上，一些会在平面之下，但没有系统性偏差使它们都在平面之上或之下。
7. 误差的正态性（通常是假设检验的假设）
假设： 误差项呈正态分布。
视觉解释： 虽然正态性假设在数据和模型的三维图中无法直观呈现，但如果我们查看误差项的直方图，我们期望看到正态分布的熟悉钟形曲线。
线性回归在机器学习宇宙中如何不同于经典统计学中的普通最小二乘法线性回归？
在经典统计学中，普通最小二乘法（OLS）可以通过最大似然估计（MLE）的视角来接近。MLE 和 OLS 都旨在找到模型的最佳参数，但它们来自不同的理论，并使用不同的方法来实现这一目标。
最大似然估计 (MLE) 方法：
MLE 基于概率。它提出的问题是：“给定一组数据点，什么是最可能的模型参数，这些参数能够生成这些数据？”MLE 假设误差具有某种概率分布（通常是正态分布），然后找到最大化观察到实际数据的可能性的参数值。在几何视觉中，这类似于调整平面（colX）的角度和位置，使得看到实际数据点（y）的概率最高。似然不仅考虑点到平面（误差）的距离，还考虑误差分布的形状。
机器学习 (ML) 中目标函数的最小化：
另一方面，ML 方法通常将回归框架作为优化问题。目标是找到那些使某个目标函数（通常是平方误差总和（SSE））最小的参数。这是一种比 MLE 更直接的方法，因为它对误差的底层概率分布假设较少。它只是试图在平方意义上将数据点到预测平面的距离尽可能小，以更严厉地惩罚较大的误差。几何解释是你在倾斜和移动平面（colX），以最小化从实际点（y）到平面的垂直距离的平方和（虚线）。
两者比较：
尽管程序有所不同——一种是基于概率的方法，另一种是优化技术——但在 OLS 的情况下，它们往往会产生相同的结果。这是因为当误差服从正态分布时，线性模型系数的 MLE 会得出与最小化平方误差总和相同的方程。从视觉上看，两种方法实际上都是在尝试在变量 X1 和 X2 的空间中定位同一个平面，以便它最好地表示与 y 的关系。
主要区别在于解释和潜在的泛化。MLE 的框架允许在建模误差结构时有更多灵活性，并且可以扩展到误差不被假设为正态分布的模型中。ML 方法通常更直接，专注于平方差的总和，而不关注底层概率分布。
总之，虽然最大似然估计 (MLE) 和最小化方法可能会得到相同的最小二乘回归系数，但它们在概念上是不同的。MLE 是概率性的，根植于在给定模型下观察数据的可能性，而最小化方法是算法性的，专注于直接减少误差。几何视觉对两者来说是相同的，但平面位置的原理不同。
附加：当你在上述解释中引入正则化时会发生什么？
正则化是一种用于防止模型过拟合的技术，过拟合可能发生在模型过于复杂，开始捕捉数据中的噪声而不仅仅是数据的真实底层模式。正则化有不同的类型，但最常见的两种是：
-
套索回归 (L1 正则化)： 这增加了一个等于系数绝对值的惩罚项。它可以将一些系数减少为零，有效地进行特征选择。
-
岭回归 (L2 正则化)： 这增加了一个等于系数平方的惩罚项。所有系数都被相同的因子缩小，没有一个被归零。
让我们用将被子（代表我们的回归模型）铺在床上（代表我们的数据）的例子来说明。在没有正则化的最小二乘法中，我们试图将被子铺得尽可能多地接触床上的点（数据点），最小化被子与床面之间的距离（误差）。
现在，想象一下床非常不平坦，而我们的被子非常柔韧。如果没有正则化，我们可能会把被子贴得非常紧，以至于它适应了每一个小凸起和凹陷，甚至是由于有人没有铺平床单的小问题——这就是过拟合。
引入正则化：
-
使用套索回归 (L1)： 就像是说：“我希望被子合适，但也希望它尽量平整，褶皱尽可能少。”每一个褶皱代表模型中的一个特征，而 L1 正则化尝试最小化褶皱的数量。最后，你会有一个合适的被子，但可能不会进入每一个缝隙，尤其是那些仅仅是噪声的缝隙。在几何视觉中，套索回归会试图让平面（colX）很好地贴合点，但可能会在不重要的变量方向上变得平坦（将系数缩小为零）。
-
使用岭回归 (L2)： 这就像是想要被子贴合得很好，但同时也希望被子均匀展开，没有任何部分离床太远。因此，即使被子仍然紧贴床面，它也不会过度扭曲以适应细微的凸起。在几何视觉中，岭回归添加了一个惩罚项，约束了系数，使它们缩小到接近零但不完全为零。这保持了平面接近点，但防止它倾斜得太厉害以适应特定点，从而保持了一定的距离（偏差）以防止对噪声的过拟合。
带有正则化的视觉解释：
当将正则化添加到几何解释中时：
-
平面（colX）可能不会像之前那样靠近每个单独的数据点（y）。正则化故意引入了一些偏差。
-
平面将趋向于更加稳定，较少倾斜于任何单个异常值点，因为对于大系数的惩罚意味着模型不能过于敏感于任何单一的数据点或特征。
-
向量（Beta1X1 和 Beta2X2）的长度可能会更短，反映了每个变量对预测的影响被故意限制的事实。
从本质上讲，正则化在让模型的训练数据拟合能力稍微打折的同时，换取了改进的模型泛化能力，这意味着模型在未见过的数据上表现更好，而不仅仅是在训练数据上。就像选择一个稍微宽松的毯子尺寸，它在所有实际用途上都足够好，而不是一个完美贴合每一个轮廓但可能不切实际或过于特定于一个床的尺寸。
结论
总结来说，线性回归的几何解释弥合了经典统计学与机器学习之间的差距，提供了对这一基本技术的直观理解。虽然经典统计学通过普通最小二乘法（Ordinary Least Squares）来处理线性回归，而机器学习则通常使用最大似然估计（Maximum Likelihood Estimation）或目标函数最小化，但这两种方法最终都旨在以视觉上易于理解的方式最小化预测误差。
引入像 Lasso 和 Ridge 这样的正则化技术进一步丰富了这种解释，突出了模型准确性和泛化能力之间的平衡。这些方法防止了过拟合，确保模型对新的、未见过的数据保持稳健和有效。
总体而言，这种几何视角不仅揭示了线性回归的奥秘，还强调了在不断发展的数据分析和机器学习领域中基础概念的重要性。这是一个强有力的提醒，表明复杂的算法可以根植于简单而深刻的几何原理中。
Python 中的地理空间数据分析
原文：
towardsdatascience.com/geospatial-data-analysis-in-python-d8fa8dd23a6c
使用 OSMnx 和 Kepler.gl 开始进行地理数据分析
·发表于Towards Data Science ·6 分钟阅读·2023 年 5 月 3 日
–
照片由Aaron Burden拍摄，来源于Unsplash
介绍
地理空间数据无处不在，并用于各种不同的应用场景（例如，根据位置计算房地产风险、设计新建筑开发、规划货物运输以及寻找不同地点之间的可能路线）。
地理空间数据通常以两种可能的格式存储：栅格和矢量：
-
栅格数据将数据表示为像素矩阵（因此具有固定分辨率）。在这种表示中，每个像素可以分配不同的值，多个栅格叠加在一起可以进一步增强同一图像。例如，同一图像可以使用 3 个通道/波段（例如 RGB——红色、绿色、蓝色）或使用单一通道进行存储。
-
矢量可以用于抽象现实世界的几何形状，使用点、线、多边形等元素，并且通常可以与表示它们的对象的一些有用元数据一起存储（例如名称、地址、所有者等）。由于它们作为数学对象存储，因此可以在不损失分辨率的情况下放大矢量。
矢量数据通常存储在 SVG 和 Shapefiles 等文件格式中，而栅格数据则存储在 TIFF、JPG、PNG 等格式中。
在处理地理空间数据时，通常需要进行许多不同形式的操作/转换，一些例子包括：
-
从非表格/原始二进制格式转换为矢量/栅格。
-
将连续数据分桶为离散类别。
-
从数据中提取多边形/特征。
-
处理无数据和离群值。
-
在不同坐标系统中重投影。
-
生成较低分辨率的数据概览，以处理不同的缩放级别，并避免图像渲染过载。
如果你对图像数据在机器学习系统中的具体应用感兴趣，可以在我的上一篇文章中找到更多信息。
演示
作为本文的一部分，我们现在将通过实际演示来分析地理空间矢量数据，以识别特定位置并计算到达该位置的最短路径。本文中使用的所有代码（以及更多！）可在我的 GitHub和Kaggle 账户上找到。
首先，我们需要确保在环境中安装了所有必要的依赖。
pip install osmnx
pip install folium
pip install keplergl
在这个例子中，我们将探索瑞士苏黎世的 Altstadt 区。利用OSMnx，只需使用两行代码就可以检索并可视化我们需要的数据（图 1）。OSMnx实际上是为了以最简单的方式提取和使用OpenStreetMap数据而设计的。OpenStreetMap 是由志愿者社区维护的全球免费地理数据库。
import osmnx as ox
import networkx as nx
import matplotlib.pyplot as plt
import folium
from keplergl import KeplerGl
place_name = "Altstadt, Zurich, Switzerland"
graph = ox.graph_from_place(place_name)
ox.plot_graph_folium(graph)
图 1：待检查区域的图表（图片来源：作者）。
我们现在可以深入探讨获取数据的不同方法。首先，我们可以获得我们所表示区域的多边形作为探索的基础。一旦获取了数据，它将以Geo Pandas数据框的形式呈现，其中包含关于该区域的所有关键信息（图 2）。Geo Pandas 是一个开源库，专门设计用于扩展 Pandas 的功能，以处理地理空间数据。
为了创建任何可视化表示，geometry列用作参考点，每一行都表示在地图上创建对象所需的所有坐标（例如多边形、线、点、多边形）。在这个例子中，Well-known text (WKT)被用作表示我们的矢量几何对象的文本标记语言，但通常也可以使用其他格式，如GeoJSON。此外，GeoSeries 中的每个值都存储在Shapely Object中，以便尽可能容易地执行操作和转换。
area = ox.geocode_to_gdf(place_name)
area.head()
图 2：GeoPandas 数据集示例（图片来源：作者）。
此时，我们可以重复相同的程序，以检索我们希望在地图上绘制的不同兴趣点。在这种情况下，假设我们正在苏黎世的 Altstadt 度假，我们想要查看所有前往餐厅的选项。为了简化我们的研究，我们可以获取所有代表街道的节点和边缘，以及所有可用的建筑物和餐馆，以便我们在地图上进行定位。
buildings = ox.geometries_from_place(place_name, tags={'building':True})
restaurants = ox.geometries_from_place(place_name,
tags={"amenity": "restaurant"})
nodes, edges = ox.graph_to_gdfs(graph)
一旦我们加载了所有数据，我们可以独立绘制每个特征（图 3）。
图 3：区域、建筑物、餐馆、节点和边缘表示（作者提供的图像）。
使用下面的code，这可以最终在下方的单一图表中漂亮地总结出来（图 4）。
fig, ax = plt.subplots(figsize=(10, 14))
area.plot(ax=ax, facecolor='white')
edges.plot(ax=ax, linewidth=1, edgecolor='blue')
buildings.plot(ax=ax, facecolor='yellow', alpha=0.7)
restaurants.plot(ax=ax, color='red', alpha=0.9, markersize=12)
plt.tight_layout()
图 4：分析摘要（作者提供的图像）。
为了使我们的分析更具互动性，我们还可以使用额外的库，例如KeplerGL。KeplerGL是由 Uber 开发的开源地理空间工具箱，用于创建高性能的基于网页的地理应用程序。
使用 KeplerGL，我们可以轻松地将地图叠加在真实的全球地图上，并实时应用各种转换和过滤（图 5）。
K_map = KeplerGl()
K_map.add_data(data=restaurants, name='Restaurants')
K_map.add_data(data=buildings, name='Buildings')
K_map.add_data(data=edges, name='Edges')
K_map.add_data(data=area, name='Area')
K_map.save_to_html()
图 5：带 KaplerGL 的互动分析摘要（作者提供的图像）。
现在我们已经构建了地图，并且有了一个互动工具来检查这些数据，我们可以最终尝试将研究范围缩小到一个餐厅。在这种情况下，我们首先将重点限制在意大利餐馆，然后选择 Antica Roma 作为我们的选择地点。
现在，我们只需要指定我们的初始位置坐标和目的地，以开始寻找最佳的步行路径（图 6）。
it_resturants = restaurants.loc[restaurants.cuisine.str.contains('italian')
.fillna(False)].dropna(axis=1, how='all')
resturant_choice = it_resturants[it_resturants['name'] == 'Antica Roma']
orig = list(graph)[993]
dest = ox.nearest_nodes(graph,
resturant_choice['geometry'][0].xy[0][0],
resturant_choice['geometry'][0].xy[1][0])
nodes[nodes.index == orig]
图 6：起始位置表示（作者提供的图像）。
为了完成这项任务，我们可以使用Networkx shortest_path function 来自动运行Dijkstra 算法 并尝试优化我们的路径以最小化其总体长度（图 7）。
route = nx.shortest_path(graph, orig, dest, weight='length')
ox.plot_route_folium(graph, route, route_linewidth=6, node_size=0)
图 7：到达目的地的最短路径（作者提供的图像）。
最后，为了验证我们的研究，我们可以尝试在 Google Maps 上执行相同的查询，如下所示，结果非常接近！（图 8）。
图 8：Google Maps 到达目的地的最短路径（作者提供的图像）。
结论
总结一下，在这篇文章中，我们首先介绍了地理空间数据如何在企业中使用，它通常是如何存储和处理的，并通过一个实际示例来识别两个不同点之间的最短路径。当然，我们也可以使用基于 UI 的工具进行类似的分析，比如OpenStreetMap 应用或Google Maps，虽然学习这些基础知识仍然非常有价值，因为它们可以应用于许多其他形式的网络问题（例如旅行推销员问题）。
联系方式
如果你想了解我最新的文章和项目，可以在 Medium 上关注我并订阅我的邮件列表。以下是我的一些联系方式：
使用 GeoPandas 进行地理空间数据分析
原文：
towardsdatascience.com/geospatial-data-analysis-with-geopandas-876cb72721cb
学习如何使用 Python 的 GeoPandas 操作和可视化矢量数据
·发布于Towards Data Science ·6 分钟阅读·2023 年 5 月 6 日
–
图片由Artem Beliaikin拍摄，发布于Unsplash
这是关于地理空间数据分析系列的第三篇文章：
-
使用 QGIS 进行地理空间数据分析
-
OpenStreetMap 入门指南
-
使用 GeoPandas 进行地理空间数据分析（本文）
-
使用 OSMnx 进行地理空间数据分析
本文是使用 QGIS 进行地理空间数据分析的实用介绍和 OpenStreetMap 入门指南的延续。在之前的教程中，我提供了地理空间数据分析的概述，这是一门无处不在且可以应用于许多领域的子领域，如物流、运输和保险。
这个学科专注于分析一种特殊类型的数据，即地理空间数据，其特点是具有位置，由一对或多对坐标描述。例子包括餐馆、道路和国家之间的边界。为了展示连续的表面，比如卫星图像，地理表格已不再足够，你需要一个包含一个或多个通道的数组。
在这篇文章中，我将重点讲解最简单的情况，即地理表格，也称为矢量数据。对于这个任务，GeoPandas 是将用于处理和可视化这种类型地理空间数据的 Python 库。正如你可能猜到的，它是 Pandas 的扩展，一个流行的 Python 包，允许你轻松且快速地处理地理空间数据。让我们开始吧！
目录：
-
导入普查数据
-
向普查数据添加几何信息
-
使用 GeoPandas 创建地图
-
从几何信息中提取质心
-
创建更复杂的地图
导入普查数据
开始地理空间数据分析之旅的最佳方法是通过处理普查数据，这些数据提供了世界各国所有人和家庭的详细信息。
在本教程中，我们将使用一个提供英国汽车或货车数量的数据集，该数据集来自于英国数据服务。数据集的链接是这里。
我将从一个不包含地理信息的数据集开始：
数据集的每一行对应于一个特定的输出区域，这是英国普查提供的最低地理级别。有三个特征：地理编码、国家以及一个或多个家庭成员拥有的汽车或货车数量。
如果我们现在想要可视化地图，我们将无法实现，因为我们没有必要的地理信息。在展示 GeoPandas 的潜力之前，我们需要进一步的步骤。
向普查数据添加几何信息
为了可视化我们的普查数据，我们需要添加一个存储地理信息的列。添加地理信息的过程，例如为每个城市添加纬度和经度，称为地理编码。
在这种情况下，不只是一个坐标对，而是有不同的坐标对，它们连接并闭合，形成输出区域的边界。我们需要从这个链接导出 Shapefile。它提供了每个输出区域的边界。
一旦数据集被导入，我们可以使用它们的公共字段 geo_code 合并这两个表：
在评估了数据框的维度在左连接后没有变化后，我们需要检查新列中是否存在空值：
df.geometry.isnull().sum()
# 0
幸运的是，没有空值，我们可以使用 GeoDataFrame 类将数据框转换为地理数据框，并将几何列设置为我们地理数据框的几何体：
现在，地理信息和非地理信息被合并到一个独特的表中。所有的地理信息都包含在一个字段中，称为 geometry。就像在普通的数据框中一样，我们可以打印出这个地理数据框的信息：
从输出中，我们可以看到我们的地理数据框是geopandas.GeoDataFrame对象的一个实例，几何体使用几何类型进行编码。为了更好地理解，我们还可以显示第一行中几何列的类型：
type(gdf.geometry[0])
# shapely.geometry.polygon.Polygon
了解几何对象中有三种常见的类是很重要的：点、线和多边形。在我们的例子中，我们处理的是多边形，这很有意义，因为它们是输出区域的边界。然后，数据集已经准备好了，我们可以从现在开始构建漂亮的可视化图。
使用 GeoPandas 创建地图
现在，我们有了使用 GeoPandas 可视化地图的所有必要要素。由于 GeoPandas 的一个缺点是它处理大量数据时存在困难，而我们有超过 20 万行数据，因此我们将专注于北爱尔兰的人口普查数据：
gdf_ni = gdf.query(‘Country==”Northen Ireland”’)
要创建地图，只需在地理数据框上调用plot()方法：
我们还希望通过根据每个输出区域的频率进行着色，查看汽车/面包车在北爱尔兰的分布情况：
从这个图中，我们可以观察到大多数区域的车辆数量在 200 辆左右，除了标记为绿色的小区域。
从几何体中提取质心
假设我们想要更改几何形状，并将坐标放在输出区域的中心，而不是多边形。可以通过使用gdf.geomtry.centroid属性来计算每个输出区域的质心来实现这一点：
gdf_ni[‘centroid’] = gdf.geometry.centroid
gdf_ni.sample(3)
如果我们再次显示数据框的信息，我们可以注意到几何体和质心都被编码为几何类型。
更好地理解我们真正获得了什么的是在一个唯一的地图中可视化几何体和质心列。要绘制质心，需要使用set_geometry()方法来切换几何体。
创建更复杂的地图
有一些高级功能可以在地图中可视化更多细节，而无需创建其他信息列。之前我们显示了每个输出区域的汽车或面包车数量，但这比信息量更混乱。基于我们的数值列创建一个类别特征会更好。使用 GeoPandas，我们可以跳过那一阶段，直接进行绘图。通过指定参数scheme=’intervals’，我们能够根据等间隔创建汽车/面包车的类别。
地图没有发生太大变化，但你可以看到图例相比之前的版本更清晰。更好的可视化地图的方法是基于分位数构建的级别来着色：
现在，在地图上可以发现更多的变化，因为每个级别包含了更多分布的区域。值得注意的是，大多数区域属于最后两个级别，对应于车辆数量最多的情况。在第一次可视化中，200 辆车似乎数量不多，但实际上存在许多高频率的离群值，这扭曲了我们的解释。
此时，我们还希望有一个背景地图来更好地使我们的结果具有上下文。最流行的方法是使用 contextily 库，它可以获取背景地图。该库需要 Web Mercator 坐标参考系统（EPSG:3857）。因此，我们需要将数据转换为该坐标参考系统。绘制地图的代码保持不变，只需增加一行代码来添加来自 Contextily 库的基础地图：
太好了！现在，我们有了一张更专业、更详细的地图！
最后思考：
这是一个入门教程，用于开始用 Python 处理地理空间数据。GeoPandas 是一个专门处理矢量数据的 Python 库。它使用起来非常简单直观，因为它具有类似于 Pandas 的属性和方法，但一旦数据量增加，特别是在绘制数据时，它会变得非常缓慢。
除了这个缺点外，还有一个事实是它依赖于 Fiona 库来读取和写入矢量数据格式。如果 Fiona 不支持某些格式，即使是 GeoPandas 也无法支持它们。一个解决方案是结合使用 GeoPandas 来处理数据和 QGIS 来可视化地图。或者尝试其他 Python 库来可视化数据，如 Folium。你知道其他的替代方案吗？如果有其他想法，请在评论中提出。
代码可以在 这里 找到。希望你觉得这篇文章有用。祝你有美好的一天！
免责声明：数据集的许可为英国开放政府许可证（OGL）
有用的资源：
使用 OSMnx 进行地理空间数据分析
原文：
towardsdatascience.com/geospatial-data-analysis-with-osmnx-8a300d77b592
学习如何使用 Python 下载、分析和可视化 OpenStreetMap 数据
·发布于 Towards Data Science ·6 分钟阅读·2023 年 6 月 16 日
–
照片由 Denys Nevozhai 提供，来源于 Unsplash
这是关于地理空间数据分析系列的第四篇文章：
-
使用 QGIS 进行地理空间数据分析
-
OpenStreetMap 入门指南
-
使用 GeoPandas 进行地理空间数据分析
-
使用 OSMnx 进行地理空间数据分析（本文）
在之前的教程中，我涵盖了地理空间数据分析的各个方面。我首先展示了不使用任何代码的地理空间数据实际示例，以帮助你深入理解这些概念。地理空间数据分析是一个普遍存在的领域，涉及处理一种特殊类型的数据——地理空间数据。
这包括将位置添加到非地理数据中。文中充满了示例。你可以想到咖啡馆、医院、道路、河流、卫星图像等。即使你用 Google Maps 搜索一个地点，你也在与地理空间数据互动。
这次我将专注于从 OpenStreetMap 下载、可视化和分析数据。OpenStreetMap 是最大的免费且可编辑的地理数据库，来自世界各地的志愿者共同参与这个庞大的项目。这个教程得益于一个名为 OSMnx 的 Python 包。让我们开始吧！
目录：
-
OSMNx 简介
-
下载和可视化 OSM 数据
-
将图形转换为 GeoDataframe
-
提取兴趣点
-
寻找最短路径
OSMNx 简介
OSMnx 是一个用于下载、分析和可视化 OpenStreetMap 网络数据的库。它依赖于两个库，NetworkX和GeoPandas。特别地，它利用 NetworkX 库中的graph模块来检索网络数据。
此外，它允许我们与两个 OpenStreetMap API 进行交互：
-
Nomatism用于地理编码，即通过名称和地址查找位置。
-
Overpass API用于提取兴趣点，如高速公路、学校和公园。
下载和可视化 OSM 数据
我们可以直接使用 OSMNx，而不是从网站或 Geofabrik 手动下载数据。
首先，我们需要导入四个库，这些库将在教程中稍后使用：
import osmnx as ox
import folium
import contextily as cx
import matplotlib.pyplot as plt
除了 OSMnx 和 matplotlib，我们还将利用folium，它以创建交互式地图而闻名，以及contextily来添加背景地图。这一点对获得现实地图非常重要。
像之前一样，我们将读取并可视化博洛尼亚的 OSM 街道网络数据，博洛尼亚是意大利最大的城市之一。
PLACE_NAME = 'Bologna, Italy'
G = ox.graph_from_place(PLACE_NAME, network_type='drive')
ox.plot_graph(G)
博洛尼亚的黑白可视化
从黑白可视化中，我们可以观察到代表节点的点和描绘线路的线条。与 OpenStreetMap 网站相比，它可能显得非常静态和基础。Folium 通过其强大的可解释地图来拯救我们：
ox.plot_graph_folium(G)
使用 Folium 获得的博洛尼亚地图
这要好得多，你不觉得吗？亮丽的颜色和与地图交互的可能性是使用 Google Maps 去未知地方时至关重要的特点。
如果你更仔细地查看 OpenStreetMap 网站，你会发现默认情况下是标准图层。除了标准图层，还有其他图层，比如 Cycle Map 和 Transport Map。我们可以根据目的利用不同的图层，这一点非常令人惊叹。
如果我们对自行车情有独钟，我们会更关注 Cycle Map。只需一行代码就能实现：
G = ox.graph_from_place(PLACE_NAME, network_type='bike')
在接下来的部分，我们将考虑标准图形。
将图形转换为 GeoDataframe
处理图形不像处理 Dataframes 和 GeoDataframes 那样直观。因此，我们可能需要将图形转换为 GeoDataframe：
area,edges = ox.graph_to_gdfs(G)
area.head()
edges.head()
你会发现我们得到了两个 GeoDataframe，一个用于节点，一个用于边。查看几何形状就会很清楚。区域 geodataframe 只有一对坐标，即纬度和经度，而包含边的 geodataframe 中有两对坐标。
提取兴趣点
在进行数据科学项目时，我们尝试通过在互联网上搜索开放数据来为数据集添加信息。从 OSM 数据中，可以提取 兴趣点（POI），这些是根据我们的分析目的我们可能会觉得有趣的地方。示例包括餐馆、教堂、博物馆和公园。
例如，我们想要分析博洛尼亚的交通情况，以优化和减少运输成本。在这种情况下，了解高速公路、加油站、停车场以及其他可能成为瓶颈的地方将非常有用。
让我们来看看城市中的所有加油站。这可以通过将 燃料 作为 amenity 键的值来实现。
fuel_stations = ox.geometries_from_place(
PLACE_NAME,
{"amenity": "fuel"},
)
fuel_stations.head()
由于我们已经提取了所有加油站，通过可视化地图来理解它们的位置会更有用。此外，我们可以添加基础地图来更好地 contextualize 我们的结果。
area_crs = area.to_crs('3857')
edges_crs = edges.to_crs('3857')
fuel_stations_crs = fuel_stations.to_crs('3857')
fig, ax = plt.subplots(figsize=(10, 14))
area_crs.plot(ax=ax, facecolor='white')
edges_crs.plot(ax=ax, linewidth=1, edgecolor='blue')
fuel_stations_crs.plot(ax=ax, color='red', alpha=0.9, markersize=12)
plt.tight_layout()
cx.add_basemap(ax,crs=area_crs.crs.to_string())
太棒了！我们可以注意到大多数加油站集中在城市边缘。此外，我们可以区分不同组的服务站，这些在测量中心外交通时需要考虑。
找到最短路线
OSMnx 库的另一个有用功能是能够计算两个点之间的最短路径。
origin = (
ox.geocode_to_gdf("Parco della Montagnola, Bologna, Italy")
.to_crs(edges.crs)
.at[0, "geometry"]
.centroid
)
destination = (
ox.geocode_to_gdf("Esso, Bologna, Italy")
.to_crs(edges.crs)
.at[0, "geometry"]
.centroid
)
origin_node_id = ox.nearest_nodes(G, origin.x, origin.y)
destination_node_id = ox.nearest_nodes(G, destination.x, destination.y)
这是通过 shortest_path() 方法实现的，该方法默认使用 Dijkstra 算法 来计算源节点和目标节点之间的路径。
route = ox.shortest_path(G, origin_node_id, destination_node_id)
route
#[400881920,
# 250763178,
# 250763179,
# 250763533, ...
# 1694666466]
我们还可以尝试在同一地图上可视化图形和最短路径：
ox.plot_route_folium(G, route, route_linewidth=6, node_size=0)
看吧！这就像我们使用了 Google Maps 来寻找路线一样，不过我们利用了 OSMnx 库的功能来进行查找。
最后的想法：
这是一份指南，让你了解如何使用 Python 操作 OSM 数据。我发现 OSMnx 是处理 OpenStreetMap 数据的最完整的 Python 库。当然，它更适合探索较小的地方，比如城市。如果数据集较大，最好使用更专业的软件，比如 QGIS 来进行可视化。如果你尝试过其他库来处理 OSM 数据，请在评论中分享。如果你想查看代码，可以点击 这里。感谢阅读！祝你有美好的一天！
有用的资源：
你喜欢我的文章吗？成为会员 并每天获得无限的数据科学文章访问权限！这是一种间接支持我的方式，对你没有额外费用。如果你已经是会员， 订阅 以便每当我发布新的数据科学和 Python 指南时接收邮件！
地理空间数据工程：空间索引
优化查询、提高运行时间和地理空间数据科学应用
·
关注 在 Towards Data Science ·6 分钟阅读·2023 年 8 月 31 日
–
图片由 Tamas Tuzes-Katai 提供，来源于 Unsplash
介绍：空间索引为何有用？
在进行地理空间数据科学工作时，优化你编写的代码是非常重要的。如何让拥有数亿行的数据库集合或联接更快？这就是空间索引等概念发挥作用的地方。在这篇文章中，我将讨论空间索引的实现方式，它的好处和局限性，并看看 Uber 的开源 H3 索引库在一些有趣的空间数据科学应用中的表现。我们开始吧！
🗺 什么是空间索引？
常规索引就像你在书籍末尾找到的内容：一个单词的列表以及这些单词在文本中出现的位置。这帮助你快速查找在特定文本中你感兴趣的单词的任何参考。如果没有这个方便的工具，你将需要手动翻阅书中的每一页，寻找你想阅读的那一个提及。
在现代数据库中，查询和搜索的问题也非常相关。索引通常比过滤更快地查找数据，你可以基于感兴趣的列创建索引。对于地理空间数据，工程师们通常需要考虑诸如“交集”或“附近”的操作。我们如何制作一个空间索引，使这些操作尽可能快速？首先，让我们看看一些地理空间数据：
两个不相交的特征（图片来源：作者）
假设我们想运行一个查询来确定这两个形状是否相交。根据构造，空间数据库创建的索引是由包含几何图形的边界框组成的：
制作一个大的边界框（图片来源：作者）
为了回答这两个特征是否相交，数据库将比较两个边界框是否有任何公共区域。如你所见，这可能会迅速导致假阳性。为了解决这个问题，像 PostGIS 这样的空间数据库通常会将这些大的边界框划分成越来越小的部分：
更小的尺寸：制作子边界框（图片来源：作者）
这些分区存储在 R 树中。R 树是一种层次数据结构：它跟踪大的“父”边界框、它的子节点、子节点的子节点，等等。每个父节点的边界框包含它子节点的边界框：
R 树可视化（图片来源：作者）
操作“交集”是从这种结构中受益的关键操作之一。在查询交集时，数据库会沿着这棵树向下查询，问“当前的边界框是否与感兴趣的特征相交？”。如果是，它会查看该边界框的子节点，并提出相同的问题。通过这种方式，它能够快速遍历树，跳过没有交集的分支，从而提高查询性能。最后，它会返回所需的交集几何图形。
🧰 实践中：尝试使用 GeoPandas 的空间索引
现在让我们具体看看使用常规行级过程与空间索引的不同。我将使用代表纽约市人口普查区和市设施的两个数据集（均通过开放数据许可提供，可在这里和这里获得）。首先，让我们在 GeoPandas 中尝试对一个人口普查区几何图形进行“交集”操作。在 GeoPandas 中，“交集”是一个逐行函数，检查每一行的列是否与我们的几何图形相交。
GeoPandas 还提供了一种使用 R 树的空间索引操作，允许我们执行交集操作。以下是对这两种方法在 100 次交集操作中的运行时比较（注意：由于默认的交集函数较慢，我只选择了原始数据集中的大约 100 个几何图形）：
💨 空间索引快了多少？（图片由作者提供）
如你所见，空间索引方法相比普通的交集方法提供了显著的性能提升。实际上，这里是每种方法运行时间的 95% 置信区间：
太好了！那么，为什么我们有时不想使用空间索引呢？是否有些情况下它没有任何好处？确实有。这些限制中的一些是由于空间索引存储数据叶子的方式。结果表明，原始数据的分布方式会影响边界框在 R 树中的位置。具体来说，如果大量数据集中在同一地理空间，它们往往会共享相同的父节点，从而被分组在同一分支中。这可能导致树的不均衡，从而在查询时无法提供太多优化。
💻 其他空间索引是什么样的？
其他公司也适应了自己的空间索引。Uber 使用 H3，这是一个将世界划分为等面积六边形的分层索引系统。六边形在建模人们在城市中的移动或解决如计算半径等问题时具有许多优点。地理空间数据被划分到这些六边形中，作为公司的主要分析单位。网格是通过将 122 个六边形单元覆盖到二十面体地图投影上来构建的，支持包括聚合、连接和机器学习应用在内的广泛功能。
这个系统以及它的许多功能是开源的，可以在 GitHub 上进行分析。H3 API 的一个功能是将经纬度点转换为表示唯一六边形的字符串，根据指定的分辨率。我们来对整个设施数据库进行这个操作，并将六边形字符串转换为多边形：
在这些空间数据分析项目中，一个常见的问题是六边形中的项目有多少被某个列，比如“机构”所分类。幸运的是，现在我们已经将数据划分到 H3 六边形中，这很容易计算和可视化：
哪些机构有最多的设施？（图片来源于作者）
在这种情况下，你可以看到 DCAS（城市行政服务部）和 PARKS（公园与娱乐部）是每个六边形中设施最多的两个机构。这很有道理，因为这两个机构通常会有更多的实体设施（比如行政建筑或像公园这样的休闲区域）。
结论
如你所见，空间索引是地理空间数据科学和分析中非常有用的优化工具。在简单的交集查询中，使用空间索引显著提高了查询性能，相比于标准的 GeoPandas 交集函数。这种索引的实现和它的影响有许多细微之处，比如拥有巨大的数据集群分支。如我们所见，各公司已经开发了自己的解决方案：一个例子是 Uber 的 H3 开源索引，它允许我们回答各种空间分析问题。虽然我演示了基于机构的设施计数操作，但 H3 为其他更复杂的机器学习应用提供了基准。
如果你喜欢这类内容，但想要更广泛地了解城市规划技术，我也写了一份名为“The Zoned Out Chronicles”的通讯。我鼓励你去看看！
感谢阅读！
地理空间数据科学：点模式分析
原文：
towardsdatascience.com/geospatial-data-science-points-pattern-analysis-a61e04a6ddb8
在 Python 中执行地理空间点模式分析的快速教程。
·发表于 Towards Data Science ·阅读时间 8 分钟·2023 年 9 月 25 日
–
图片由 Bernard Hermant 在 Unsplash 提供
介绍
地理空间数据科学是数据领域的一个子领域，处理考虑事件发生在空间中的位置的数据点分析。
假设我们拥有一家销售智能手机的零售连锁店。我们的连锁店有几个配送中心，我们即将开设几家新店。我们可以在哪里开设这些新店？
这样的洞察可能来自地理空间分析，这种分析会向我们展示销售的集中地点，如果销售在某些地方更高或更低，以及其他洞察。
当我们想确保查看的是一个地理上集中的数据集时，点模式分析就会进入这个游戏。就像我们数据科学家的工作一样，点模式涉及到创建一个假设并通过应用于数据的统计方法来消除许多不确定性，以确认这个假设。在这种情况下，也不例外。这里有几个统计测试需要完成，本文将展示这些测试。
顺便提一下，我们最近在我的博客中研究了地理空间数据科学。如果你对这个主题了解不多，这里有两篇不错的读物可以在你深入阅读之前参考。
一个包含 Python 代码的实用数据分析帖子。
[towardsdatascience.com ## 使用 Python 分析地理空间数据（第二部分 — 假设检验）
学习阿什维尔 AirBnb 房源的地理空间假设检验。
[towardsdatascience.com
编码
包
让我们从这次练习中使用的包开始。如果你的环境中没有安装它们，别忘了使用pip install或conda install（对于 Anaconda 用户），然后跟上包名。
import pandas as pd
import numpy as np
import geopandas as gpd
import seaborn as sns
import matplotlib.pyplot as plt
import contextily
# Spatial Stats
from pointpats import distance_statistics, QStatistic, random, PointPattern
数据集
使用的数据集再次是AirBnb在美国北卡罗来纳州阿什维尔市的房源数据。数据可以从网站 insideairbnb.com/ 的一个独立项目中获取，任何人都可以前往下载数据集进行分析。这些数据在 创意共享署名 4.0 国际许可协议 下开放。
我已经下载了文件 listings.csv.gz。
要将数据加载到 Python 会话中，这里是代码。第一个代码片段是来自 Pandas 的简单read_csv()，我们预先确定了希望从原始数据中提取哪些列。然后，我们使用gpd.GeoDataFrame将数据集转换为Geopandas对象类型，指定用作 X 轴和 Y 轴的列，以及地理坐标系统（crs — 使用 4326，它与 GPS 系统相同，是最常见的参考系统之一）。
# Import the file to this exercise
# Open listings file
listings = pd.read_csv('/content/listings.csv',
usecols=['id', 'property_type', 'neighbourhood_cleansed',
'bedrooms', 'beds', 'bathrooms_text', 'price',
'latitude','longitude'])
# Convert the file to GeoPandas
points_gpd = gpd.GeoDataFrame(listings,
geometry= gpd.points_from_xy(
x=listings.longitude,
y=listings.latitude),
crs= "EPSG:4326")
基本上，将数据转换为 Geopandas 就是创建这个列geometry和对象类型。
Geopandas 数据集。图片由作者提供。
很好。完成这些后，我们来快速看看地图。代码很简单，创建了一个图形fig和两个坐标轴ax，因为会有两个图：一个是底图，另一个是点图。
#Quick check gpd dataframe
fig, ax = plt.subplots(figsize=(8,8))
# zorder=1 is the plot below
asheville.plot(ax=ax, color=None, zorder=1)
# zorder=2 is the top layer
points_gpd.plot(ax=ax, zorder=2, color='black', markersize=8)
这段代码生成了下一个图表。
阿什维尔，北卡罗来纳州的 Airbnb 房源列表。图片由作者提供。
很酷。地图看起来不错，但没有提供太多信息。现在让我们增强分析。
点模式
当考虑点模式分析时，我们可以做的第一个分析是检查这些数据点在地理上的集中程度。
这里，seaborn库可以提供帮助。方法jointplot将散点图与边缘的直方图结合起来。这是分析中的一个绝妙补充，因为它通过查看条形图较高的地方，给我们提供了点的集中情况和位置的见解。
要创建它，我们可以传递来自geometry列的 x 和 y，加上数据集（data）、点的大小（s）、color和图形的height。下一个片段是将基图添加到此联合图中，可以使用contextily，将联合图变量传递给该方法。
# Check concentration of the points
plot2 = sns.jointplot(
x= points_gpd.geometry.x,
y= points_gpd.geometry.y,
data= points_gpd,
s=5, height=7, color='k')
# Add a basemap to the jointplot
contextily.add_basemap( plot2.ax_joint,
crs="EPSG:4326",
source= contextily.providers.Stamen.TonerLite)
结果是我们看到了这个漂亮的地图。
在北卡罗来纳州阿什维尔的基础地图上的联合图。图片由作者提供。
我喜欢这个图。通过它，我们可以快速获得一些有价值的见解。用红色标记的市中心区域（及其周边）确实是列表集中所在的地方。然后我们还可以注意到，随着离该区域越来越远，列表数量减少并变得更加稀疏。
这完全有意义。我们来想一想：Airbnb 是一个让人们租赁房屋或卧室的平台。这些租赁房产通常位于住宅区，因为它们的主要用途是作为居住场所，而非商业场所。而房屋社区的地点通常位于城市基础设施周围的区域，比如商场、超市、药店、银行等。由于山区偏远地区的房屋建设更具挑战，因此预计那里也会有较少的位置点。
统计检验
现在我们已经绘制了联合图并获得了一些有用的见解，但我们仍然需要测试点的分布模式，以确定它们是否在统计上存在聚类。我是说，这些点也可能只是偶然聚集在一起的。谁知道呢。
所以，为了测试并确保你在处理的是一个模式，有两个很好的检验：
-
Ripley 的 G： 这个检验会检查一个点到其最近邻的距离的累积分布。因此，该检验测量从给定房屋到邻居 1、2、3、……、n 的距离，并将这些距离的分布与模拟的随机点分布进行比较。如果我们观察到观察到的数据在特定空间中的行为与模拟不同，我们可以得出结论：数据有一定模式，因此是聚类的。
-
Ripley 的 K： 这个检验进行类似的测试，将观察到的数据与随机模拟分布进行比较。与 G 检验的不同之处在于，K 检验考虑了数据中的所有距离，而不仅仅是最近的邻居。
Ripley 的 G 检验距离最近邻的分布。Ripley 的 K 检验整个数据集的距离分布。
Ripley 的 G
现在我们来学习如何执行这些检验。首先是 G 检验。这个检验写起来相对简单，但运行时间可能较长，取决于数据集的大小。在这里，检验获取了 40 个邻居的距离。运行大约花了 6 分钟。
# Coding Ripley's G (6 mins to run)
ripley_g = distance_statistics.g_test(points_gpd[['longitude', 'latitude']].values,
support=40,
keep_simulations= True)
要绘制结果，代码片段如下。我们为每次模拟的中位数绘制一条黑线，为每个数据点绘制一条红线与统计数据。
# Plot G test
plt.figure(figsize=(20,7))
# Simulated Data line plot
plt.plot(ripley_g.support,
np.median(ripley_g.simulations, axis=0),
color='k', label= 'Randomly Simulated Data')
# Ripley Stat plot for Observed data
plt.plot(ripley_g.support,
ripley_g.statistic, marker='o',
color='red', label= 'Observed Data')
# Plot setup
plt.legend(loc=4)
plt.xlabel('Distance')
plt.xticks( np.arange(0.0, 0.023, 0.001) )
plt.ylabel('Ripley G function statistic')
plt.title('Rypleys G Test')
plt.show()
结果是，代码显示了下一个图形。
Ripley 的 G 检验。图片由作者提供。
我们可以看到，对于距离 0 和 0.003，观察到的数据增长速度快于模拟数据，确认了数据集中存在显著的空间模式。
Ripley 的 K
第二个检验是 K 检验。这个检验使用了完全空间随机性的原假设和数据中存在空间模式的备择假设。
运行 K 检验和其他检验一样简单。不过，该检验的计算开销更大。在 Google Colab 上运行一个会话花费了超过 30 分钟。
# Coding Ripley's K
ripley_k = distance_statistics.k_test(points_gpd[['longitude', 'latitude']].values,
keep_simulations= True)
# Plot K test
# Simulated Data line plot
plt.plot(ripley_k.support,
ripley_k.simulations.T,
color='k', alpha=.1)
# Ripley Stat plot for Observed data
plt.plot(ripley_k.support,
ripley_k.statistic, marker='x',
color='orange')
# p<0.05 = alternative hypothesis: spatial pattern
plt.scatter(ripley_k.support,
ripley_k.statistic,
cmap='viridis', c=ripley_k.pvalue < .05,
zorder=4)
# Plot setup
plt.xlabel('Distance')
plt.ylabel('Ripley K function')
plt.title('Rypleys K Function Plot')
plt.show()
上面的代码将给我们下一个图表。
Ripley 的 K 检验。图片由作者提供。
再次强调，观察到的数据远高于模拟数据，确认了我们的空间模式。
在你离开之前
在这篇文章中，我们学习了如何为地理空间点模式确认创建一些统计检验。
一旦我们将一些数据点绘制在地图上，它们可能仅仅是偶然地被聚集在一起，而没有明确的模式。但确认地理空间模式的好方法是运行Ripley 的 G和Ripley 的 K检验。
这些检验可以帮助确认在例如北卡罗来纳州阿什维尔市的高价和低价簇中是否存在模式。因此，如果你要出租物业，你将知道在该市每个街区的最佳竞争价格。
现在，你可以在我的 GitHub 仓库中下载代码并应用到你的数据中。
## Studying/Python/Geospatial/Points_Pattern_Python.ipynb at master · gurezende/Studying
这是一个包含我对新软件包的测试和研究的仓库 - Studying/Python/Geospatial/Points_Pattern_Python.ipynb…
如果你喜欢我的内容，不要忘记关注我或在LinkedIn上找到我。
在 Medium 上阅读 Gustavo Santos 的文章。他是一名数据科学家。我从数据中提取见解，以帮助个人和公司…
参考文献
JORDAN, David S. [2023]. Applied Geospatial Data Science with Python. 1 ed. Pactk Publishing.
Ripley 函数是表征点模式空间结构的经典工具。它被广泛应用于…
www.hindawi.com](https://www.hindawi.com/journals/isrn/2013/753475/?source=post_page-----a61e04a6ddb8--------------------------------) [## pointpats/notebooks/distance_statistics-numpy-oriented.ipynb at main · pysal/pointpats
PySAL 中的平面点模式分析。通过在 GitHub 上创建帐户来贡献于 pysal/pointpats 的开发。
github.com](https://github.com/pysal/pointpats/blob/main/notebooks/distance_statistics-numpy-oriented.ipynb?source=post_page-----a61e04a6ddb8--------------------------------)
地理空间索引 102
一个实际示例，展示如何应用地理空间索引
·
关注 发表在 Towards Data Science · 7 分钟阅读 · 2023 年 4 月 11 日
–
介绍
地理空间索引是一种索引技术，它提供了一种优雅的方式来管理基于位置的数据。它使得地理空间数据可以高效地搜索和检索，从而系统能够为用户提供最佳体验。本文将通过将地理空间索引应用于实际数据，并展示这样做所带来的性能提升，来演示其实际操作。让我们开始吧。（注意：如果你从未听说过地理空间索引或希望了解更多信息，可以查看 这篇文章)
数据
本文使用的数据是芝加哥犯罪数据，这是Google Cloud 公共数据集计划的一部分。任何拥有 Google Cloud Platform 账户的人都可以免费访问此数据集。数据集包含约 800 万行数据（总量为 1.52 GB），记录了自 2001 年以来发生在芝加哥的犯罪事件，每条记录都有指示事件位置的地理数据。
平台
我们不仅会使用 Google Cloud 的数据，还会使用 Google Big Query 作为数据处理平台。Big Query 提供了每个执行查询的作业执行细节。这包括使用的数据量和处理的行数，这将非常有助于展示优化后的性能提升。
实现
我们将通过优化基于位置的查询来展示地理空间索引的强大功能。在这个例子中，我们将使用Geohash作为索引，因为它的简洁性和 Google BigQuery 的原生支持。
我们将检索发生在距离芝加哥联合车站2 公里范围内的所有犯罪记录。在优化之前，让我们看看在原始数据集上运行此查询时的性能表现：
-- Chicago Union Station Coordinates = (-87.6402895591744 41.87887332682509)
SELECT
*
FROM
`bigquery-public-data.chicago_crime.crime`
WHERE
ST_DISTANCE(ST_GEOGPOINT(longitude, latitude), ST_GEOGFROMTEXT("POINT(-87.6402895591744 41.87887332682509)")) <= 2000
以下是作业信息和执行细节的样子：
作业信息（作者提供的图片）
执行细节（作者提供的图片）
从处理的字节数和读取的记录数中，你可以看到查询扫描了整个表并处理了每一行，以获取最终结果。这意味着数据越多，查询所需时间越长，处理成本也越高。这可以更高效吗？当然可以，这就是地理空间索引发挥作用的地方。
上述查询的问题在于，尽管许多记录距离关注点（芝加哥联合车站）较远，但仍然必须处理这些记录。如果我们能排除这些记录，那么查询将更高效。
Geohash 可以是解决这个问题的方案。除了将坐标编码为文本外，geohash 的另一个优点是哈希还包含地理空间属性。哈希之间的相似性可以推测它们所代表位置的地理相似性。例如，wxcgh 和 wxcgd 代表的两个区域相近，因为这两个哈希非常相似，而 accgh 和 dydgh 之间则相距较远，因为这两个哈希非常不同。
我们可以利用这一属性与 聚簇表 进行配合，通过预先计算每一行的 geohash。然后，我们计算芝加哥联合车站的 geohash。这样，我们可以在预先筛选掉所有与芝加哥联合车站的 geohash 不够 接近 的记录。
实现方法如下：
- 创建一个新表，并添加一个新列来存储坐标的 geohash。
CREATE TABLE `<project_id>.<dataset>.crime_with_geohash_lv5` AS (
SELECT *, ST_GEOHASH(ST_GEOGPOINT(longitude, latitude), 5) as geohash
FROM `bigquery-public-data.chicago_crime.crime`
)
2. 使用 geohash 列作为聚簇键创建聚簇表
CREATE TABLE `<project_id>.<dataset>.crime_with_geohash_lv5_clustered`
CLUSTER BY geohash
AS (
SELECT *
FROM `<project_id>.<dataset>.crime_with_geohash_lv5`
)
通过使用 geohash 作为聚簇键，我们创建了一个表，其中共享相同哈希的行被物理地存储在一起。如果你考虑一下，实际发生的情况是数据集按地理位置进行分区，因为地理位置越接近的行更有可能具有相同的哈希。
3. 计算芝加哥联合车站的 geohash。
在这篇文章中，我们使用了 这个网站，但有很多不同编程语言的库可以让你以编程方式实现这一功能。
芝加哥联合车站的 geohash（作者提供的图像）
4. 将 geohash 添加到查询条件中。
SELECT
*
FROM
`<project_id>.<dataset>.crime_with_geohash_lv5_clustered`
WHERE
geohash = "dp3wj" AND
ST_DISTANCE(ST_GEOGPOINT(longitude, latitude), ST_GEOGFROMTEXT("POINT(-87.6402895591744 41.87887332682509)")) <= 2000
这次查询应该只扫描位于 dp3wj 的记录，因为 geohash 是表的聚簇键。这应该能节省很多处理时间。让我们看看会发生什么。
创建聚簇表后的作业信息（作者提供的图像）
创建聚簇表后的执行细节（作者提供的图像）
从作业信息和执行细节中，你可以看到处理的字节数和扫描的记录数显著减少（从 1.5 GB 减少到 55 MB，7M 减少到 260k）。通过引入 geohash 列并将其用作聚簇键，我们可以仅通过查看一列就提前排除所有明显不满足查询的记录。
但是，我们还没有完成。仔细查看输出的行数，你会发现它只有 100k 条记录，而正确的结果应该有 380k 条。我们得到的结果仍然不正确。
5. 计算邻近区域并将它们添加到查询中。
在这个例子中，所有的邻近哈希是 dp3wk、dp3wm、dp3wq、dp3wh、dp3wn、dp3wu、dp3wv 和 dp3wy。我们使用在线的 geohash 探索器 来查看这些，但这也可以完全通过代码实现。
dp3wj 的邻居（作者提供的图像）
为什么我们需要在查询中添加邻域区域？因为地理哈希只是位置的一种近似。虽然我们知道芝加哥联合车站在dp3wj区域内，但我们仍然不知道它在区域中的具体位置。是在顶部、底部、左侧还是右侧？我们无从得知。如果它在顶部，那么在dp3wm中的一些数据可能比 2 公里更接近。如果它在右侧，则dp3wn区域中的一些数据可能距离更近。依此类推。这就是为什么所有邻域哈希必须包含在查询中以获得正确结果的原因。
请注意，地理哈希级别 5 的精度为 5 公里。因此，除上图中的区域外，所有其他区域距离芝加哥联合车站都太远。这是另一个需要做出的重要设计选择，因为它有巨大的影响。如果精度过粗，我们得到的收益很少。另一方面，使用过细的精度级别会使查询变得复杂。
这是最终查询的样子：
SELECT
*
FROM
`<project_id>.<dataset>.crime_with_geohash_lv5_clustered`
WHERE
(
geohash = "dp3wk" OR
geohash = "dp3wm" OR
geohash = "dp3wq" OR
geohash = "dp3wh" OR
geohash = "dp3wj" OR
geohash = "dp3wn" OR
geohash = "dp3tu" OR
geohash = "dp3tv" OR
geohash = "dp3ty"
) AND
ST_DISTANCE(ST_GEOGPOINT(longitude, latitude), ST_GEOGFROMTEXT("POINT(-87.6402895591744 41.87887332682509)")) <= 2000
这就是执行查询时发生的情况：
添加邻域哈希后的工作信息（图片由作者提供）
添加邻域哈希后的执行细节（图片由作者提供）
现在结果是正确的，查询处理了 527 MB 的数据，总共扫描了 250 万条记录。与原始查询相比，使用地理哈希和聚簇表节省了大约 3 倍的处理资源。然而，任何东西都不是免费的。应用地理哈希增加了数据预处理和检索的复杂性，例如需要提前选择的精度级别和额外的 SQL 查询逻辑。
结论
在本文中，我们已经看到地理空间索引如何帮助改善地理空间数据的处理。然而，它有一定的成本，需要提前仔细考虑。归根结底，这不是免费的午餐。要使其正常工作，需要对算法和系统要求有良好的理解。
最初发布于 https://thanakornp.com。
地理空间索引和评分：释放基于位置的数据分析的力量
使用 Python 和 H3 进行地理空间索引的实用指南
·发表于Towards Data Science ·阅读时长 7 分钟·2023 年 1 月 25 日
–
图片由Antoine Merour提供，来源于Unsplash
在本文中，我们将讨论地理空间索引和评分，以及如何使其对开发人员变得简单。地理空间索引是创建一个与地理位置（如纬度和经度）相关的数据集的索引的过程。评分则是根据特定标准对数据进行排名或排序的过程。地理空间索引和评分结合在一起，可以创建强大且高效的应用程序，能够快速检索和显示地理数据。
地理空间索引和评分可能复杂且耗时，但借助合适的工具和技术，可以使其变得简单。本文将探讨如何使用 Uber 的 H3 库进行索引和评分任务。
为了更好地理解，让我们定义一下假设的问题陈述。
作为零售公司的一名数据科学家，你被分配了一个重要的项目。公司需要组织和了解其运营的不同区域。
这个项目的最终目标是将区域划分为不同的部分，并识别每个区域中的关键热点。这些信息将帮助公司了解哪些区域对他们的业务最为重要，以及他们应集中精力的地方。
但你的工作并没有到此为止。公司还需要你为每个网点分配一个“重要性评分”。这个评分将基于之前提到的热点机制。那些位于潜在热点附近的网点将获得更高的重要性评分。
好的，让我们开始构建这个解决方案。
需要遵循的步骤
– 了解你的数据（Understand the data）
– 将区域拆分成小六边形
– 为区域内的每个六边形推导一个分数
– 定义出口重要性分数
1️⃣ 了解你的数据
根据问题陈述，我们有两个数据集。
-
兴趣点数据集 — 用于推导重要性分数
-
出口位置数据集
注意：提供的数据是随机生成的，不反映实际情况。你可以从 这里. 访问这些数据。
import pandas as pd
pois = pd.read_csv("pois.csv")
outlets = pd.read_csv("outlets.csv")
让我们绘制兴趣点类别分布。
此外，我们还需要该区域的地理信息。我已经通过geojson.io网站创建了区域的 JSON 文件。
{'type': 'FeatureCollection',
'features': [{'type': 'Feature',
'properties': {},
'geometry': {'coordinates': [[[80.19371259270753, 6.307037132469091],
[80.19371259270753, 6.072771744744614],
[80.4434199500783, 6.072771744744614],
[80.4434199500783, 6.307037132469091],
[80.19371259270753, 6.307037132469091]]],
'type': 'Polygon'}}]}
让我们把这些整合在一起。
目前我们所拥有的 | 作者提供的图片
在上图中，你可以看到两个图层。
-
矩形 geojson 图层 — 代表我们正在分析的区域
-
兴趣点图层 — 提供该区域的兴趣点数据集
2️⃣ 将区域拆分成小六边形
为了在特定区域内确定热点，有必要将区域划分为较小的部分。一个方法是使用六边形。六边形在这个背景下特别有用，因为它们具有一些优点，如能够均匀地划分区域和准确表示地理形状。如果你对为何使用六边形或它们提供的好处有进一步问题，可以参考 H3 的官方文档获取更多信息。
简而言之，通过选择正确的分辨率，可以通过一个六边形单元表示多个位置。
将多个点映射到一个六边形 | 作者提供的图片
[## H3 | H3
六边形层级地理空间索引系统 H3 提供了一个易于使用的 API，用于将坐标索引到六边形网格中……
h3geo.org.](https://h3geo.org/?source=post_page-----b77a7599b4c--------------------------------)
我们可以使用 H3 的*polyfill*函数将给定区域拆分成六边形。但在此之前，我们必须将区域 geojson 文件转换为 shapely 多边形格式。
from h3 import h3
from shapely.geometry import Polygon,mapping
def dict_to_shapely(d):
coords = d['features'][0]['geometry']['coordinates'][0]
return Polygon(coords)
def polygon_to_h3(polygon,resolution):
polygon = mapping(polygon)
hexas = h3.polyfill(polygon,res=resolution,geo_json_conformant=True)
return list(set(hexas))
#read geojson file
with open("sample-bbx.geojson",'r') as f:
polygon_dict = json.loads(f.read())
shapely_polygon = dict_to_shapely(polygon_dict)
f.close()
h3_idx = polygon_to_h3(shapely_polygon,resolution=8)
h3_idx = pd.DataFrame(pd.Series(h3_idx,name='h3'))
在上述代码片段中，polygon_to_h3函数将输入的多边形转换为六边形列表。你还需要提供分辨率参数。分辨率参数定义了六边形的大小。根据提供的分辨率参数，我们可以将给定区域拆分成六边形，如下所示。
将 geojson 瓦片转换为六边形列表 | 作者提供的图片
在这种情况下，使用了 8 的分辨率。根据文档，分辨率为 8 的六边形面积为 0.737327598 平方公里。在这个特定的情况下，区域被划分为 912 个六边形。
3️⃣ 为区域内每个六边形推导得分
由于我们已经将区域划分为六边形，我们现在可以继续为每个六边形推导得分。
该过程的方法包括以下内容：
-
为数据集中每个兴趣点分配唯一的六边形 ID。
-
建立评分标准，对每个类别给予加权。
-
将带分数的数据集与区域数据集合并。
scores = 1- (pois['category'].value_counts(normalize=True)) * 3
scores_dict = dict(zip(scores.index,scores.values * 100))
pois['weight'] = pois['category'].map(scores_dict)
grouped_df = pois.groupby("h3").agg(
counts=('category','count'),
score = ('weight','sum')
).sort_values(by='counts',ascending=False)
grouped_df.head()
让我们将评分层叠加在基础区域六边形上。右侧的图示展示了如何将得分分配给每个六边形。条形的高度反映了重要性得分。
图片来源于作者
让我们检查推导出的得分如何分布。
grouped_df.reset_index(inplace=True)
full_scored_idx = h3_idx.merge(grouped_df,on='h3',how='left').fillna(0)
full_scored_idx['score'].plot.hist(color='skyblue',edgecolor='black',bins=5,figsize=(10,8))
plt.title("Score Distibutin")
plt.xlabel("Score")
plt.ylabel("Frequency")
plt.grid(visible=True, linestyle='-.')
plt.show()
得分分布 | 图片来源于作者
现在区域内的所有六边形都已评分，是时候创建一个新的数据层进行可视化了。
区域的热点 | 图片来源于作者
请注意上面的颜色编码，灰蓝色的六边形表示重要性低的六边形，而粉色和米色的六边形表示该区域内相对重要的六边形。
4️⃣ 定义出口重要性分数
我们现在有一个包含所选区域所有六边形得分的数据集。这使我们能够识别该区域最关键（热点）六边形。这解决了问题的一半。接下来，我们需要为出口定义重要性得分。
将遵循的步骤包括以下内容：
-
为出口数据集分配唯一的六边形 ID。
-
对于每个出口，识别最近的 K 个六边形。
-
根据 K 个六边形的得分定义出口的重要性。
h3_to_score = dict(zip(full_scored_idx['h3'],full_scored_idx['score']))
get_krings = lambda h:h3.k_ring(h,k=3)
outlets_h3['rings'] = outlets_h3['h3'].apply(get_krings)
print(outlets_h3[['category','h3','rings']])
在上面的代码片段中，k_ring 函数用于检索每个出口最近的 K 个六边形。你还会注意到在这个例子中，K 被设置为 3。
让我们可视化数据，以更好地理解我们在这一步中所完成的工作。
图片来源于作者
现在我们差不多完成了。由于我们已经定义了该区域内每个六边形单元的得分，我们需要做的就是汇总每个出口附近的六边形得分。
outlets_exploded = outlets_h3.explode(column='rings')
outlets_exploded['weight'] = outlets_exploded['rings'].map(h3_to_score)
#can contains hexagons not relavent to current region
outlets_exploded.dropna(inplace=True)
#final scoring for outlets
outlets_exploded.groupby("category").agg(outlet_importance=("weight",'sum'))
最终评分 | 图片来源于作者
很好，我们现在已经成功建立了每个出口的重要性分数。这样，两个目标都已经实现了。做得好！ 😃
你可以在 GitHub 仓库下访问所有代码。
[## GitHub - Ransaka/Spatial-indexing-with-h3
目前无法执行该操作。你在另一个标签页或窗口中登录了。在其他标签页或…
github.com](https://github.com/Ransaka/Spatial-indexing-with-h3.git?source=post_page-----b77a7599b4c--------------------------------)
结论
这个简化的例子突出了 H3 库的一些宝贵功能。从测量两个地理位置之间的距离到进行最近邻分析，H3 库提供了广泛的选项。通过在适当的时候纳入空间索引，不仅可以节省资源，还能降低成本。原因在于，对点数据进行数学运算可能会消耗大量内存，但通过索引可以大大减少内存使用。
感谢阅读！可以在LinkedIn与我联系。
geotiff.js: 如何获取纬度-经度坐标的投影 GeoTIFF 数据
作者提供的图片
使用 Javascript 将纬度和经度坐标重新投影到 GeoTIFF 的坐标系统中，并使用 geotiff.js 查询数据
·
关注 发表在 Towards Data Science · 10 min 阅读 · 2023 年 3 月 12 日
–
Javascript 已被使用多年来提供交互式网页地图，这些地图通常由矢量数据和 RGB 图像切片组成。这些简单的前端 historically 需要更强大的语言和服务器技术来实际提供和渲染正在可视化的地理空间数据，以及查询或分析它们的机制。
随着 Javascript 语言和生态系统的成熟，曾经几乎不可能的事情现在变得简单而且性能不错。多亏了专注的开源开发者，现在在前端（浏览器）或后端（NodeJS）直接处理地理空间光栅数据变得相当容易。甚至还有多线程支持！
让我们看看在纯 Javascript 中曾经非常困难的一个功能：查询特定坐标的地理空间光栅数据。这些数据集通常以 GeoTIFF 格式提供，为了处理这些数据集，我们将使用一个更流行的处理文件格式的库：geotiff.js。
使用这个库从 GeoTIFF 中提取像素数据很简单。但是，给定一个纬度-经度坐标，你怎么知道读取哪个像素？如果 GeoTIFF 在一个投影坐标系统中怎么办？
与GDAL（一个流行的开源地理空间库，具有 NodeJS 绑定）不同，geotiff.js没有重投影功能：它仅仅是解析 GeoTIFF 文件。这意味着我们的解决方案将比GDAL中的等效方法稍微复杂一些，后者在最简单的情况下只需一个命令：
gdallocationinfo -valonly -wgs84 <geotiff filename> <longitude> <latitude>
我们只需编写一些代码来镜像 GDAL 在后台进行的额外工作：从光栅中提取投影信息，将纬度-经度坐标重新投影到光栅的坐标系统中，并确定该点适用于光栅的哪个像素。
为什么不首先重新投影整个 GeoTIFF？
通过将我们的 GeoTIFF 重新投影到像 EPSG:4326 这样的纬度-经度坐标系统，整个过程大部分可以被规避。但对于许多投影坐标系统，这会引入严重的不准确性。
例如，考虑查询 NBM 天气模型以获取科罗拉多州小熊峰顶的预测温度：
图片来源：作者
选择这个例子是因为它靠近数据集中一个像素的边缘。
现在让我们使用最近邻重采样将数据集重新投影到 EPSG:4326。相同的坐标现在返回了不同的值：
图片来源：作者
我们可以通过使用不同的重采样方法来稍微改进重投影，比如双线性插值：
图片来源：作者
但显然，我们丢失了一些原始天气模型的信息。我们可以通过使重新投影的光栅图像变得非常非常大来减少丢失的信息量，但即使在巨大的尺寸下，仍然会有一些边缘情况，其中靠近像素边界的地点无法正确提取原始天气模型数据。
这是否重要取决于数据的用途和分辨率。如果你使用的是高分辨率的数字高程模型来为长达数英里的路线构建高程剖面图，那么如果高程数据在某些地方略有偏差可能不是什么大问题。相反，如果你使用的是分辨率较粗的气象模型数据，那么处于错误的网格单元可能会导致该区域的预测结果有显著不同。
总体而言：最好使用原始投影来查询来自地理空间光栅的原始数据。
将纬度和经度重投影到光栅坐标系统
第一步是获取 GeoTIFF 的投影信息，并使用该信息将我们期望的坐标重新投影到光栅的坐标系统中。
使 GeoTIFF 成为地理空间光栅而非普通 TIFF 图像的一个重要因素是使用 GeoKeys（参见OGC GeoTIFF 规范）。我们可以使用 geotiff.js 读取给定光栅的 GeoKeys：
// assume the variable gtiff is the geotiff object
// that was created by loading data with geotiff.js
const image = await gtiff.getImage();
const geoKeys = image.getGeoKeys();
我们将把这些 GeoKeys 传递给一个投影库，它将处理繁重的计算，因为投影往往是相当复杂的！
例如，北美的气象模型数据往往使用 Lambert 等角圆锥投影。如果我们从一个源自椭球基准（通常是 WGS84）的纬度经度坐标开始，那么我们会面临一个相当复杂的公式：
Snyder, John (1987). “地图投影：工作手册（USGS 专业论文：1395）” — 公有领域
不需要重新发明轮子，因为开源社区已经投入了多年的工作，创建了全面的重投影库，这些库可以轻松处理各种复杂的地图投影之间的转换。
对于 JavaScript，最强大、最流行且功能齐全的库之一是Proj4js，它是常见的proj库的分支，支持GDAL。我们只需使用光栅的 GeoKeys 生成一个proj字符串，并将最终字符串传递给库即可。
然而，有很多有效的 GeoKeys：
/**
* Geokeys. If you're working with `geotiff` library, this is result of `image.getGeoKeys()`.
* @typedef {Object} module:geokeysToProj4.GeoKeys
* @property {number} GeographicTypeGeoKey See GeoTIFF docs for more information
* @property {number} GeogGeodeticDatumGeoKey See GeoTIFF docs for more information
* @property {number} GeogPrimeMeridianGeoKey See GeoTIFF docs for more information
* @property {number} GeogLinearUnitsGeoKey See GeoTIFF docs for more information
* @property {number} GeogLinearUnitSizeGeoKey See GeoTIFF docs for more information
* @property {number} GeogAngularUnitsGeoKey See GeoTIFF docs for more information
* @property {number} GeogAngularUnitSizeGeoKey See GeoTIFF docs for more information
* @property {number} GeogEllipsoidGeoKey See GeoTIFF docs for more information
* @property {number} GeogSemiMajorAxisGeoKey See GeoTIFF docs for more information
* @property {number} GeogSemiMinorAxisGeoKey See GeoTIFF docs for more information
* @property {number} GeogInvFlatteningGeoKey See GeoTIFF docs for more information
* @property {number} GeogPrimeMeridianLongGeoKey See GeoTIFF docs for more information
* @property {number} ProjectedCSTypeGeoKey See GeoTIFF docs for more information
* @property {number} ProjectionGeoKey See GeoTIFF docs for more information
* @property {number} ProjCoordTransGeoKey See GeoTIFF docs for more information
* @property {number} ProjLinearUnitsGeoKey See GeoTIFF docs for more information
* @property {number} ProjLinearUnitSizeGeoKey See GeoTIFF docs for more information
* @property {number} ProjStdParallel1GeoKey See GeoTIFF docs for more information
* @property {number} ProjStdParallel2GeoKey See GeoTIFF docs for more information
* @property {number} ProjNatOriginLongGeoKey See GeoTIFF docs for more information
* @property {number} ProjNatOriginLatGeoKey See GeoTIFF docs for more information
* @property {number} ProjFalseEastingGeoKey See GeoTIFF docs for more information
* @property {number} ProjFalseNorthingGeoKey See GeoTIFF docs for more information
* @property {number} ProjFalseOriginLongGeoKey See GeoTIFF docs for more information
* @property {number} ProjFalseOriginLatGeoKey See GeoTIFF docs for more information
* @property {number} ProjFalseOriginEastingGeoKey See GeoTIFF docs for more information
* @property {number} ProjFalseOriginNorthingGeoKey See GeoTIFF docs for more information
* @property {number} ProjCenterLongGeoKey See GeoTIFF docs for more information
* @property {number} ProjCenterLatGeoKey See GeoTIFF docs for more information
* @property {number} ProjCenterEastingGeoKey See GeoTIFF docs for more information
* @property {number} ProjCenterNorthingGeoKey See GeoTIFF docs for more information
* @property {number} ProjScaleAtNatOriginGeoKey See GeoTIFF docs for more information
* @property {number} ProjScaleAtCenterGeoKey See GeoTIFF docs for more information
* @property {number} ProjAzimuthAngleGeoKey See GeoTIFF docs for more information
* @property {number} ProjStraightVertPoleLongGeoKey See GeoTIFF docs for more information
* @property {number[]} GeogTOWGS84GeoKey Datum to WGS transformation parameters, unofficial key
*/
使用类似geotiff-geokeys-to-proj4的辅助库来处理geotiff.js的 geokeys 对象到proj字符串的正确转换可能是值得的。实际上，上面的注释块来自于那个库。
安装了这两个库后，我们无需编写太多代码就能设置我们的重投影功能：
import proj4 from 'proj4';
import * as geokeysToProj4 from 'geotiff-geokeys-to-proj4';
... // Not shown: importing geotiff.js and loading a geotiff file.
// for this example, the variable `gtiff` references the file.
const image = await gtiff.getImage();
const geoKeys = image.getGeoKeys();
const projObj = geokeysToProj4.toProj4( geoKeys );
const projection = proj4( `WGS84`, projObj.proj4 );
我们现在创建了一个名为“projection”的 proj 投影对象。使用起来非常简单：
const { x, y } = projection.forward( {
x: -105, // the longitude
y: 40 // the latitude
} );
上面的代码将纬度经度坐标重新投影到光栅坐标系统（可能是米或英尺）。
现在我们需要获取与这些新坐标关联的实际栅格像素。首先，让我们获取有关栅格大小及其边界框的一些信息：
const width = image.getWidth();
const height = image.getHeight();
const [ originX, originY ] = image.getOrigin();
const [ xSize, ySize ] = image.getResolution();
const uWidth = xSize * width;
const uHeight = ySize * height;
上述代码建立了栅格像素大小与投影坐标系统中栅格的实际坐标/大小之间的关系。
然后我们将投影坐标相对于栅格的原点坐标进行处理，查看它在栅格的宽度和高度中的位置。
// x and y come from the projection.forward example earlier
const percentX = ( x - originX ) / uWidth;
const percentY = ( y - originY ) / uHeight;
const pixelX = Math.floor( width * percentX );
const pixelY = Math.floor( height * percentY );
现在我们有了原始经纬度坐标的像素值！
你可能会想确保这个像素确实在栅格内。如果原始坐标没有与栅格相交，那么 pixelX 或 pixelY 可能会是负值或大于栅格的宽度/高度。
提取像素的 GeoTIFF 数据
如果像素在栅格内，那么考虑到一个单一波段的 GeoTIFF，我们有两种提取数据的方法。
第一种方法是将一个 1x1 像素的窗口传递给 geotiff.js 的“readRasters”函数：
const [ value ] = await image.readRasters( {
interleave: true,
window: [ pixelX, pixelY, pixelX + 1, pixelY + 1],
samples: [ 0 ]
} );
这种方法简单明了，但调用“readRasters”有一个小的开销，如果你查询数千个坐标，这个开销可能会累积。
在这种情况下，采用不同的方法可能会带来轻微的性能提升——但代价是 RAM，这对于非常大的 GeoTIFF 可能会相当可观。我们可以将整个栅格读取到内存中作为 ArrayBuffer：
const data = await image.readRasters( {
interleave: true,
samples: [ 0 ]
} );
然后我们只需要获取相应的索引：
const value = data[ width * pixelY + pixelX ];
实际上，geotiff.js 在内部缓存瓦片和条带，因此性能提升充其量是微乎其微的，并且对于适中的数据集通常不明显。
因此，如果查询的数据分布在 GeoTIFF 中的大多数条带或瓦片上，这种方法只会提供（轻微的）速度提升。否则，读取整个图像的开销，包括未使用的瓦片或条带，将超过潜在的性能提升，不过这可以通过读取一个大小调整到最小和最大 x 和 y 或坐标集合的窗口来缓解。
使用上述技术，可以在交互式网页地图中实现以下功能，无论是完全在浏览器中还是通过基本的 NodeJS API：
-
点击网页地图以获取特定点的栅格数据
-
使用底层栅格数据对点特征进行标注（例如，显示天气预报温度的城市名称的天气地图）
-
对栅格数据执行基本的几何操作（例如，使用点集合或边界框进行交集操作），并将该信息显示在表格中或导出为 CSV
-
为二维横截面或线生成图表或数据列表（例如，海拔剖面）
-
对原始 GeoTIFF 数据进行统计分析和数据转换
这种功能（更不用说实际将地理空间栅格数据渲染到地图上）传统上需要专用的 GIS 服务器软件和专用基础设施，交互式网络地图的作用仅仅是向服务器发送请求并显示响应。随着 JavaScript 生态系统的发展，浏览器变得越来越能够解析地理空间数据，从而减少了对专用后端技术的依赖。
传统上，使 GeoTIFF/地理空间栅格数据可供网络地图使用的一种流行方法是拥有一个提供通用开放地理空间联盟接口的 GIS 服务器应用程序——如 WMS（Web 地图服务）或 WMTS（Web 地图瓦片服务）。这些接口允许浏览器请求图像文件或数据的瓦片，这些图像在服务器端渲染后作为 RGB 图像文件发送到浏览器。
服务器还可以提供 WCS（Web 覆盖服务）接口，允许浏览器查询特定坐标的底层地理空间栅格数据。如果使用专有的 GIS 服务器技术（如 Esri 的 ArcGIS 产品），端点和语法可能会有所不同，但最终结果是相同的。
然而，随着云优化 GeoTIFF 标准的成熟以及使用 openlayers + geotiff.js 直接在浏览器中访问和渲染数据的能力，这种技术在某些场景下已经消除了对提供 WMS 和 WMTS 端点的 GIS 服务器的需求。本文描述的查询技术展示了一种去除对 WCS 或查询端点需求的方法。
取而代之的是，可以仅仅利用现有的云提供商（例如，一个 S3 存储桶）来托管 COGs，并让网络浏览器负责解析和处理地理空间数据。具有 NodeJS 后端的组织可以将 GIS 功能整合到这些堆栈中，从而减少对专用软件和额外硬件的需求。
最后，能够在浏览器中处理更多地理空间栅格数据的另一个巨大好处是减少或消除了对互联网连接的需求。对于智能手机地图应用程序，特别是那些旨在远程地区使用的应用程序，这为已经用 JavaScript 编写的应用程序提供了一种离线地理空间工具的途径。
这些应用程序中有几款非常流行，提供了诸如坡度角计算、地形剖面图、保存的天气地图等功能。这些应用允许你保存地图或区域以供离线使用——但一旦用户离线，许多功能会停止工作，因为它们依赖于服务器来实际解析和分析地理空间数据。
这些功能在失去互联网连接后不一定会中断——本文提到的技术只是众多可以让应用在离线时仍然保持功能的技术之一。
对 AUC 和 Harrell’s C 的直观理解
图形化方法
·
关注 发布于 Towards Data Science ·6 分钟阅读·2023 年 9 月 15 日
–
作者拍摄的照片
每个进入机器学习或预测建模领域的人都会遇到模型性能测试的概念。教科书通常仅在读者首先学习的内容上有所不同：回归及其均方误差（MSE），还是分类及其众多性能指标，如准确率、灵敏度或精确度等。尽管准确率可以简单地通过正确/错误预测的比率来计算，因此非常直观，但 ROC AUC 初看可能会让人感到困惑。然而，它也是评估预测器质量的常用参数。让我们首先解开其机制，以了解细节。
首先了解 AUC
假设我们建立了一个二分类器来预测样本属于某个类别的概率。我们的测试数据集具有已知类别，得到了以下结果，这些结果可以总结在混淆矩阵中，并在表格中详细报告，其中样本按预测为 P（正例）类别的概率排序。
混淆矩阵和详细预测表，显示了各个样本的概率。图片由作者提供。
ROC AUC 定义为 ROC（接收者操作特征）曲线下的面积。ROC 曲线是将真实正率（TPR）与假正率（FPR）绘制在一起的图 [Wikipedia]。TPR（也称为敏感性）是正确识别的正例数与所有正例数的比率。在我们的例子中，TPR 计算为 4/5（五个样本中有四个被正确分类为正例）。FPR 计算为错误分类为正例的负例数与实际负例总数的比率。在我们的例子中，FPR 计算为 2/6（如果我们将“阳性”阈值设置为 0.5，那么六个负样本中有两个被误分类为正例）。
我们可以根据 TPR 和 FPR 值绘制 ROC 曲线，并计算 AUC（曲线下面积）：
基于预测概率的 ROC 曲线。图片由作者提供。
个别 TPR/FPR 值是如何得到的？为此，我们考虑我们的概率表并计算每个样本的 TPR/FPR，将样本被视为正例的概率设置为表中给出的值。即使我们超过了通常的 0.5 阈值（样本通常被视为“负例”），我们仍继续将其视为正例。让我们在例子中遵循这一过程：
图片由作者提供
在阈值 0.81 时，五个正例中的一个被正确分类为正例，没有样本被预测为负例。我们继续，直到遇到第一个负例：
图片由作者提供
在这里，我们的 TPR 停滞在之前的值（五个正例中的三个被正确预测），但 FPR 增加，我们错误地将六个负样本中的一个分配给了正类。我们继续直到最后：
图片由作者提供
Voilà：我们得到了用于创建 ROC 曲线的完整表格。
为什么 Harrell 的 C 只是 AUC
那么 Harrell 的 C 指数（也称为一致性指数或 C 指数）呢？考虑预测特定疾病（例如癌症）发生后的死亡任务。最终，所有患者都会死，无论是否患有癌症——一个简单的二分类器帮助不大。生存模型考虑了结果（死亡）发生的持续时间。事件发生越早，个体遇到结果的风险越高。如果你要评估生存模型的质量，你会查看 C 指数（也称为一致性、也称为 Harrell 的 C）。
为了理解 C 指数的计算，我们需要引入两个新概念：允许对和一致对。允许对是指在观察期间有不同结果的样本对（例如：患者），即在实验进行时，其中一个患者经历了结果，而另一个则被删失（即尚未达到结果）。然后分析这些允许对，查看具有较高风险评分的个体是否经历了事件，而删失的个体没有。这些情况称为一致对。
简化一点，C 指数计算为一致对数量与允许对数量的比率（为了简化，我省略了风险相等的情况）。让我们通过一个例子来说明，假设我们使用了一个计算风险而不是概率的生存模型。以下表格仅包含允许对。如果具有较高风险评分的患者经历了事件（属于我们的“积极”组），则“Concordance”列设置为 1。id 只是之前表格中的行号。特别注意个体 4 与 5 或 7 的比较。
作者提供的图片
这使我们在 30 个允许对中得到 27 个一致对。比率（简化的 Harrell 的 C）是 C = 0.9，这让我们不禁想起之前计算的 AUC。
我们可以构建一个一致性矩阵，来可视化 C 统计量的计算过程，如 Carrington et al 所建议的那样。该图显示了实际正例的风险评分与实际负例的风险评分，并显示了所有对（绿色 + 红色）中正确排名对（绿色）的比例，如果我们将每个网格方块视为样本的表示：
Harrell 的 C 计算的一致性矩阵。作者提供的图片
一致性矩阵显示了底部右侧的正确排名对、顶部左侧的错误排名对，以及正好对应于我们之前看到的 ROC 曲线的边界。
在拆解构建 ROC 曲线和一致性矩阵的过程中，我们发现了一个相似点：在这两种情况下，我们都根据样本的概率/风险得分对其进行排序，并检查排序是否与实际情况相符。我们设置的分类概率阈值越高，假阳性越多。实际阳性病例的风险越低，实际阴性病例被误分类为阳性的可能性就越大。根据我们的排名数据绘制曲线，会得到形状和面积相同的曲线，这就是我们称之为 AUC 或 Harrell’s C 的曲线，具体取决于上下文。
我希望这个例子有助于培养对 AUC 和 Harrell’s C 的直观理解。
致谢
比较这两个参数的想法来源于一次富有成效的讨论，讨论发生在高级机器学习学习小组的聚会中，致敬 Torsten！
参考文献: Carrington, A.M., Fieguth, P.W., Qazi, H. 等人。用于评估机器学习算法中的不平衡数据的新一致部分 AUC 和部分 c 统计量。BMC Med Inform Decis Mak 20, 4 (2020)。https://doi.org/10.1186/s12911-019-1014-6
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
