逻辑回归的几何解释
理解如何使用几何解释导出逻辑回归的成本函数
在 Unsplash 上拍摄的 ThisisEngineering RAEng
逻辑回归是一种统计模型，它使用逻辑函数来模拟二元因变量。在几何解释术语中，逻辑回归试图找到一条线或一个平面来最好地区分这两类。逻辑回归处理几乎或完全线性可分的数据集。
术语线性可分是什么意思？
对于二元分类数据集，如果一条直线或一个平面可以几乎或完全分隔两个类，那么这样的数据集称为线性可分数据集。否则，如果这两个类不能被线或平面分开，则数据集不是线性可分的。
左图:线性分离，右图:非线性分离
从上面的二维样本数据集中，左边的样本数据集几乎可以用一条线线性分离，而对于右边的样本数据集，没有一条线可以将这两类点分开。
图片 2
在(图 2)中，该线最好地分隔了两类点，错误分类了 3 个点(红色圆圈)。
深入推导算法的几何解释:
图 3: Xi 和 Xj 是正确分类的点
对于上面的样本数据集，假设我们需要找到一个平面‘P’来分隔这两个类。
平面的一般方程由下式给出:
所以最后我们需要找到一个平面 P。设这两类点为 y_i = {1，-1}。对于任意两个随机点 x_i 和 x_j。
d_i =平面和 x_i 之间的距离
d_j =平面和 x_j 之间的距离
这里 W 垂直于平面。当一个点的方向垂直于平面时，那么距离为正，否则为负。
对于正确分类的点(图 3):
d_i > 0 和 d_j < 0
Now for a positive class point: **y_i = +1**
**(d_i * y_i) > 0**, since **d_i > 0** and **y_i > 0**For a negative class point: **y_i = -1**
**(d_i * y_i) > 0,** since **d_i < 0** and **y_i < 0**
所以对于正确分类的点， (y_i * d_i)总是正的
对于错误分类的点(图 4):
图 4: Xi 和 Xj 是错误分类的点
d_i < 0 和 d_j > 0
For a positive class point: **y_i = +1**
**(d_i * y_i) < 0**, since **y_i > 0 and d_i < 0 (because on wrong side)**For a negative class point: **y_i = -1**
**(d_i * y_i) < 0,** since **d_i > 0** and **y_i < 0 (because on wrong side)**
对于错误分类的点， (y_i * d_i)总是负的
因此，为了得到最优解，我们需要最大化 **(y_i * d_i)。**我们需要找到最优的 W，w0，它使下面的等式最大化。
乙状结肠挤压:
图 4
由于存在负类的离群点，因此上面获得的成本函数将获得平面“P2”作为最佳平面，但这不是真的。平面“P1”最好地分隔了这两类点。
异常值或极值点的存在会在很大程度上影响平面。为了避免这种情况，我们需要找到一个函数，使 (y_i * d_i) 值变小，如果它太大，id (y_i * d_i) 值变小，它应该保持小。
我们需要这样一个函数
图 5
如果 x 的值很大，那么 f(x)逐渐变小，如果 x 的值很小，那么它仍然很小。
乙状结肠功能:
sigmoid 函数的图形清楚地定义了它满足我们的条件。sigmoid 函数的数学方程:
来源:Google Plots，图片 6，sigmoid 函数的绘图
所以我们的方程式可以归结为:
最小化方程式:
如果 G(x)是单调递增函数，那么 G(F(x))也是单调递增函数。
把 F(x)作为上面推导的方程，G(x)作为 log(x)作为 log(x)是一个单调递增的函数。
来源:Google Plots，图 7，log_e(x)的绘图
方程式归结为:
正规化:
因为 log_e(z)的最小值是 0。因此，优化器将尝试将上述等式的值最小化为 0，这使得 log_e(z)的‘z’= 1。
所以 W，w0 趋向于无穷大来满足方程，这就趋向于进行逻辑回归。
添加 L2 正则化:
上面导出的等式是逻辑回归算法的成本函数。我们使用优化器来计算使上述成本函数最小化的最佳值 W，w0。*【lambda】*上式中的‘lambda’是超参数。
偏差方差权衡:
如果λ= 0，则上述方程将不包含任何会使模型过拟合的正则化项。
如果 lambda -> infinite(大值)，则正则化项的权重非常高，并且它将遮蔽该项的其余部分，这导致欠拟合模型。
预测查询点的目标类别:
图 7:查询
w，w0 通过最小化成本函数来计算
对于查询’ q1 '，y _ pred>0(W 方向的 q1)
对于查询’ q2 '， y_pred < 0 (与 W 方向相反的 q2)
感谢您的阅读！
计算机视觉中的几何变换:Python 示例的直观解释
旋转，平移和缩放转换解码。
佩顿·塔特尔在 Unsplash 上拍摄的照片
几何变换是任何图像处理流水线中最常见的变换操作之一。在今天的帖子中，我们将看看其中的三种变换: 旋转 ， 平移 和 缩放 ，然后仅使用 Numpy 从头开始构建它们。图 1 显示了我们想要实现的视觉效果。所以， 洛戈特！
图 1 我们的目标是将图像 a 顺时针旋转 45 度，以创建图像 b .来自带有 CC 许可证的 COCO 数据集的参考图像
OpenCV 方式
如果你使用任何图像库，如 OpenCV 或 PIL 内置函数，实现上述功能是非常简单的。使用 OpenCV，我们可以用两行代码来完成，如下所示。
使用 OpenCV 旋转图像
图 2 固定边界的旋转
我们得到的结果如图 2 所示。请注意，OpenCV 不会自动扩展图像的边界。为了确保我们看到整个旋转的图像，我们需要做另外两件事。首先，计算目标图像的大小，其次，由于新图像的中心不同于原始中心，我们需要考虑旋转矩阵中中心值的差异。有了这些补充，我们就是黄金！
使用动态边界旋转
这一次旋转包括了我们预期的整个图像！让我们来分析一下这些函数调用背后到底发生了什么。
旋转矩阵
我们使用上面的getrotationmatrix 2d()方法(片段 1 第 5 行)来创建一个旋转矩阵，我们稍后使用它来扭曲原始图像(片段 1 第 6 行)。这个函数以图像中心、旋转角度和缩放因子作为参数，并给出一个旋转矩阵。这个 旋转矩阵 到底是什么？嗯，它源于线性代数。让我们考虑一个任意的 2D 点[x，y]，那么旋转运算可以用下面的矩阵运算来表示。
2D 旋转变换
其中，R 是旋转矩阵
旋转矩阵 R
更简单地说，旋转矩阵给出了“函数 f”x’，y’= f(x，y) 将输入点映射到其旋转的对应点。矩阵 R 由两个列向量组成，表示我们的初始基向量在转换后的最终位置(如果你像我一样是一个 3blue1brown 的粉丝，这将立即敲响警钟！).
点旋转示例
让我们用单位圆把事情说得更清楚些。
在这个单位圆中，考虑头部在 B = [1，0]的两个向量 u 和头部在 D [0，1]的两个向量 w 。然后，这些向量围绕中心 A(=原点= [0，0])旋转某个角度θ、(=φ)，之后它们分别落在 C 点和 E 点。**
经过变换后，向量 v 可以用其分别以 F 和 H 为头的正交投影向量来表示。使用基本的三角学并记住向量 v 和 a 具有单位长度，可以容易地看出，以 F 为头的向量的长度是 b = cos(theta) ，以 H 为尾的向量的长度是 d= sin(theta) 。
导出向量 v 的正交分量
向量 v 因此可以表示为列向量*【cos(theta)sin(theta)】。同样，我们可以证明，向量 a 可以表示为列向量[-sin(theta)cos(theta)】。负号表示方向。*
**将这些列向量放在一起，我们得到了旋转矩阵 r。**为了简化矩阵乘法，通常在旋转矩阵中添加第三个轴。直觉上，这将是旋转 3D 结构的旋转轴。该轴上的所有点在变换后将保持不变，因此该附加轴对剩余的变换没有净影响。
旋转变换
类似地，平移操作可以由下面所示的平移矩阵来表示，其中 tx 和 ty 是 X 和 Y 方向上的平移量。
平移变换
还有我们的缩放矩阵，其中SxansSy是我们在 x 和 y 方向的缩放因子
缩放变换
注意:通常，旋转和平移被组合成如下所示的单个变换矩阵。这与绕原点旋转θ然后平移 tx 和 ty 的效果相同
旋转和平移操作打包到一个矩阵中
试试吧！取上面任意一个矩阵，给它们一些值乘以任意一个点【x，y】看看变换后的点值是不是你所期望的！
我们现在有了表达几何变换的所有必要成分。请记住几件重要的事情:
- 矩阵乘法是不可交换的也就是说，如果你将两个矩阵 A 和 B 相乘，A.B！因此，我们变换的顺序很重要。旋转一个点然后平移和先平移那个点然后旋转相同的因子是不一样的！
- 变换表达式的矩阵乘法是从右向左进行的。
- 注意你的坐标轴，因为方向很重要。像 OpenCV 这样的大部分图像框架都是考虑左上角的原点。这不是我们写方程的经典的“第一象限”，而是“第四象限”。实质上，y 轴的方向是反向的。我们可以将它直接放入旋转矩阵中，或者在计算过程中添加一个负号(我做第二个是为了保持转换操作与文献一致)。
- 旋转一般是围绕图像中心进行的。
很好，现在让我们回到我们的原始图像，看看我们需要做什么来得到想要的转换图像。操作的确切顺序如下:
- 平移图像，使图像的中心为原点。这是因为我们希望通过图像的中心而不是左上角来旋转图像，左上角通常是图像中的原点像素/坐标。让我们把这些翻译因子设为【tx1，tx2】。****
- 将图像旋转我们想要的角度θ。****
- 将图像平移回其原始中心。让我们把这个概括为【tx2，ty2】。****
- 计算新的中心，并利用这些新的和旧的中心之间的差异来平移图像。还要考虑新图像的大小。(我们已经为 OpenCV 计算完成了这一步)。设这些平移因子为【CX _ shift，cy_shift】。****
用矩阵的方式，我们可以将上述四个步骤表达如下。从右到左，我们将中心平移到原点(右数第一个矩阵)，绕中心旋转并平移回中心(中心矩阵)，最后调整新维度的中心(左数第一个矩阵)。
我们的转换矩阵
我们在这里没有使用缩放变换，但是如果你也想缩放你的图像，只需要把缩放变换加到上面的等式中(在正确的地方！).简化上述内容并替换a = cos(θ)**，b = sin(θ)**
简化变换矩阵
这是我们最终的旋转矩阵，我们将在下一节中使用 Numpy 来旋转我们的图像。
数字之路
调用 OpenCV 方法既快速又简单，但是一点也不好玩！因此，我们将把我们讨论的所有内容放入代码中，并使用 Numpy 旋转图像！
由于这个脚本有点太长，无法粘贴到这里，请点击下面的链接查看完整代码。
成像实用程序脚本。在 GitHub 上创建一个帐户，为 borarak/imutils 的发展做出贡献。
github.com](https://github.com/borarak/imutils/blob/master/geometric/rotate.py)
在这之后，我们可以使用自己的函数旋转和平移图像！
今天就到这里，希望你喜欢。一如既往的感谢阅读！**
GeoPandas:实用指南
入门
测绘 1965 年至 2016 年的地震
由 Andrew Buchanan 在 Unsplash 上拍摄的照片
GeoPandas 是一个 Python 库，旨在处理地理空间数据。这使得基于地理位置创建可视化变得相当容易。
在本帖中，我们将可视化 1965 年至 2016 年间发生的重大地震。数据集可在 Kaggle 上获得。
GeoPandas 有两种主要的数据结构，即地理数据框架和地理系列。它们可以被认为是熊猫系列和数据框的子类。
让我们从安装和导入 GeoPandas 以及我们将使用的其他库开始。
pip install geopandas
import geopandasimport pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
可以通过修改熊猫数据框来创建地理数据框。因此，我们首先将数据集读入熊猫数据帧。
eq = pd.read_csv("/content/earthquakes.csv")eq.shape
(23412, 21)
该数据集包含超过 23412 个事件，其中大部分是地震。我们将过滤数据帧，使其仅包含地震数据。
eq = eq[eq['Type'] == 'Earthquake']
我们的分析中还有一些冗余的列，因此我也将过滤掉这些列。
eq = eq[['Date', 'Time', 'Latitude', 'Longitude', 'Depth', 'Magnitude']]eq.head()
(图片由作者提供)
我们有一个数据框架，包含了两万多次地震的数据、位置、深度和震级。为了使用 GeoPandas，我们需要将此 Pandas 数据框架转换为地理数据框架。
我们将按如下方式使用地理数据框架函数:
gdf = geopandas.GeoDataFrame(eq, geometry=geopandas.points_from_xy(eq.Longitude, eq.Latitude))gdf.head()
(图片由作者提供)
地理数据框架和 pandas 数据框架之间的区别是一个称为“几何”的地理系列。当空间方法应用于地理数据框架时，它将作用于几何列。
可以将“geometry”列看作是纬度和经度值的重新格式化版本。
我们现在将地震数据存储在地理数据框架中。下一步是绘制世界地图，使用“世界”地理数据框架可以轻松完成。
world = geopandas\
.read_file(geopandas.datasets.get_path('naturalearth_lowres'))world.columns
Index(['pop_est', 'continent', 'name', 'iso_a3', 'gdp_md_est', 'geometry'], dtype='object')
它包含关于国家及其位置的基本信息。我们现在画一张空的世界地图。
world.plot(color='white', edgecolor='black', figsize=(12,8))
(图片由作者提供)
为了绘制地震图，我们将创建一个世界地图的轴对象，然后根据“几何”列绘制地震。
ax = world.plot(color='white', edgecolor='black', figsize=(16,12))
gdf.plot(ax=ax, color='red', markersize=2)
plt.show()
(图片由作者提供)
这张地图包含了 1965 年至 2016 年间发生的所有重大地震。如果你在谷歌上快速搜索地震断层线，你会发现它们与上面的地图重叠。
markersize 参数调整定位地震的标记的大小。您还可以传递一个列名，标记的大小将根据该列中的值进行调整。我想过用大小来调整标记的大小，但是差别似乎并不明显。
我们还可以画出特定地点的地震地图。例如，日本发生过多次地震。我不确定，但它可能是世界上地震最多的国家。
关注特定国家的一种方法是根据纬度和经度值过滤地震。日本的纬度和经度值如下所示:
- 纬度= 36.204824
- 经度= 138.252924
我们可以围绕这些值创建一个范围，用作过滤范围。
japan_lat = 36.204824
japan_long = 138.252924japan_eq = eq[(eq.Latitude > 30) & (eq.Latitude < 42) & (eq.Longitude > 130) & (eq.Longitude < 145)]japan_eq = japan_eq.reset_index(drop=True)
我调整了范围，使位置占据了日本的面积。请注意，这些数值并不是日本的国界。
让我们创建一个仅包含日本境内或周边发生的地震的地理数据框架。
japan_gdf = geopandas.GeoDataFrame(japan_eq, geometry=geopandas.points_from_xy(japan_eq.Longitude, japan_eq.Latitude))
我们将绘制日本地图，并在日本 gdf 中标记地震。
ax = world[world.name == 'Japan'].plot(color='white', edgecolor='black', figsize=(12,8))japan_gdf.plot(ax=ax, color='blue', markersize=japan_gdf['Magnitude']*4)plt.title("Earthquakes, 1965-2016", fontsize=16)plt.show()
为了过滤地图，我们使用了世界地理数据框的“名称”列。
这是由此产生的日本地震地图。
(图片由作者提供)
日本发生过多次地震。标记的密度表明大部分都在东海岸。
结论
GeoPandas 是一个函数库，可以加速地理空间可视化的创建过程。它提供了许多功能和方法来丰富地图。
如果你正在或者计划使用地理空间数据，我强烈推荐你查看 GeoPandas 的文档。它们还提供了一些示例，有助于更容易地调整函数和方法。
感谢您的阅读。如果您有任何反馈，请告诉我。
Geopandas 安装 Windows 的简便方法！
如果你想做一些房地产分析， GeoPandas 软件包提供了一种处理地理信息的惊人方法。它扩展了 pandas 使用的数据类型，允许对几何类型进行空间操作。
空间操作的例子有地图叠加(将两张地图结合在一起)、简单缓冲，但更普遍的是，GeoPandas 可用于地理可视化。
安装
虽然 GeoPandas 是一个强大的空间操作包，但安装对一些人来说有点困难。它需要难以协调的依赖性。
我发现了一种在我的笔记本电脑上有效安装它的方法，下面的步骤可能对你也有用。
这些步骤假设你已经安装了 wheel ( *pip install )。cmd 上的 whl)。
- 转到Python 扩展包的非官方 Windows 二进制文件。
- 在一个特定的文件夹中下载以下二进制文件:GDAL、Pyproj、Fiona、Shapely 和 Geopandas 匹配 Python 的版本，以及你的笔记本电脑上安装的是 32 位还是 64 位操作系统。
(例如对于 Python 3.7x(64 位)，GDAL 包应该是GDAL-3 . 1 . 2-cp37-cp37m-win _ amd64 . whl。)
下载的依赖于一个文件夹
3.转到下载二进制文件的文件夹。
使用命令提示符并转到下载二进制文件的文件夹
4.重要提示:使用 pip 安装的以下安装顺序是必要的。注意文件名。如果文件名是正确的，它应该工作:(提示:键入“pip install”后跟一个空格，并键入二进制文件的前两个字母，然后按 Tab 键。(例如 pip install gd(按 Tab))
- pip 安装。\ GDAL-3 . 1 . 1-cp37-cp37m-win _ amd64 . whl
- pip 安装。\ pyproj-2 . 6 . 1 . post 1-cp37-cp37m-win _ amd64 . whl
- pip 安装。\ Fiona-1 . 8 . 13-cp37-cp37m-win _ amd64 . whl
- pip 安装。\ Shapely-1 . 7 . 0-cp37-cp37m-win _ amd64 . whl
- pip 安装。\ geo pandas-0 . 8 . 0-py3-无-任何
5.就是这样！检查 GeoPandas 是否已正确安装，并通过导入它和查看模块上的帮助选项进行浏览。
注意:如果以下步骤不起作用，可能是因为系统上安装的其他不同的软件包与当前版本不兼容。我重新安装了 Anaconda，并确保首先安装 GeoPandas。
使用 Python 和自然语言处理进行地理解析
提取地名和指定坐标-教程
纳蕾塔·马丁在 Unsplash 上的照片
大量的可用文本数据具有可以从自动信息提取中受益的位置特征。自然语言处理(NLP)在过去的五年中取得了显著的进步。然而，从文本中提取地理信息的研究仍处于起步阶段。在本教程中，我们使用 Python 和 NLP 来 Geoparse twitter 数据集。
地质公园
Geoparsing 是一个地名解析过程，将对地点的自由文本描述(如“伦敦以东两公里”)转换为地理标识符(带经纬度的坐标)。
地理解析在许多应用中是必不可少的，包括地理信息检索(GIR)、地理信息提取(GIE)和地理信息分析(GIA)任务。
我们可以使用地理解析来确定文档的地理范围，解码灾难响应、商业新闻分析以及其他多个领域的位置信息。
为了说明什么是地质公园，让我们考虑这个讽刺标题示例。
“抗议者偷走纽约市的环卫车，用它们来阻挡特朗普大厦”
通常，地理区划包含两个部分:地名识别和地名解析。首先，是地名提取或识别[纽约市，特朗普大厦]。下一步是将地名链接到地理坐标[(40.768121，-73.981895)，(40.762347，-73.973848)]。
在下一节中，我们将使用 Python 地理解析库 Mordecai 对简单文本进行地理解析。
Python 地理解析示例
对于本教程，我们将使用 Mordecai 库进行地理解析。Mordecai 是全文地理解析 Python 库。使用这个库，您可以从一段文本中提取地名，将它们解析到正确的位置，并返回它们的坐标和结构化的地理信息。
让我们从一个简单的地质分析例子开始。Mordecai Python Geoparsing 库具有 Geoparse 函数，该函数接收文本并从文本中返回结构化的地理信息。
from mordecai import Geoparsergeo = Geoparser()
geo.geoparse(“Eiffel Tower is located in Paris”)
对于任何文本输入，末底改返回文本中存在的位置特征。在这个例子中，它正确地预测了埃菲尔铁塔和巴黎城。有趣的是，与这两个位置相关联的纬度和经度是不同的。预测的埃菲尔铁塔坐标比巴黎这座城市还要具体。
[{'word': 'Eiffel Tower',
'spans': [{'start': 0, 'end': 12}],
'country_predicted': 'FRA',
'country_conf': 0.611725,
'geo': {'admin1': 'Île-de-France',
'lat': '48.85832',
'lon': '2.29452',
'country_code3': 'FRA',
'geonameid': '6254976',
'place_name': 'Tour Eiffel',
'feature_class': 'S',
'feature_code': 'MNMT'}},
{'word': 'Paris',
'spans': [{'start': 27, 'end': 32}],
'country_predicted': 'FRA',
'country_conf': 0.9881995,
'geo': {'admin1': 'Île-de-France',
'lat': '48.85339',
'lon': '2.34864',
'country_code3': 'FRA',
'geonameid': '2988506',
'place_name': 'Paris',
'feature_class': 'A',
'feature_code': 'ADM3'}}]
Mordecai Python 库采取了不同的步骤来实现这一结果。首先，它使用 spaCy 的命名实体识别从文本中提取地名。然后，它使用 Geonames gazetteer 查找地名的潜在坐标。最终过程使用神经网络从地名录条目中预测国家和地名。
地质公园推文
为了抓取 tweets，首先，我们设置了 Tweepy API 来抓取标签。下面这段代码使用 Tweepy 抓取一个标签(#BlackLivesMatter)，并将标签中的所有 tweets 保存到本地 CSV 文件中。
用 Tweepy 抓取推文
让我们读一下关于熊猫的 CSV 推文，看看前几个专栏。
df = pd.read_csv(“tweets.csv”, header=None, names=[“date”, “tweet”])
df.head()
推特数据帧
数据框现在保存了推文和文本的日期。让我们使用 Mordecai Geoparsing 功能来提取位置信息并分配坐标。我们设置了这个函数，它获取一个数据帧，并生成一个干净的数据帧，其中包含来自地理解析的附加位置信息。
我们的干净数据集现在已经提取了地名，并通过预测和预测的置信度为 tweet 文本中的每个地名分配了坐标。
清理 Geoparsed 数据帧
为了绘制使用 Mordecai Geoparsing Python 库提取的#BlackLivesMatter hashtag 的地理范围，我们现在可以使用任何您喜欢的地理空间数据可视化 Python 库。我用 Plotly Express 来绘制数据。
fig = px.scatter_mapbox(df_clean, lat=”lat”, lon=”lon”, size_max=15, zoom=1, width=1000, height=800, mapbox_style=”dark”)
fig.data[0].marker = dict(size = 5, color=”red”)
fig
#BlackLivesMatter 标签的地理范围-地理公园
结论
地理句法分析是从文本中自动提取位置特征的重要组成部分。在本教程中，我们了解了如何使用 Mordecai Geoparsing Python 库对文本进行 Geoparse 处理。
要使用 Mordecai 运行 Geoparsing，您需要安装它。您还需要有一个正在运行的 docker 容器。你可以在这里找到安装说明。
本教程的代码可以在这个 Jupyter 笔记本中找到。
permalink dissolve GitHub 是超过 5000 万开发人员的家园，他们一起工作来托管和审查代码，管理…
github.com](https://github.com/shakasom/geoparsing/blob/master/Geoparsing%20with%20Python.ipynb)
佐治亚理工的 MS Analytics 项目值得吗？
我在佐治亚理工学院 MS 分析项目的经历
2017 年末，我开始关注不同的在线硕士项目，重点关注数量。这包括统计学、数据科学、分析学、经济学和各种 MBA 课程。不幸的是，许多顶级项目，如卡内基梅隆大学的商业分析或伯克利的数据科学项目，都有 50，000 美元或更多的疯狂价格。虽然高薪不是我想要硕士学位的唯一原因，但高昂的学费大大降低了潜在的投资回报。在研究了几个星期的程序后，我终于找到了我要找的东西。这个分析程序的价格低于 15，000 美元，在不同的专业领域具有灵活性，并且来自一所著名的大学。这是佐治亚理工学院最近成立的在线分析科学硕士项目 (OMSA)。
快进到今天，我现在正在参加他们的金融建模课程。完成后，这将是我完成这个项目所需的 36 学分中的第 21 学分。如果一切顺利，我将在明年年底前完成这个项目。因为我的全职工作，我采取了保守的方法，每学期只报一门课。如果你对注重数量的在线硕士项目感兴趣，请继续阅读我对这个项目的诚实评论。
承认
选择乔治亚理工学院项目的便利之一是最低入学要求。他们不需要 GRE 或 GMAT 成绩，这在申请过程中节省了大量时间。缺乏需求既可以被视为优势，也可以被视为劣势。我对这个项目最大的抱怨之一是入学要求和课堂期望之间的巨大脱节。的要求声明你只有需要至少一门大学水平的课程或者“Python 中的计算机编程，达到 Python 中的计算入门的水平”的同等知识。我的计算机编程背景包括几门本科课程、在大会上为期两周的训练营以及一些工作经验；然而，我在课程中的两门必修课上挣扎了很久。
上层社会
数据分析计算简介
压力。愤怒。挫败感。正如我当时的任何一个室友可以告诉你的那样，在这个项目的第一堂课上，我一直处于混乱状态。数据分析计算入门让我想起了我的大学经历。这类似于我们在马里兰大学描述 STEM 专业的微积分 I 和 II。清除课程。根据佐治亚理工学院的分布报告，自 2018 年以来，约有 17%的学生最终退出了这门课程。
本课程的讲座只是浅尝辄止，所以大部分 python 编程都是我们自己学的。幸运的是，在在线可以找到大量的 python 教程和指南。虽然每周的家庭作业还不错，但考试却一点也不。我们进行了多次限时两天的考试。这听起来不错，只是一般人需要 10-12 个小时来完成每项考试。为什么需要 10 个小时才能充分检验一个人是否掌握 python？这在现实世界中是完全不现实的。我从这门课上学到了很多，但考试让我尝到了苦涩的滋味。
数据和可视化分析
在一个学期中，专注于软件工具的最佳数量是多少？如果你回答了 10 分或更高，那么这就是为你准备的课程！这个高级核心需求结合了不同的可视化和大数据工具的大杂烩，将它们投入到四个广泛的项目中，然后就到此为止。对于非分析人员或高管来说，这可能看起来合乎逻辑；但是技术岗位需要相当的深度才能成为专家。从概念的角度来看，这是真的，对于软件工具也是如此，因为概念可以很好地转换，BI 工具通常是相似的。如果你只花一周时间学习 Tableau，我保证几个月后你会忘记几乎所有的事情。然而，如果你花一整个学期的时间在一个老式的 MySQL 数据库上，这些概念和技能中的许多将无缝地转化为云中的现代数据库环境。我和我的一个同事一起亲眼目睹了这一点，他因为拥有 10 多年的 SQL 数据库经验，已经非常快地掌握了 Azure 中的数据工程。这个课程需要大量的努力，平均每周 14 个小时，根据 OMS 中央课程评估网站。这意味着大约 50%的学生每周需要超过 14 个小时。对于有家庭和/或全职工作的人来说，这是一项艰巨的工作。我会争辩说，这是一个不公平的的时间来期待 3 学分，但我跑题了。不幸的是，这门课的糟糕设计使得学生付出的努力与学到的知识不均衡。
模拟
模拟是我一生中上过的最好的课程之一。戈德曼教授对所有模拟事物的极大热情可以让任何讨厌数学的人变成书呆子。到目前为止，他是唯一一个在这个节目的预先录制的视频讲座中加入一些角色的教授。不断有老笑话，youtube 音乐视频的链接，闪烁的颜色或噪音。任何学生都知道，在学习方面，热情和创造力可以创造一个不同的世界。
虽然这门课程比大多数课程更重数学，但它的主题对许多类型的工作都非常有用。从交通模拟到医院的病人流动，每个行业都有模拟实验。这门课涵盖了概率分布、随机变量生成、泊松过程和在 ARENA 中设计模拟实验等主题。如果你对这些主题更详细的内容感兴趣，我把课堂讲稿整合到这里。
社区
有句老话说，你知道什么并不重要，重要的是你知道谁是 T2。虽然面对面交流或通过在线项目建立长期关系很难，但我认为这种经历最有价值的部分是在线社区。OMSA 项目管理着一个松散的网络，包括每个班级的特定频道、学生发布的工作信息等等。每时每刻你都会收到来自世界各地成千上万聪明的学生和助教的有趣帖子。我在我们的 Slack 频道上发布了一个调查，看看在开始这个项目之前，每个人最感兴趣的话题是什么。有趣的是，这个游泳池非常多样化。大约 60%的学生在数学、工程、商业或经济方面有很强的背景。知识和经验的多样性促进了难以置信的对话，以及在数据分析中对思想、链接和笔记的奇妙的、经过过滤的选择。它比任何谷歌搜索都好！
OMSA 松弛频道民意调查
虽然社区是在线的，但还是有一些机会可以见到人。对于拥有大量学生追随者的城市，比如纽约，有专门的渠道，在那里人们会安排面对面的会面。我去年在纽约参加了一个，遇到了几个同学。
替代方案和成本
随着数据在过去几年中变得越来越受欢迎，在线硕士选项比我在 2017 年时要多得多。然而，一旦你把学费考虑进去，选择过程就变得容易多了。如果你需要一个 30，000 美元以下的知名项目，那么剩下的选择就不多了。这里有几个符合这个标准的项目，按照美国新闻研究生计算机科学或统计学排名的顺序排列。
选项#1*
MS Analytics，佐治亚理工
学费:～1 万美元(未收费)
CS 排名:#8
选项#2*
德克萨斯大学奥斯汀分校 MS 数据科学
学费:～10000 美元(未收费)
CS 排名:#10
选项#3
宾夕法尼亚州立大学世界校区 MS 应用统计学 学费:～27000 美元(未收费)
统计排名:#20 (CS 排名:#30)
选项#4
印第安纳大学 MS 数据科学
学费:~ $ 23400(未收费)
CS 排名:#55
选项#5
MS Data Analytics，科罗拉多州立大学
学费:~ 19000(未收费)
CS 排名:#75
*GT MS Analytics 与德克萨斯大学 MS 数据科学项目的比较
在线硕士课程的两种选择是新兵训练营和自学。我在大会上有过很好的经历，但是训练营的变化很大。证书主义的重要性正在改变，但如果你打算参加训练营，我会做你的尽职调查。至于自学，我在节目中学到的东西在网上都能很容易找到。你在 Coursera 和麻省理工学院有大量的在线公开课。我最喜欢的学习数据科学、分析和数学的资源可以在这里找到。自学需要大量的训练和时间来设计课程。就我个人而言，我知道我必须参加一个正式的项目来推动自己学习这些材料。
回首
当我回顾我参加 OMSA 项目的决定时，我不禁要分析我是否做出了正确的决定。在你做出选择之前，你只能在网上找到这么多信息。以我现在所知，我还会选择这个项目吗？老实说，我不知道，但是这山望着那山高。
更新:第二部分的评论在这里发表。
Apache Sedona 在大规模处理地理空间数据方面表现突出
空间数据洪流
在过去十年中，可用的地理空间数据量大幅增加。这些数据包括但不限于:天气地图、社会经济数据和地理标记的社交媒体。例如，美国国家航空航天局的宇宙飞船一直在监测地球的状况，包括陆地温度、大气湿度。截至今天，美国宇航局已经发布了超过 22PB 的卫星数据。今天，我们在全球拥有近 50 亿台移动设备。因此，移动应用程序产生了大量的网络视频数据。例如，Lyft、优步和摩拜单车每天从数百万骑行者那里收集万亿字节的 GPS 数据。事实上，我们在移动设备上所做的一切都会在地球表面留下数字痕迹。此外，配备 GPS 的移动设备和物联网(IoT)传感器的前所未有的普及，导致不断产生结合周围环境状态的大规模位置信息。例如，一些城市已经开始在十字路口安装传感器，以监测环境、交通和空气质量。
理解隐藏在数据中的丰富的地理空间属性可能会极大地改变我们的社会。这包括许多正在进行深入研究的课题，如气候变化分析、森林砍伐研究、人口迁移、疫情扩散分析、城市规划、交通、商业和广告。这些数据密集型地理空间分析应用高度依赖底层数据管理系统(DBMSs)来高效地检索、处理、争论和管理数据。
阿帕奇塞多纳(前 GeoSpark)概述
Apache Sedona(以前的 GeoSpark)(http://sedona.apache.org)是一个集群计算框架，可以大规模处理地理空间数据。GeoSpark 扩展了 Apache Spark 中的核心数据结构弹性分布式数据集(RDD ),以在集群中容纳大地理空间数据。SpatialRDD 由分布在 Spark 集群中的数据分区组成。空间 RDD 可以通过 RDD 变换创建，也可以从永久存储的文件中加载。该层提供了许多 API，允许用户从各种数据格式中读取异构空间对象。
GeoSpark 允许用户使用现成的空间 SQL API 和 RDD API 进行查询。RDD API 提供了一组用操作编程语言编写的接口，包括 Scala、Java、Python 和 r。空间 SQL 接口为用户提供了一个声明性语言接口，因此他们在创建自己的应用程序时可以享受更多的灵活性。这些 SQL API 实现了 SQL/MM Part 3 标准，该标准广泛用于许多现有的空间数据库，如 PostGIS(在 PostgreSQL 之上)。接下来，我们将展示如何使用 GeoSpark。
Apache Sedona 中支持的空间数据源
过去，研究人员和从业人员开发了许多用于不同目的的地理空间数据格式。然而，异构数据源使得将地理空间数据集成在一起极其困难。例如，WKT 格式是一种广泛使用的空间数据格式，它将数据存储在人类可读的制表符分隔值文件中。Shapefile 是一个空间数据库文件，它包括多个子文件，如索引文件和非空间属性文件。此外，地理空间数据通常具有不同的形状，如点、多边形和轨迹。
目前，Sedona (GeoSpark)可以从本地磁盘、亚马逊 S3 和 Hadoop 分布式文件系统(HDFS)等不同的外部存储系统中读取 WKT、WKB、GeoJSON、Shapefile 和 NetCDF / HDF 格式的数据到空间 RDDs。空间 RDDs 现在可以容纳七种类型的空间数据，包括点、多点、面、多面、线串、多线串、几何集合和圆。此外，具有不同形状的空间对象可以共存于同一空间 RDD 中，因为 Sedona 采用了灵活的设计，该设计概括了不同空间对象的几何计算接口。
**空间 RDD 内置几何库:**空间数据科学家需要在 Apache Sedona 中挖掘空间对象的一些几何属性，比如周长、面积、交集，这是相当常见的。空间 RDD 配备了一个内置的几何库来执行大规模的几何运算，因此用户不会涉及复杂的计算几何问题。目前，系统在该库中提供了 20 多种不同的功能，并将它们分为两个不同的类别
常规几何函数应用于空间 RDD 中的每个空间对象。对于每个对象，它都会生成相应的结果，如周长或面积。输出必须是常规 RDD 或空间 RDD。
几何聚合函数应用于空间 RDD 以生成聚合值。它只为整个空间 RDD 生成单个值或空间对象。例如，系统可以计算整个空间 RDD 的边界框或多边形并集。
使用 RDD API 运行查询
这里，我们概述了使用 GeoSpark RDD API 创建空间 RDD 和运行空间查询的步骤。示例代码是用 Scala 编写的，但也适用于 Java。
**设置依赖关系:**在开始使用 Apache Sedona(即 GeoSpark)之前，用户必须将相应的包作为依赖关系添加到他们的项目中。为了便于管理依赖项，GeoSpark 的二进制包托管在 Maven 中央存储库中，该存储库包含了全世界所有基于 JVM 的包。只要项目由 Apache Maven 和 sbt 等流行的项目管理工具管理，用户就可以通过在 POM.xml 和 build.sbt 等项目规范文件中添加工件 id 来轻松添加 Apache Sedona。
**初始化 Spark Context:**Spark 或 Apache Sedona 中的任何 RDD 都必须由 Spark Context 创建。因此，GeoSpark 应用程序的第一项任务是启动 SparkContext。下面的代码片段给出了一个例子。为了使用自定义空间对象和索引序列化程序，用户必须在 SparkContext 中启用它们。
val conf = new SparkConf()conf.setAppName(“GeoSparkExample”)// Enable GeoSpark custom Kryo serializer
conf.set(“spark.serializer”, classOf[KryoSerializer].getName)conf.set(“spark.kryo.registrator”, classOf[GeoSparkKryoRegistrator].getName)val sc = new SparkContext(conf)
**创建一个空间 RDD:**Spatial rdd 中的空间对象并不局限于某种几何类型，而是对更多的场景开放。它允许包含混合类型几何图形的输入数据文件。例如，WKT 文件可能包括三种类型的空间对象，如线串、多边形和多重多边形。目前，该系统可以加载许多不同数据格式的数据。这是由一组文件阅读器完成的，如 WktReader 和 GeoJsonReader。例如，用户可以调用 ShapefileReader 来读取 ESRI 形状文件。
val spatialRDD = ShapefileReader.readToGeometryRDD(sc, filePath)
变换坐标参考系统: Apache Sedona 不控制空间 RDD 中对象的坐标单位(即基于度数或基于米)。计算两个坐标之间的距离时，GeoSpark 只需计算欧几里得距离。在实践中，如果用户想要获得精确的地理空间距离，他们需要将坐标从基于度数的坐标参考系统(CRS)，即 WGS84，转换到平面坐标参考系统(即 EPSG: 3857)。GeoSpark 为用户提供了这一功能，使他们可以对空间 RDD 中的每个对象执行这一转换，并使用集群扩展工作负载。
// epsg:4326: is WGS84, the most common degree-based CRSval sourceCrsCode = “epsg:4326" // epsg:3857: The most common meter-based CRSval targetCrsCode = “epsg:3857"objectRDD.CRSTransform(sourceCrsCode, targetCrsCode)
**构建空间索引:**用户可以调用 API 在空间 RDD 上构建分布式空间索引。目前，系统提供两种类型的空间索引，四叉树和 R 树，作为每个分区上的本地索引。这一步的代码如下:
spatialRDD.buildIndex(IndexType.QUADTREE, false) // Set to true only if the index will be used join query
**编写一个空间范围查询:**空间范围查询返回位于一个地理区域内的所有空间对象。例如，范围查询可能会找到菲尼克斯大都市地区的所有公园，或者返回用户当前位置一英里范围内的所有餐馆。就格式而言，空间范围查询以一组空间对象和一个多边形查询窗口作为输入，并返回位于查询区域内的所有空间对象。空间范围查询将范围查询窗口和空间 RDD 作为输入，并返回与查询窗口相交/被查询窗口完全覆盖的所有几何。假设用户拥有空间 RDD。他或她可以使用以下代码对该空间 RDD 发出空间范围查询。空间范围查询的输出格式是另一种空间 RDD。
val rangeQueryWindow = new Envelope(-90.01, -80.01, 30.01, 40.01) /*If true, return gemeotries intersect or are fully covered by the window; If false, only return the latter. */val considerIntersect = false// If true, it will leverage the distributed spatial index to speed up the query executionval usingIndex = falsevar queryResult = RangeQuery.SpatialRangeQuery(spatialRDD, rangeQueryWindow, considerIntersect, usingIndex)
**编写空间 K 最近邻查询:**将 K、查询点和空间 RDD 作为输入，并在 RDD 中找到与查询点最近的 K 个几何。如果用户拥有空间 RDD，他或她可以执行如下查询。空间 KNN 查询的输出格式是包含 K 个空间对象的列表。
val geometryFactory = new GeometryFactory()val pointObject = geometryFactory.createPoint(new Coordinate(-84.01, 34.01)) // query pointval K = 1000 // K Nearest Neighborsval usingIndex = falseval result = KNNQuery.SpatialKnnQuery(objectRDD, pointObject, K, usingIndex)
**编写一个空间连接查询:**空间连接查询是用一个空间谓词组合两个或更多数据集的查询，比如距离和包含关系。生活中也有一些真实的场景:告诉我所有有湖的公园，告诉我所有 500 英尺内有杂货店的加油站。空间连接查询需要两组空间对象作为输入。它从这两个数据集的叉积中找到一个子集，使得每个记录都满足给定的空间谓词。在 Sedona 中，空间连接查询将两个空间 RDDs A 和 B 作为输入。对于 A 中的每个对象，从 B 中查找被它覆盖/相交的对象。a 和 B 可以是任何几何类型，并且不必具有相同的几何类型。空间 RDD 空间分区可以显著提高连接查询的速度。有三种空间划分方法:KDB 树、四叉树和 R 树。两个空间 rdd 必须由同一个空间分区网格文件进行分区。换句话说，如果用户首先对空间 RDD A 进行分区，那么他或她必须使用 A 的数据分割器对 b 进行分区。示例代码如下:
// Perform the spatial partitioningobjectRDD.spatialPartitioning(joinQueryPartitioningType)queryWindowRDD.spatialPartitioning(objectRDD.getPartitioner)// Build the spatial indexval usingIndex = truequeryWindowRDD.buildIndex(IndexType.QUADTREE, true) // Set to true only if the index will be used join queryval result = JoinQuery.SpatialJoinQueryFlat(objectRDD, queryWindowRDD, usingIndex, considerBoundaryIntersection)
使用 SQL APIs 运行空间查询
这里，我们概述了使用 GeoSpark 的空间 SQL 接口管理空间数据的步骤。SQL 接口遵循 SQL/MM Part3 空间 SQL 标准。具体来说，GeoSpark 将可用的空间 SQL 函数分为三类:(1)构造函数:创建一个几何类型列(2)谓词:评估一个空间条件是真还是假。谓词通常用于 WHERE 子句、HAVING 子句等(3)几何函数:对给定的输入执行特定的几何运算。这些函数可以生成几何图形或数值，如面积或周长。
为了使用该系统，用户需要添加 GeoSpark 作为他们项目的依赖项，如前一节所述。
**启动 SparkSession:**Spark 或 Sedona 中的任何 SQL 查询都必须由 SparkSession 发出，Spark session 是集群的中央调度器。要启动 SparkSession，用户应使用如下代码:
var sparkSession = SparkSession.builder().appName(“GeoSparkExample”)// Enable GeoSpark custom Kryo serializer.config(“spark.serializer”, classOf[KryoSerializer].getName).config(“spark.kryo.registrator”, classOf[GeoSparkKryoRegistrator].getName).getOrCreate()
注册 SQL 函数: GeoSpark 在 Spark 的 catalyst 优化器中增加了新的 SQL API 函数和优化策略。为了启用这些功能，用户需要使用如下代码将 GeoSpark 显式注册到 Spark 会话。
GeoSparkSQLRegistrator.registerAll(sparkSession)
创建一个几何类型列: Apache Spark 提供了一些格式解析器，可以将数据从磁盘加载到 Spark DataFrame(一个结构化的 RDD)中。在获得数据帧之后，想要运行空间 SQL 查询的用户必须首先在该数据帧上创建几何类型列，因为在关系数据系统中每个属性都必须有一个类型。这可以通过一些构造函数来完成，例如 ST_GeomFromWKT。在这一步之后，用户将获得一个空间数据框架。以下示例显示了该函数的用法。
SELECT ST_GeomFromWKT(wkt_text) AS geom_col, name, address
FROM input
**转换坐标参考系:**与 RDD API 类似，空间 SQL APIs 也提供了一个函数，即 ST_Transform，用于转换空间对象的坐标参考系。它的工作原理如下:
SELECT ST_Transform(geom_col, “epsg:4326", “epsg:3857") AS geom_col
FROM spatial_data_frame
编写一个空间范围查询: GeoSpark 空间 SQL APIs 有一组谓词，用于评估空间条件是真还是假。ST_Contains 是一个经典函数，它将两个对象 A 作为输入，如果 A 包含 B，则返回 true。在给定的 SQL 查询中，如果 A 是单个空间对象，B 是列，则这将成为 GeoSpark 中的空间范围查询(请参见下面的代码)。
SELECT *
FROM spatial_data_frame
WHERE ST_Contains (ST_Envelope(1.0,10.0,100.0,110.0), geom_col)
**编写一个空间 KNN 查询:**要使用 SQL APIs 执行空间 KNN 查询，用户需要首先计算查询点和其他空间对象之间的距离，按升序排列距离，并取前 K 个对象。以下代码查找点(1，1)的 5 个最近邻点。
SELECT name, ST_Distance(ST_Point(1.0, 1.0), geom_col) AS distance
FROM spatial_data_frame
ORDER BY distance ASC
LIMIT 5
**编写一个空间连接查询:**空间 SQL 中的空间连接查询也使用前面提到的评估空间条件的空间谓词。但是，要触发连接查询，空间谓词的输入必须至少包含两个几何类型的列，这两个列可以来自两个不同的数据帧，也可以来自同一数据帧。以下查询涉及两个空间数据帧，一个面列和一个点列。它查找每一对可能的$ < p o l y g o n ， p o i n t polygon，point polygon，point > $以使多边形包含该点。
SELECT *
FROM spatial_data_frame1 df1, spatial_data_frame2 df2
WHERE ST_Contains(df1.polygon_col, df2.point_col)
执行几何运算: GeoSpark 提供超过 15 个 SQL 函数。用于几何计算。用户可以在其空间 SQL 查询中轻松调用这些函数，GeoSpark 将并行运行该查询。例如，获取每个空间对象的面积的一个非常简单的查询如下:
SELECT ST_Area(geom_col)
FROM spatial_data_frame
系统还提供空间对象的聚集功能。它们通常将数据帧中的所有空间对象作为输入，并产生单个值。例如，以下代码计算数据框中所有面的并集。
SELECT ST_Union_Aggr(geom_col)
FROM spatial_data_frame
通过 Zeppelin 笔记本与 GeoSpark 互动
尽管 Spark 在每个版本中都捆绑了交互式 Scala 和 SQL shells，但这些 shell 并不用户友好，并且无法进行复杂的分析和图表。数据科学家倾向于使用图形界面交互式地运行程序和绘制图表。从 1.2.0 开始，GeoSpark (Apache Sedona)提供了一个为 Apache Zeppelin 基于 web 的笔记本量身定制的氦插件。用户可以在 Zeppelin web 笔记本上执行空间分析，Zeppelin 会将任务发送到底层 Spark 集群。
用户可以在 Zeppelin 笔记本上创建一个新段落，并用 Scala、Python 或 SQL 编写代码，与 GeoSpark 进行交互。此外，用户可以点击界面上的不同选项，并要求 GeoSpark 在查询结果上呈现不同的图表，如条形图、折线图和饼图。例如，Zeppelin 可以将以下查询的结果可视化为条形图，并显示美国每个县的地标数量。
SELECT C.name, count(*)
FROM US_county C, US_landmark L
WHERE ST_Contains(C.geom_col, L.geom_col)
GROUPBY C.name
另一个例子是找到美国每个县的面积，并将其可视化在条形图上。相应的查询如下。这实际上利用了 GeoSpark 中提供的几何功能。
SELECT C.name, ST_Area(C.geom_col) AS area
FROM US_county C
Apache Sedona 如何处理空间数据洪流
空间数据划分
此外，空间 rdd 配备了分布式空间索引和分布式空间分区来加速空间查询。所采用的数据分区方法适合集群中的空间数据处理。空间 rdd 中的数据根据空间数据分布进行分区，并且附近的空间对象很可能被放入同一个分区中。空间分区的效果是双重的:(1)当运行以特定空间区域为目标的空间查询时，GeoSpark 可以通过避免对空间上不接近的分区进行不必要的计算来加速查询。(2)它可以将空间 RDD 分割成多个数据分区，每个分区具有相似数量的记录。这样，系统可以确保负载平衡，并避免在集群中执行计算时掉队。
空间索引
Sedona 使用分布式空间索引来索引集群中的空间 rdd。这个分布式索引由两部分组成(1)全局索引:存储在主机上，在空间分区阶段生成。它索引空间 rdd 中分区的包围盒。拥有这样一个全局索引的目的是修剪那些保证没有合格的空间对象的分区。(2)局部索引:建立在空间 RDD 的每个分区上。因为每个本地索引只作用于它自己分区中的数据，所以它可以有一个很小的索引大小。给定一个空间查询，空间 RDD 中的本地索引可以加速并行查询。
空间 RDD 定制串行器
Sedona 为空间对象和空间索引提供了定制的序列化程序。所提出的序列化器可以将空间对象和索引序列化为压缩的字节数组。该序列化程序比广泛使用的 kryo 序列化程序更快，并且在运行复杂的空间操作(例如空间连接查询)时占用的内存更少。将空间对象转换为字节数组时，序列化程序遵循 Shapefile 的编码和解码规范。
序列化器还可以序列化和反序列化本地空间索引，如四叉树和 R 树。对于序列化，它使用深度优先搜索(DFS)按照预先排序策略遍历每个树节点(首先写入当前节点信息，然后写入其子节点)。对于反序列化，它将遵循序列化阶段使用的相同策略。反序列化也是一个递归过程。当序列化或反序列化每个树节点时，索引序列化程序将调用空间对象序列化程序来处理单个空间对象。
结论
总之，Apache Sedona 为数据科学家处理大规模地理空间数据提供了一个易于使用的界面。目前，该系统支持 SQL、Python、R 和 Scala 以及许多空间数据格式，例如 ShapeFiles、ESRI、GeoJSON 和 NASA 格式。以下是 GitHub 资源库的链接:
Apache Sedona(孵化)是一个用于处理大规模空间数据的集群计算系统。塞多纳扩展阿帕奇…
sedona.apache.org](http://sedona.apache.org)
GeoSpark 有一个活跃的小型社区，由来自工业界和学术界的开发人员组成。您还可以在此尝试更多编码示例:
Apache Sedona(孵化)是一个用于处理大规模空间数据的集群计算系统。塞多纳扩展阿帕奇…
github.com](https://github.com/apache/incubator-sedona/)
如果您有更多问题，请随时在 Twitter 上给我发消息
地理空间冒险。第一步:匀称。
使用 KeplerGL 生成的图像
快速浏览使用 Shapely 库在 Python 中处理几何对象的基础知识。
这是一系列文章中的第一篇，总结了过去几年中使用 PropTech 处理地理空间数据的一些关键成果。这些将会有所有的东西——地理空间数据集、几何形状、栅格文件、地图、可视化。从最基础的开始，在后面的文章中，朝着更有趣和更有挑战性的方向努力…
让我们从介绍身材匀称的开始吧。毫无疑问，这是我最喜欢的 Python 库之一——非常重要，对于你最终从事的任何几何/地理相关的工作都是绝对必要的。
该库允许您使用三种主要类型的几何对象:点、线串和多边形+几何图形集合(如果您想要组合它们)。还有一堆其他的——线性环、多点、多多边形等。，但是现在这些就够了，这些方法是可以移植的。
安装
非常标准的安装，使用 pip。如果你像我一样使用 jupyter，你可以直接运行
!pip install shapely
请注意，GeoPandas 正在使用 Shapely，因此如果您已经安装了 GeoPandas，您可能已经有了 Shapely 的最新版本。您可以在导入 Geopandas 后执行一些 Shapely 操作，而无需单独导入，但是，如果您想要直接处理点和面对象，您仍需要首先加载它们。关于 Geopandas 已经说得够多了，但是我们将在下一篇文章中更详细地讨论它。
安装完成后，将其导入笔记本并加载主要几何图形类型:
import shapely
from shapely.geometry import Point, Polygon, LineString, GeometryCollection
import numpy as np
我也在进口 numpy，因为我太喜欢了。不过说真的，我经常发现自己在 shapely 对象和等价的 numpy 坐标数组之间来回跳跃，因为 numpy 允许您以显式矢量化形式更快地进行一些操作，所以从一开始就了解这两者之间的联系是一个好主意。
点状物体
顾名思义，这只是二维平面上的一个点，以一对坐标为特征。
Shapely 的一个超级方便的功能是——它允许您查看所有的几何对象，而不必求助于任何图形包。请注意，无论对象在坐标系中的位置如何，当您想要查看它时，它总是以对象为中心。
pt = Point(10, 10)
pt1 = Point(100, 101)
您还可以显示对象的字符串表示，只需将 str()括起来，或者将其转换为坐标的 numpy 数组。正如我提到的，我发现后者特别有用，因为我经常发现自己在处理大型几何对象阵列(例如来自 OSM 的超过 600 万个建筑多边形)，如果我想以矢量格式进行计算，numpy 是绝对不可替代的。
还有一个将这个字符串表示加载回几何格式的方法，当您必须加载以非几何格式存储的数据(例如从 csv 文件)时，这将非常方便。
如果您想要快速查看几个对象，并查看它们如何相互缩放，您所要做的就是将它们转换为几何体集合:
其他一些简便的方法——距离、坐标采集。
在[8]中:
pt.distance(pt1)
Out[8]:
127.98828071350908
在[9]中:
pt.x, pt.y, pt.xy
Out[9]:
(10.0, 10.0, (array('d', [10.0]), array('d', [10.0])))
在继续学习线条之前，还有最后一件事。所有 shapely 对象都有一个. name 属性。例如，当您要将 GeoPandas 或 Pandas 数据帧中存储的大型集合中的每个多边形转换为一组较小的多边形(如格网)时，并且希望有一种简单的方法将它们与原始多边形相关联时，这将非常有用。
在[10]中:
pt.name = 'My Point'
pt.name
Out[10]:
'My Point'
线串
LineStrings 以非常相似的方式启动，只是这次我们有一个元组列表，而不是单个元组。它们可以交叉并多次通过相同的点，但是，不建议使用后者，因为它会对性能产生负面影响，您最好将它们拆分为单独的组件。请注意，点的顺序很重要，因为它决定了您通过它们的顺序(这同样适用于多边形，您将在下面看到)。
ln = LineString([(0, 1), (20, 100), (100, 3), (120, 102), (200, 5)])
与点一样，可以将 LineString 对象转换为点坐标数组。这里保留了顺序，因此这可用于快速获得第一个和最后一个点的坐标——例如，这对于构建表示道路网络的树对象很方便。
在[13]中:
np.array(ln)
Out[13]:
array([[ 0., 1.],
[ 20., 100.],
[100., 3.],
[120., 102.],
[200., 5.]])
或者，如果您希望首先使用与创建 LineString 相同的元组列表表示:
在[14]中:
list(ln.coords)
Out[14]:
[(0.0, 1.0), (20.0, 100.0), (100.0, 3.0), (120.0, 102.0), (200.0, 5.0)]
您也可以只分离出 X 坐标或 Y 坐标(当然，您也可以使用 numpy 来完成):
在[15]中:
list(ln.xy[0]), list(ln.xy[-1])
Out[15]:
([0.0, 20.0, 100.0, 120.0, 200.0], [1.0, 100.0, 3.0, 102.0, 5.0])
快速浏览一下图形表示以及我们之前创建的点:
计算点到线的距离、线上的点投影(从起点沿着线的距离)和线的长度非常简单:
在[17]中:
pt.distance(ln)
Out[17]:
8.01980198019802
在[18]中:
ln.project(pt), ln.length
Out[18]:
(10.801980198019802, 453.46769176178475)
在[19]中:
list(ln.interpolate(ln.project(Point(1, 1))).coords)
Out[20]:
[(0.039211841976276834, 1.1940986177825703)]
请注意，如果点在线上的投影恰好位于定义的区域之外，则距离将计算到最近的线端。如果你在实际的投影之后，你将需要做一些额外的几何。最简单的方法是延伸直线的第一段和最后一段，并仍然使用相同的投影方法。
线交点也很简单，即使你有多个交点。
结果是一个多点对象，您可以对其进行迭代，就像这是一个常规列表一样，将点对象作为可迭代对象:
在[21]中:
str(ln.intersection(LineString([(0, 0), (200, 100)])))
Out[22]:
'MULTIPOINT (72.55474452554745 36.27737226277372, 110.561797752809 55.28089887640449, 144.5255474452555 72.26277372262774)'
在[22]中:
[np.array(a) for a in ln.intersection(
LineString([
(0, 0),
(200, 100)
])
)]
Out[22]:
[array([72.55474453, 36.27737226]),
array([110.56179775, 55.28089888]),
array([144.52554745, 72.26277372])]
多边形
毫不奇怪，创建一个与创建一个 LineString 非常相似，并且与 LineString 一样，列出点的顺序很重要。
多边形内部可以有洞，这些洞的定义遵循一个简单的规则:多边形([多边形坐标列表]，[洞列表])，其中每个洞本身就是一个多边形。请注意，表示洞的多边形必须完全位于原始多边形内部，或者只能在一个地方接触到原始多边形。
请注意，多边形不能像线串和点那样直接转换成点的集合。相反，我们必须处理它们的外部和内部轮廓(它们本身就是线条，也就是一个自身循环的线条)。对于内部边界，我们得到一个迭代器:
在[25]中:
np.array(poly.exterior)
Out[25]:
array([[0., 0.],
[0., 1.],
[1., 1.],
[1., 0.],
[0., 0.]])
在[26]中:
[np.array(a) for a in poly.interiors]
Out[26]:
[array([[0.5, 0.5],
[0.5, 0.6],
[0.6, 0.6],
[0.6, 0.5],
[0.5, 0.5]]), array([[0.1, 0.1],
[0.1, 0.5],
[0.5, 0.5],
[0.5, 0.1],
[0.1, 0.1]])]
几个额外的方法作为一个追赶者:
多边形的质心(或者，类似地，线对象)
在[27]中:
np.array(poly.centroid)
Out[27]:
array([0.53795181, 0.53795181])
检查点是否在多边形内的两种不同方法:
在[28]中:
Point(0.7, 0.7).within(poly), poly.contains(Point(0.7, 0.7))
Out[28]:
(True, True)
你可以很容易地将线串和点对象转换成多边形通过在它们周围添加一个缓冲区
注意-这仍然是一个多边形(不是严格意义上的圆)，在这种情况下，由 66 个点组成。
在[30]中:
len(np.array(Point(0.7, 0.7).buffer(1).exterior))
Out[30]:
66
如果需要，您可以通过在缓冲区中传递分辨率参数来提高默认分辨率(16)的精度:
在[31]中:
len(np.array(Point(0.7, 0.7).buffer(1, resolution=32).exterior))
Out[31]:
130
缓冲区允许您轻松地关联多边形和其他几何体对象。例如，如果你想用一条特定宽度的线切割一个多边形
或者切割它而不损失任何表面积:
最后，超级有用的 bounds metod，可以应用于多边形、线串和多对象，返回最小边界矩形的左下角和右上角的坐标，与 XY 轴对齐。
在[34]中:
Point(10, 10).buffer(5).bounds
Out[34]:
(5.0, 5.0, 15.0, 15.0)
今天到此为止。下一站熊猫大战地质熊猫！
本系列其他帖子:
地理空间历险记。第四步。魔法的色彩或者如果我没有看到它——它就不存在
地理空间冒险。第二步:熊猫大战地质熊猫
使用 KeplerGL 生成
探索熊猫和 GeoPandas 之间的差异，并介绍一些真实世界的地理空间数据集，我们将在本系列的稍后部分进行分析。
在我的 第一篇帖子 中介绍了 shapely 之后，是时候看看一些有趣的地理数据集了，为了做到这一点，我们不可能没有熊猫，更具体地说是地理熊猫。我假设你以前遇到过熊猫，会试着强调两者之间的一些区别，更重要的是，你如何将一个转换成另一个。
如果您以前从未使用过 GeoPandas，我们将从头开始——继续安装它，特别是如果您想自己按照这篇文章中的步骤操作数据。你只需要一个值得信任的老皮普。
!pip install geopandas
如果你愿意(你绝对应该这样)，你可以在这里阅读这些文件。
让我们从导入熊猫和 GeoPandas 开始。
import pandas as pd
import geopandas as gpd
现在，在我们对这两个进行任何操作之前，我们需要一个数据集。在本文的最后有一些有趣的地理特定数据集的链接(一定要查看它们)，但让我们从相对简短但足够现实的东西开始:英国地方当局边界。我喜欢这个数据集的原因是，一方面，它只有 380 条记录，另一方面，它使用的边界和形状一点也不简单，所以它很好地说明了您可能会在现实生活中使用的数据集。我还认为这是对英国数据进行任何地理空间分析的一个很好的起点，因为它允许你将所有东西分成易于管理的小块。相信我——有时候真的很有帮助。
浏览数据集。英国地方当局界限。
所以，事不宜迟，数据集在这里存放。点击你右边的下载按钮。您将获得一个压缩文件，文件名简短而甜美，大约 35mb(在撰写本文时，他们偶尔会更新)。继续解压缩，您将得到一个如下所示的文件夹:
我们要的是形状文件，带*的那个。shp 扩展。
据我所知，普通的熊猫不能处理这种格式，所以这是熊猫和地理熊猫的第一个区别，当然后者可以。让我们打开它:
la = gpd.read_file(
'Downloads/Local_Authority_Districts__December_2017__Boundaries_in_Great_Britain-shp/Local_Authority_Districts__December_2017__Boundaries_in_Great_Britain.shp'
)
确保您使用的是与您的系统相关的路径。如果你不确定你的笔记本默认使用的是哪个文件夹，只需在其中一个单元格中运行即可——你将获得当前的目录内容，因此你应该能够从这里找到下一步。如果你需要进入一个文件夹，使用<…/>。
让我们来看看数据是什么样子的。所有标准的 Pandas 命令都有效，所以我们可以这样做:
la.head()
得到
不错！好吧，这到底是什么意思？
让我们一栏一栏地来看:
objectid 相当简单，坦率地说，对我们来说不是很有用。
lad17cd 为当地权限代码。这实际上非常方便，因为它与 ONS 的超级输出区域等一起使用，所以这给了我们一个很好的参考系统，可以在不同的地理分区和相应的数据集之间来回移动。非英国人——我很抱歉在这里表现得如此党派化，事实是，大多数国家看待事物的方式都有些相似，所以即使你发现这不是直接的兴趣所在，这仍然是可以转移和有用的。
lad17nm: 地方机关名称。其中有些很诡异很奇妙，绝对值得一探究竟。没有吗？那就只有我…
lad17nmw，里面好像都是看起来友好的<没有> s(抱歉，忍不住)。这一列也有名称，而“ w ”有点泄露了这一点——这是威尔士的地方当局名称，因此前五个记录中没有值，因为这些都是英格兰的地方当局。如果你想知道的话，我们要到 317 号才能离开英国。
bng_e 和bng _ n——邪恶的双胞胎。说真的，这些是东/北——我稍后会再多谈一点。基本上，坐标显示这个地区在地图上的位置。
长和横——这次多了一些坐标——经纬度。
st_lengths —边界长度(也以米为单位)。
最后，这是我们真正需要的！请敲鼓。
几何图形-这是实际的多边形，或者在某些情况下，是多重多边形。不要让我开始谈论苏格兰群岛…Geometry 是 GeoPandas 表的关键属性，许多使用 geo pandas 表的应用程序基本上都假设这个列在那里，并且它被准确地称为“geometry”。所以，不管出于什么原因，如果你想给它取一个更友好的名字，请三思。
在我们进一步深入之前，因为我们要比较 Pandas 和 GeoPandas，我们可以用常规的 Pandas 格式创建这个表的副本。
la_pd = pd.DataFrame(la)
la_pd.head()
你会注意到这张新桌子看起来完全一样。不仅如此，几何表中的对象仍然非常完整:
到目前为止一切顺利。我们还可以像在 Pandas 中一样查看数据帧的所有统计数据，这两个表看起来是一样的，所以我将只在这里显示 GeoPandas 的结果(您必须相信我，不是吗？好的，在您自己的笔记本上运行它，以便检查)。
注意，这个函数只给出数字数据列的统计信息。
在[9]中:
la.info()
Out[9]:
<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 380 entries, 0 to 379
Data columns (total 11 columns):
objectid 380 non-null int64
lad17cd 380 non-null object
lad17nm 380 non-null object
lad17nmw 22 non-null object
bng_e 380 non-null int64
bng_n 380 non-null int64
long 380 non-null float64
lat 380 non-null float64
st_areasha 380 non-null float64
st_lengths 380 non-null float64
geometry 380 non-null geometry
dtypes: float64(4), geometry(1), int64(3), object(3)
memory usage: 32.7+ KB
最后:
在[10]中:
la.memory_usage()
Out[10]:
Index 80
objectid 3040
lad17cd 3040
lad17nm 3040
lad17nmw 3040
bng_e 3040
bng_n 3040
long 3040
lat 3040
st_areasha 3040
st_lengths 3040
geometry 3040
dtype: int64
深部热疗
到目前为止，一切顺利。是时候发现另一个不同之处了。此命令不适用于您的熊猫数据帧:
在[11]中:
la.crs
Out[11]:
<Projected CRS: EPSG:27700>
Name: OSGB 1936 / British National Grid
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E
- bounds: (-9.2, 49.75, 2.88, 61.14)
Coordinate Operation:
- name: British National Grid
- method: Transverse Mercator
Datum: OSGB 1936
- Ellipsoid: Airy 1830
- Prime Meridian: Greenwich
哇哦。！！那是什么意思？！还记得我答应过要扩展北/东吗？是时候看看坐标参考系统了。你可以在这里阅读一个更专业的 GeoPandas 作为具体的解释。简而言之，定义坐标系并不像听起来那么简单。这里的主要问题是地球有不平坦的大胆。不仅如此——它甚至不是球形的，所以在表面上定义一个精确的点可能很棘手(如果你不相信我——试着用一张方面纸包裹橄榄球)。更有趣的是，因为它太大了——有时假装它是平的是有道理的，因为在规模上甚至大到足以覆盖整个国家，如英国，这不会引入足够的误差来真正担心。因此，您可能会遇到两个主要的参考系统——经度和纬度(伪 3D)和北距/东距(2D)。我不得不承认，在大多数情况下，我更喜欢和后者一起工作。原因是——这很容易，而且我很懒。不过说真的，北距和东距都是用米来表示的，并且测量的是直接的距离。所以，如果你想知道事物之间的距离，你所要做的就是得到 Xs(东距)和 Ys(北距)之间的差异，应用毕达哥拉斯，就大功告成了。请注意，我们在之前的帖子中谈到的 Shapely 库也内置了 CRS 的概念，但是，我们并不真的需要使用它，因为 GeoPandas 为我们提供了所有可能需要的高精度弹药。
CRS 转换
GeoPandas 中有许多 CRS 系统。上表中使用的这个特别的是“epsg:27700”——如果你计划使用英国数据，这是你会经常使用的东西，所以你要记住它，要小心。在英国，使用 lat/long 的另一种方法是“epsg:4326”。GeoPandas 为我们提供了一种简单的转换方法。我们要做的就是跑:
la_4326 = la.to_crs("epsg:4326")
la_4326.head()
获得:
眼尖的人会发现，这个表中唯一的变化是我们的多边形/多多边形对象列，现在每个对象都由引用纬度和经度的点组成。您仍然可以看到与我们之前看到的相同的图形表示，它看起来应该是相同的。我会让你来验证的。
保存和读取非地理特定格式(csv)。几何对象转换。
在我们在这些之上进行一些计算之前，让我们看一下保存到 csv 和从 Pandas 转换到 GeoPandas。在 Pandas 的较新版本中，您应该能够将数据帧直接保存到 csv 文件中，尽管最后一列中有几何对象，即
la.to_csv('Downloads/la.csv', compression='gzip')
会成功的。更重要的是，结果文件只有 39.9mb，而 shapefile 只有 63.7mb(当然，我们已经在这里应用了压缩)。然而，这是有代价的。据我所知，您不能直接从 GeoPandas 读取 csv 文件，所以您必须将其作为普通的数据帧加载回来。只需运行:
la_new = pd.read_csv('Downloads/test.csv', compression='gzip')
乍一看，没什么变化:
la_new.head()
本质上，我们有一个单独的额外列，通过运行
la_new = la_new[la_new.columns[1:]]
然而，这还不是全部，如果我们像以前一样尝试查看我们的一个几何对象，我们会得到非常不同的结果:
在[17]
la_new['geometry'].iloc[0]
Out[17]:
'MULTIPOLYGON (((447213.8995000003 537036.1042999998, 447228.7982999999 537033.3949999996, 447233.6958999997 537035.1045999993, 447243.2024999997 537047.6009999998, 447246.0965 537052.5995000005, 447255.9988000002 537102.1953999996, 447259.0988999996 537108.8035000004, 447263.6007000003 537113.8019999992, 447266.1979 537115.6015000008, 447273.1979999999 537118.6007000003, 447280.7010000004 537120.1001999993, 447289.3005999997 537119.6004000008, 447332.2986000003 537111.5026999991, 447359.5980000002 537102.3953000009, 447378.0998999998 537095.0974000003, 447391.0033999998 537082.9009000007, 447434.6032999996 537034.5046999995, 447438.7011000002 537030.9956999999, 447443.7965000002 537027.6966999993...
好吧，我在这里作弊，在现实中，输出比这长得多。实际上，我们所有的几何对象都被转换成了字符串。然而，这一切并没有失去，因为 shapely 为我们提供了一种将它们转换回来的方法。除了在这里加载 shapely，我还需要更快的库。通过使用直接应用方法，您可以不使用它，但除此之外，swifter 为您提供了一个很好的进度条和时间估计(它也可以提高性能，但这超出了本文的范围)。
import swifter
from shapely import wkt
la_new['geometry1'] = la_new['geometry'].swifter.apply(lambda x: wkt.loads(x))
这里发生的事情是，我将转换函数应用于几何列的每个元素，并将输出存储在新列中。在我的笔记本电脑上，在我们的 380 记录长的数据集上，这需要 12 秒。我在非常非常大的数据集上做过，公平地说，还不算太差。
然后，我们可以通过运行以下命令来抽查新元素是否确实是几何图形对象
la_new['geometry1'].iloc[0]
但是我们还没有完成。事实上，如果我们运行 la_new.info()，我们会得到以下结果:
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 380 entries, 0 to 379
Data columns (total 12 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 objectid 380 non-null int64
1 lad17cd 380 non-null object
2 lad17nm 380 non-null object
3 lad17nmw 22 non-null object
4 bng_e 380 non-null int64
5 bng_n 380 non-null int64
6 long 380 non-null float64
7 lat 380 non-null float64
8 st_areasha 380 non-null float64
9 st_lengths 380 non-null float64
10 geometry 380 non-null object
11 geometry1 380 non-null object
dtypes: float64(4), int64(3), object(5)
memory usage: 35.8+ KB
新列显示为对象类型，而不是几何图形类型。但是，这不会阻止我们将其转换为地理数据框架。我们所要做的就是去掉我们不再需要的额外的 geometry 列，并将 geometry1 重命名为 geometry(记住—必须这样命名),然后(！)这一点很重要，在转换为地理数据框架后，我们必须为其设置 crs，
la_new_geo = la_new.drop(columns=['geometry']).rename(columns={'geometry1': 'geometry'})
la_new_geo = gpd.GeoDataFrame(la_new_geo)
la_new_geo.crs = 'epsg:27700'
瞧啊。la_new_geo.info()为我们提供了:
<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 380 entries, 0 to 379
Data columns (total 11 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 objectid 380 non-null int64
1 lad17cd 380 non-null object
2 lad17nm 380 non-null object
3 lad17nmw 22 non-null object
4 bng_e 380 non-null int64
5 bng_n 380 non-null int64
6 long 380 non-null float64
7 lat 380 non-null float64
8 st_areasha 380 non-null float64
9 st_lengths 380 non-null float64
10 geometry 380 non-null geometry
dtypes: float64(4), geometry(1), int64(3), object(3)
memory usage: 32.8+ KB
基本计算示例
最后，对我们新的闪亮地理数据框架进行一些快速计算:
验证长度和面积:
la_new_geo['area']=la_new_geo['geometry'].swifter.apply(lambda x: x.area)
la_new_geo['length'] = la_new_geo['geometry'].swifter.apply(lambda x: x.length)
la_new_geo[['lad17cd', 'lad17nm', 'st_areasha', 'st_lengths', 'area', 'length']].head()
或者，如果我们想要验证整个列，我们可以这样做:
(la_new_geo.st_lengths == la_new_geo.length).all()
您会注意到，由于舍入误差，这实际上会返回 False。例如，在第一条记录中，st_lengths 的值是 71，707.455227013，length 的值是 7 . 1676767771
因此，要进行检查，我们可以以设定的精度运行它:
(la_new_geo.st_lengths.round(4)==la_new_geo.length.round(4)).all()
在任何情况下，当处理地理对象时，考虑到现有数据的质量，移动超过 25 厘米的精度可能是不合理的，对于大多数应用程序，您可能可以将其限制在最近的米。
让我们也看看哪个区域具有最复杂的多边形/多多边形。我们将通过描述它所需的点数来判断。本质上，这意味着计算外部和内部边界的点:
la_new_geo['point_count'] = la_new_geo['geometry'].swifter.apply(
lambda x: np.sum([len(np.array(a.exterior)) + np.sum([len(np.array(b)) for b in a.interiors]) for a in x]) if x.type == 'MultiPolygon'
else len(np.array(x.exterior)) + np.sum([len(np.array(a)) for a in x.interiors])
)
好了，这个有点复杂，我们来看看是怎么回事。首先，我们正在寻找多多边形对象，如果我们有一个，我们必须迭代通过每个多边形，将其转换为外部边界，然后 numpy 数组获得一个坐标数组。我们取数组的长度。我们还获得了内部边界的列表，因此我们还必须遍历它们，将它们转换为数组并取长度。类似地，如果我们处理一个简单的多边形，我们减少一级迭代，所以只有一个外部边界，但是我们仍然需要迭代潜在的多个洞。
然后我们对这些值进行排序，瞧，我们得到了前十名:
la_new_geo[
['lad17cd', 'lad17nm', 'point_count']
].sort_values('point_count', ascending=False).head(10)
这就是为什么每个人都如此热爱苏格兰。
没有吗？那就只有我了…
下一个帖子是关于多边形相互匹配，人类最酷的发明——R 树(Erm…算是)。嗯……好吧，也许还有一两张照片。
最后是承诺的…
数据集
世界各地的 OSM 数据集:http://download.geofabrik.de/
英国的 OSM 数据集:http://download.geofabrik.de/europe/great-britain.html
ONS 开放数据:https://www . ordnance survey . co . uk/opendata download/products . html
NOMIS(劳动力市场/人口普查数据):
-这些通常包含特定区域的数据，然后您可以通过使用它们的边界形状文件(通过 ID 匹配它们)将这些数据与地图相关联
我相信还有很多很多，请在评论中随意添加任何有趣的地理空间数据集(最好是免费使用的)。
回头见…
还在这个系列:
地理空间历险记。第四步。魔法的颜色或者如果我没有看到它——它就不存在
地理空间冒险。第四步。魔法的颜色，或者如果我看不到它，它就不存在。
使用 KeplerGL 生成
快速浏览几何/地理对象(建筑多边形、道路、水路等)的可视化工具。).
从 继续第三步 。现在我们已经有了约克郡东骑马区的所有(OSM 已知的)建筑的数据集，让我们看看可视化数据的方法。
简单回顾一下——在上一篇文章中，我们观察了英国的地方政府多边形，并把它们与 OSM 建筑多边形结合起来。然后，我们讨论了一种将所有建筑多边形归属于相应的地方当局多边形的方法，并集中讨论了约克郡的 East Riding，这是由它的形状之美和与赫尔的邻近性所激发的。然后，我们继续创造一个非常不恰当的结果视觉化，我们即将纠正它。这篇文章是关于在地图上查看几何/地理数据的不同方法。
在我们能正确开始之前，我们确实需要一些小步骤。
首先，让我们加载我们在上一篇文章中使用的所有库:
import geopandas as gpd
import swifter
import numpy as np
from shapely.geometry import Polygon, GeometryCollection
在此之上，我们仍然有我们的数据框架和属性化的建筑，我们称之为“bd_ery”。让我们快速浏览一下:
bd_ery.head()
用 GeoPandas 绘图
我们将逐步进行，首先，让我们看看地理数据框架本身可以做些什么，因为除了显示几何对象之外，它还具有一些额外的绘图功能。如果你是 GeoPandas 的新手——看看本系列的第二篇文章
我们仍然需要使用 MatPlotLib，所以让我们继续加载它。
import matplotlib.pyplot as plt
%matplotlib inline
还有一件事——策划整个地方当局将使我们所有美丽的建筑形状变得非常小，这将是不可原谅的。因此，相反，我们将首先利用我们创建的图块来对它们进行属性化(详见上一篇文章)。我们不想处理它们的多边形，然而，我们需要的只是一种简单的方法来选择位于特定瓦片中的建筑物。为此，我们将把多边形转换成 id。最简单的方法是获取左下角的坐标(整数形式)并将它们合并成一个 12 位的 id(您也可以考虑散列多边形的字符串表示，但是我们不需要走那么远)。您可以将我们的 12 位 id 保存为字符串或整数(后者占用的内存更少)。我将把它们保存为字符串，因为无论如何它们都必须留在列表中，所以 DataFrame 最终会把它们视为对象。
转换只需一行代码，使用一些基本的列表理解即可完成:
bd_ery['ids']=bd_ery['grid'].swifter.apply(
lambda x: [
''.join(np.array(a.exterior)[0].astype(int).astype(str)) for a in x
]
)
我将使用瓷砖’ 499323426360 '没有任何具体原因，除了它有相当多的建筑，所以它非常适合我们的目的。
fig, ax = plt.subplots(1, figsize=(15, 15))
ax.axes.get_xaxis().set_visible(False)
ax.axes.get_yaxis().set_visible(False)
bd_ery[bd_ery['ids'].apply(lambda x: 1 if '499323426360' in x else 0) == 1].plot(ax=ax)
好吧，这还是有些乏味，不是吗？你肯定能看到一堆建筑形状，但也就这些了。当然，我们可以做得更多一点。
让我们这次使轴可见，让我们也根据建筑的类型给每个建筑分配颜色，让我们在看的时候放一张色图。
网格线，有人吗？哦，好吧，网格线也是…那图例呢？是的，是的，好的…
fig, ax = plt.subplots(1, figsize=(15, 15))
colormap="RdYlBu"
ax.grid()
bd_ery[bd_ery['ids'].apply(lambda x: 1 if '499323426360' in x else 0) == 1].plot(
column='type',
ax=ax,
cmap=colormap,
legend=True
)
那更好…这实际上在合适的人手里是有用的，尤其是考虑到制作一个这样的东西是多么容易。把它贴在你的网络应用或交互式仪表盘上，或者至少把截图粘贴到一个中型博客上，假装你知道你在说什么。
然而，我们还没有完全达到目标。我们可以做得更好。
输入…
叶子
关于 folium 的伟大之处在于，虽然它仍然完全由您的应用程序控制，并且输出仍然可以用于 webapps 和仪表盘，但输出是一个 html 文件，这使用户能够对生成的地图进行大量控制。是的，这是一张地图，我们现在要在实际的地图上画出我们可爱的多边形！
在我们开始多边形之前，让我们看看是否可以自己创建地图。奔跑
!pip install folium
如果你还没有，让我们加载库和底图。
import folium
state_geo = '[http://geoportal1-ons.opendata.arcgis.com/datasets/01fd6b2d7600446d8af768005992f76a_4.geojson'](http://geoportal1-ons.opendata.arcgis.com/datasets/01fd6b2d7600446d8af768005992f76a_4.geojson')
我们还需要传递一个起点，这样 leav 就知道我们要找的是地图的哪个部分。这需要以纬度/经度格式完成，因此我将创建一个迷你地理数据框架，以我们的切片中心作为唯一的点，然后使用 crs 转换方法来获得纬度/经度。这是一种有点复杂的方式，但它确实有效。曾经有一篇 Hannah Fry(BBC4 科学播客 fame——很棒，看看吧)的优秀博客文章，其中有转换代码，你仍然可以在 GitHub 上找到用 Python 写的版本，但我相信博客本身已经不在了。我觉得这是值得一提的，因为我仍然偶尔使用这些代码，并且知道你仍然可以在那里找到它的版本。
from shapely.geometry import Point
tile = gpd.GeoDataFrame(
{'geometry':[Point(499323+1000,426360+1000)]}
)
tile.crs = 'epsg:27700'
tile = tile.to_crs(epsg=4326)
这为我们处理了转换，现在我们要做的就是绘制:
m = folium.Map(location=[tile.geometry.iloc[0].y, tile.geometry.iloc[0].x], zoom_start=15)
m
看起来有些眼熟？
现在让我们画一些多边形。
首先，我们将把图块的内容转换成一个单独的数据帧，并将其转换成 epsg:4326，这样整个事情就更容易理解了:
bd_tile = bd_ery[
bd_ery['ids'].apply(lambda x: 1 if '499323426360' in x else 0) == 1
].to_crs('epsg:4326').reset_index(drop=True)
现在我们可以开始了:
folium.GeoJson(bd_tile['geometry'].to_json()).add_to(m)
m
啊哈！我们快到了！请注意，我们只是在之前的地图上添加了一些东西，并没有重新绘制整个地图。
让我们给它添加一些样式。我要从 Seaborn 那里偷调色板，把它变成一个字典，然后把它传到 leav。
import seaborn as sns
types = list(bd_tile['type'].unique())
color_dict = dict(
zip(
types,
sns.color_palette("Paired", len(types)).as_hex()
)
)
m = folium.Map(
location=[
tile.geometry.iloc[0].y,
tile.geometry.iloc[0].x
],
zoom_start=15
)
folium.GeoJson(
bd_tile[['osm_id', 'type', 'geometry']].to_json(),
style_function = lambda feature: {
'fillColor': color_dict[feature['properties']['type']],
'color': 'black',
'weight': 0.2,
'fillOpacity': 0.7
}
).add_to(m)
m
然后，您可以像下面这样简单地将其保存为 html:
m.save_to('map.html')
还不算太差，但还是有一些弊端。除非你准备创建一个自定义的 html 图例，就像这里的一样，否则你会陷入没有图例的困境。标准方法只允许添加自动生成的颜色图，这对于分类变量来说一点也不好玩。还有另一种方法使用 choropleth 方法创建相同的地图，你可以在这里选择，但是同样没有可定制的图例和工具提示/弹出窗口。后者你可以通过添加圆形物体到多边形质心并混合它们。这些家伙可以有弹出窗口和工具提示，所以它让你的生活变得更容易一点。不过，你可能明白我的意思——还有另一个选择，我认为这可能是最好的方法。
进入
开普勒格尔
你可以在这里阅读。它可以用在 jupyter 笔记本上，这里有一篇非常详细的文章(包括最后的安装说明)。它不是没有一些缺点，因为它不能通过笔记本本身的代码完全控制，地图需要一些手动定制。尽管如此，它仍然可以保存为 html 格式。
让我们看看开普勒的例子。
from keplergl import KeplerGl
map_1 = KeplerGl(height=500)
Height 参数是指 jupyter 中窗口的大小。
现在我们要做的就是将地理数据框添加到地图对象中。请确保不要传递任何带有几何对象的其他列，因为这会使其混淆。
map_1.add_data(
data=bd_tile[['osm_id', 'type', 'geometry']],
name='Building types'
)
map_1
看起来很熟悉，第一眼看上去没什么印象。关键是——从现在开始，您可以非常轻松地进行定制，并将其放在需要的地方。
你所需要做的就是点击左上角带箭头的小方块，进入下面的菜单:
在第一个(层)菜单中，您可以调整对象的颜色和边界，颜色可以由您上传的数据帧的列决定:
交互菜单允许你控制工具提示的内容。右边的菜单允许你切换到 3D 模式，添加图例等等。同样，3D 表示的高度可以由数据框列决定，也可以设置为特定值。
你还可以放大和缩小，移动到其他位置，改变视角等。
这是一个 3D 效果的示例，根据类型进行颜色编码，并显示带有 osm_id 和类型的工具提示:
最后一步——将其保存为 html 格式并分享出去…
map_1.save_to_html(file_name='kepler_example.html')
还不错…唯一的主要缺点是必须手动将最终图像调整到正确的状态，所以这只对 webapp 中的一次性演示或静态数据表示有好处。尽管如此，注意这个领域——这些缺点可能会相对较快地得到纠正。
哦，对了……差点忘了——标题图片是威斯敏斯特所有的建筑，都是由同一个开普勒格尔建造的。这实际上指出了它的另一个优势——它可以处理相当多的数据，因此在地图上绘制超过 20000 个多边形对象根本不是问题。
今天到此为止。下次—我们来看看光栅文件及其分析/解释。
还在这个系列:
地理空间冒险。第五步。离开平地或飞越多边形的海洋。
使用 KeplerGL 生成
使用激光雷达数据(栅格文件)估计 OSM 多边形对象的高度。
顾名思义，是时候让我们做一些真正酷的事情，用“激光”去飞翔了。字面上。我们将了解如何解释激光雷达数据并将其连接到我们的多边形。
在第 4 步中，我们集中精力以各种方式绘制约克郡 East Riding 的特定瓷砖。我们将继续使用同一个 bd_tile 数据集，它是约克郡某个特定位置的建筑物集合(具体来说，是一个 2 公里长的正方形，左上角位于 N:426360，E:499323)。如果您计划遵循本文的步骤，我建议您查看以前的帖子(参见本文末尾的完整链接列表)来重新创建相同的数据集。
我们在这里要尝试和解决的问题是:在公开的街道地图上确定建筑物的高度。
事实上，它不一定是建筑物——它可以是任何我们知道位置和形状的物体。
所以事不宜迟:
激光雷达
首先，我们必须找到数据，因此我们将向环境署寻求帮助。他们进行定期调查，在夜间飞越英国的特定地区，并向地面发射激光。测量返回的光束，结合平面的位置，我们可以计算出光束反射的物体的高度。数据然后被分成数字表面模型 DSM(包括一切)和数字地形模型 DSM(去除地面物体)。这些以光栅文件的形式出现。asc ),分辨率为 0.25 米、0.5 米、1 米、2 米。
并不是所有的英国都在所有的分辨率下可用，存在差距，并且数据是在不同的时间点获取的，它不是每年都更新，所以当然存在(在某些情况下是显著的)缺点。尽管如此，这可能是我们目前能做的最好的事情了，此外，建筑物往往会在一段时间内保持不变，因此对于它们中的大多数来说，2-3 年的旧数据不是问题。
数据可以在这里找到。我们现在正在寻找两块瓷砖:SE92NE 和 TA02NW。原来我们到目前为止处理的区域正好位于多个瓷砖之间，这虽然有点烦人，但为我们提供了一个方便的机会来完成将这些瓷砖拼接在一起的过程。
我们的目标是合成 DTM 和合成 DSM 数据集，在这些数据集内，我们将达到 50 厘米的分辨率。每个 asc 文件的前六行包含有关其位置和分辨率的信息。随后是一个 2000 乘 2000 的数字方阵，每个数字代表一个 0.5 米乘 0.5 米的像素，并给出它的高度。对于相同的像素，DSM 和 DTM 数之间的差异给出了相应物体距地平面的高度，单位为米。
准备数据
该切换到 jupyter 了，让我们先加载库。除了 geopandas 和 numpy，我们还需要编解码器来读取我们的 asc 文件和 matplotlib/seaborn 以进行一些可视化。
import geopandas as gpd
import numpy as np
import codecsimport matplotlib.pyplot as plt
import seaborn as sns%matplotlib inline
接下来会有一大段代码，做好准备。加载相关的图块，进行一些滤光，然后将它们组合在一起
with codecs.open('LIDAR-DSM-50CM-SE92ne/se9926_DSM_50CM.asc', encoding='utf-8-sig') as f:
X00 = np.loadtxt(f, skiprows=6)X00[X00 < 0] = 0with codecs.open('LIDAR-DSM-50CM-SE92ne/se9927_DSM_50CM.asc', encoding='utf-8-sig') as f:
X01 = np.loadtxt(f, skiprows=6)X01[X01 < 0] = 0with codecs.open('LIDAR-DSM-50CM-SE92ne/se9928_DSM_50CM.asc', encoding='utf-8-sig') as f:
X02 = np.loadtxt(f, skiprows=6)X02[X02 < 0] = 0X_ = np.vstack([X02, X01, X00])with codecs.open('LIDAR-DSM-50CM-TA02nw/ta0026_DSM_50CM.asc', encoding='utf-8-sig') as f:
X10 = np.loadtxt(f, skiprows=6)X10[X10 < 0] = 0with codecs.open('LIDAR-DSM-50CM-TA02nw/ta0027_DSM_50CM.asc', encoding='utf-8-sig') as f:
X11 = np.loadtxt(f, skiprows=6)X11[X11 < 0] = 0with codecs.open('LIDAR-DSM-50CM-TA02nw/ta0028_DSM_50CM.asc', encoding='utf-8-sig') as f:
X12 = np.loadtxt(f, skiprows=6)X12[X12 < 0] = 0X_1 = np.vstack([X12, X11, X10])with codecs.open('LIDAR-DSM-50CM-TA02nw/ta0126_DSM_50CM.asc', encoding='utf-8-sig') as f:
X20 = np.loadtxt(f, skiprows=6)X20[X20 < 0] =0with codecs.open('LIDAR-DSM-50CM-TA02nw/ta0127_DSM_50CM.asc', encoding='utf-8-sig') as f:
X21 = np.loadtxt(f, skiprows=6)X21[X21 < 0] = 0with codecs.open('LIDAR-DSM-50CM-TA02nw/ta0128_DSM_50CM.asc', encoding='utf-8-sig') as f:
X22 = np.loadtxt(f, skiprows=6)X22[X22 < 0] = 0X_2 = np.vstack([X22, X21, X20])X = np.hstack([X_, X_1, X_2])
好了，这里实际发生了什么:我们正在读取每个文件，跳过前 6 行，因为它们包含描述性信息。您可以在文本编辑器中打开它们进行查看(或者使用 print 在笔记本中显示它们)。我们现在正在追求 2000 乘 2000 矩阵。我们用零代替所有的负值。这是因为数据中的缺口是用-9999 值填充的，我们肯定要忽略这些。然后，我们垂直堆叠得到的矩阵，得到 2 公里乘 6 公里的区域，最后水平堆叠这些矩阵，得到一个 6 公里乘 6 公里的区块。
让我们来看看:
dim = (15, 15)
fig, ax = plt.subplots(figsize=dim)
sns.heatmap(X, cmap="binary", square=True, ax=ax)
还不是很鼓舞人心。请注意左上角的大块缺失数据。恐怕对此无能为力。我们还需要加载 DTM 文件(可以在上面的代码块中用 DTM 替换 DSM)。我将结果保存为 X_t。
让我们直接看结果:
看起来更糟，不是吗？让我们来看一下不同之处，将地图的高度比例限制为 20 米，以使所有内容更加突出:
Xdiff = X-X_t
Xdiff[Xdiff < 0] = 0
dim = (15, 15)
fig, ax = plt.subplots(figsize=dim)
sns.heatmap(Xdiff, vmax=20, cmap="binary", square=True, ax=ax)
啊哈！我们有所进展，这实际上开始看起来像一张地图。现在，这个地图左下角的坐标由我们加载的第一个图块给出— N:426000，E:499000。我们需要在左下角的 N:426360，E:499323 处切出一个 2 公里乘 2 公里的正方形。请记住，我们矩阵的元素是从左上角开始索引的，而不是从底部，所以我们需要进一步调整，然后还要记住，我们矩阵中的每个点都是 0.5 米宽，所以我们需要将索引乘以 2，使它们等于北距/东距(因此，我们需要移动 720 和 646，而不是 360 和 323，2 公里的瓷砖边变成 4000 像素)
Xtile = Xdiff[6000-4720: 6000-720, 646: 4646]
画出它，我们会得到与原始瓷砖非常相似的东西:
和…相对
fig, ax = plt.subplots(1, figsize=(15, 15))
ax.axes.get_xaxis().set_visible(True)
ax.axes.get_yaxis().set_visible(True)
ax.grid()
bd_tile.plot(ax=ax)
匹配多边形
下一步是匹配两组值，点和面。我们将要采取的方式与步骤 3 中详细描述的非常相似。我们取每个多边形，用 0.5 米的步长创建一个包围盒的网格。然后，我们可以遍历网格单元，检查它们中的哪些与我们的多边形相交。它们的坐标将与我们的高度矩阵的元素直接关联，因此我们可以简单地选择相关的元素，并获得我们的建筑物的高度轮廓。听起来合理吗？让我们试一试。
加载 shapely Polygon 和 GeometryCollection，我将使用数据集中的第一个多边形:
from shapely.geometry import Polygon, GeometryCollection
bounds = bd_tile['geometry'].iloc[0].bounds
创建网格并将每个元素转换为多边形:
X, Y = np.mgrid[int(bounds[0]): int(bounds[2]) + 1: 0.5,
int(bounds[1]): int(bounds[3]) + 1: 0.5]
grid = list(
map(
list,
list(
zip(
list(
zip(X.flatten(), Y.flatten())
),
list(
zip(X.flatten(), Y.flatten() + 0.5)
),
list(
zip(X.flatten() + 0.5, Y.flatten() + 0.5)),
list(
zip(X.flatten() + 0.5, Y.flatten())
)
)
)
)
)
grid_poly = [Polygon(a) for a in grid]
这与我们在步骤 3 中使用的方法完全相同，只是现在我们使用 0.5 米的步长，而不是 2000 米。
让我们看看我们对多边形的近似程度:
poly = [
a for a in grid_poly if a.intersects(
bd_tile['geometry'].iloc[0]
)
]
g = GeometryCollection(poly + [bd_tile['geometry'].iloc[0]])
g
这是非常整洁的。
现在我们需要从我们的矩阵中得到相应的点。
让我们抓住坐标并把它们转换成矩阵指数。
poly_coords = np.array([np.array(a.exterior)[0] for a in poly])
这给了我们每个左下角的坐标对。现在把它们分成 x 数组和 y 数组:
inds_array = np.moveaxis(poly_coords, -1, 0)
调整以使它们与我们的矩阵坐标一致:
inds_array[0]=(inds_array[0] - 499323) * 2
inds_array[1]=(2000 - inds_array[1] + 426360) * 2
最后，使用花式索引(老实说，这是一个技术术语):
heights = Xtile[
inds_array[1].astype(int),
inds_array[0].astype(int)
]
让我们检查一下我们所选择的实际上看起来像我们的原始多边形:
dim=(15,15)
fig, ax = plt.subplots(figsize=dim)
sns.heatmap(
Xtile[
np.min(inds_array[1].astype(int)): np.max(inds_array[1].astype(int)),
np.min(inds_array[0].astype(int)): np.max(inds_array[0].astype(int))
],
cmap="binary",
square=True,
ax=ax
)
嘣！我会说——完全正确。
现在，理论上，我们得到了一个很好的唯一的数字来显示我们建筑的高度。当然，在现实中，这是极不可能发生的。你可以从上面的图片中看到，我们得到了相当多的分布——倾斜的屋顶，多个建筑块，烟囱，等等。等等。一个不错的猜测是，使用这样的中值，在这个特定的例子中，我们得到了不到 3.6 米。更好的是，我们可以丢弃选择中的零值，这将使我们的中值不到 4 米。对于某些应用程序，您可能想要最大高度，尽管在这种情况下，我们得到 9 米，这很可能是由附近的树驱动的。要做出更明智的选择，您可以查看总体分布并选择特定的百分位数，或者创建与不同身高值对应的总面积相关的条件。另一种选择是排除边界网格销售(通过迭代相关网格销售并检查哪些销售与边界相交来识别它们)。外景法)。
我们可以使高度更谨慎一点，例如将它们四舍五入到最接近的半米，并获得计数，以给我们一个稍微好一点的画面:
h, c = np.unique(np.round((heights * 2), 0)/2, return_counts=True)
dict(zip(h, c))
生产
{0.0: 184,
0.5: 30,
1.0: 19,
1.5: 22,
2.0: 16,
2.5: 70,
3.0: 118,
3.5: 43,
4.0: 56,
4.5: 79,
5.0: 98,
5.5: 92,
6.0: 89,
6.5: 38,
7.0: 7,
7.5: 2,
8.0: 2,
8.5: 1,
9.0: 1}
现在剩下的就是选择你的规则，把我们的步骤包装成一个函数，然后把它应用到数据集中的每个多边形。我要把它留给你去玩，否则我会得到所有的乐趣，那只是贪婪。
本系列到此结束，希望你喜欢它，或者至少发现它在某些方面是有用的。当然，这还不是全部。例如，我们可以看看当没有 OSM 作为后盾时，如何将栅格数据转换为实际的建筑物多边形(他们缺少许多建筑物，因此这不仅仅是一个理论练习)。这是一个很好解决的问题，我可能会回来。但是还有很多很多地理空间(不仅仅是地理空间)项目值得探索，所以请继续关注。
再见，
鐽
本系列的前几篇文章
地理空间冒险。第四步。魔法的颜色或者如果我没有看到它——它就不存在
Python 和 Jupyter 笔记本中的地理空间分析
使用 geopandas 和 kepler.gl 对巴塞罗那自行车租赁服务(bicing)进行地理空间分析。
世界上大多数国家的首都都在使用公共城市自行车服务，这减少了燃料消耗、排放和市中心的拥堵。自行车共享还鼓励体育活动，有助于城市居民的健康。
照片由 Jasser Gómez 在 Unsplash 拍摄
数据摄取
幸运的是，它有一个开放的 CityBikes API ，可以用来实时检查自行车站点的状态(例如，我们可以检查任何站点的空闲位置的数量)。此外，还有一个用于查询 CityBikes API 的 python 包[python-citybikes](https://github.com/eskerda/python-citybikes)。
巴塞罗那是使用 CytyBikes API 公开数据的城市之一。使用上面提到的包，巴塞罗那的自行车服务(又名 bicing )可以使用以下代码进行查询:
client = citybikes.Client()
bicing = citybikes.Network(client, uid='bicing')
bicing.stations
前一个调用以 JSON 格式返回工作站的数据:
{
"company": ["Barcelona de Serveis Municipals, S.A. (BSM)", "CESPA", "PBSC"],
"href": "/v2/networks/bicing",
"id": "bicing",
"location": {
"city": "Barcelona",
"country": "ES",
"latitude": 41.3850639,
"longitude": 2.1734035
},
"name": "Bicing",
"stations": [
{
"empty_slots": 17,
"extra": {
"ebikes": 0,
"has_ebikes": true,
"normal_bikes": 16,
"online": true,
"uid": 361
},
"free_bikes": 16,
"id": "ed25291d0f5edd91615d154f243f82f9",
"latitude": 41.376433,
"longitude": 2.17871,
"name": "PG. DE COLOM (LES RAMBLES)",
"timestamp": "2020-10-16T18:19:06.097000Z"
},
...
]
}
处理完数据后，我们将数据从 JSON 模式转换成表格模式:
自行车站数据框
该 API 不提供历史服务，但是您可以开发自己的脚本来每 5 分钟下载一次数据(Bicing refresh rate)。配套的笔记本显示如何下载和处理数据；在data/bicing_saturday_morning.parquet中，处理后的数据样本也作为拼花文件提供。
巴塞罗那——自行车站
第一个空间操作将是使用从 CityBikes API 下载数据时可用的经度和纬度信息在地图上绘制车站:
df_stations = df[["name","latitude","longitude"]].drop_duplicates()
第一张地图会将所有这些站放在地图上:
stations_map = KeplerGl(height=600, data={"stations": df_stations})
stations_map
巴塞罗那自行车站
kepler.gl很聪明地推断出列latitude和longitude代表了车站的地理坐标。该地图可以显示在 Jupyter 笔记本中，并且可以使用 UI 配置该地图，以在不同的地图样式(平面或卫星)和图层(道路、建筑物等)之间切换。这些 UI 调整可以导出为 python 字典，然后在下次加载地图时使用。查看 kepler.gl 文档了解更多关于定制地图的信息。
细川玉子区—自行车站
CityBikes API 不提供按地区划分的车站信息(只有车站的名称，也就是车站的地址)。但是不用担心，我们仍然可以使用社区提供的公开数据。在这种情况下，我们将使用巴塞罗那开放地理数据存储库，以及那里提供的地区地理位置。
Geopandas 可以阅读。geojson 文件使用:
df_districts = geopandas.read_file("data/districtes.geojson")
对象df_districts类似于一个常规的 Pandas dataframe，但是包含一个名为geometry的列。Geopandas 数据帧可以根据需要包含任意多的地理列，但是我强烈建议每个数据帧只使用一个，并将其命名为geometry。在这种情况下，文件包含 10 个地区(行)和 47 列(带有地区名称、网站、区域等信息)。
巴塞罗那地区数据框架
geopandas 数据帧可以按原样绘制:
巴塞罗那区
或者，我们可以在另一个地图实例中沿着自行车站投影该地区的边界:
district_map = KeplerGl(height=600, data={"stations": df_stations, "districts": df_districts})
巴塞罗那自行车站和区
在上图中，我们可以直观地查看哪些车站位于细川玉子区，但是有没有办法以表格的方式提取这些信息呢？是的，您可以:欢迎使用空间连接。
首先，我们需要将包含站点位置的数据帧转换为 Geopandas 数据帧:
df_stations = geopandas.GeoDataFrame(df_stations.name, geometry=geopandas.points_from_xy(df_stations.longitude, df_stations.latitude, crs='epsg:4326'))
在上面的代码片段中，您可以观察到一个名为crs的参数。总结一下:经纬度数字是有规律的数字，但是当你对它们应用一个坐标参考系 ( CRS )的时候，这些数字就可以投影到一张地图上。参考系比较多，课题比较复杂，我就把我的研究缩减到只有两个 CRS 的(请 GIS 专家原谅我的无礼):
- 通用 WGS 86 (EPSG: 4326) 以经度/纬度表示世界上的位置
- 墨卡托 (EPSG: 3395) 当与距离一起工作时
当我们使用两个坐标系来投影世界地图时，可以理解两个坐标系之间的差异:
使用 WGS86 和墨卡托投影的世界地图
如果您需要有关投影和使用不同参考系统的更多信息，我建议查看管理投影 GeoPandas 文档。
空间连接将返回包含在每个区内的点(=自行车站):
df_stations_districte = geopandas.sjoin(df_stations, df_districts[["NOM","geometry"]], how="inner", op='intersects')
从那里，过滤掉细川玉子区的电台就这么简单:
df_stations_district_gracia = df_stations_district[df_stations_district.NOM=='Gràcia']
格拉西亚区的自行车站
细川玉子区——自行车站“影响区”
影响区域可以显示为以每个站为中心的圆。圆的半径可以是用户定义的参数(例如，300 米)。正如我们之前提到的，为了处理距离，我们需要暂时将参考系统切换到墨卡托:
def buffer_in_mercator(geometry, distance_x_meters):
crs = geometry.crs
return geometry.to_crs(epsg=3395).buffer(distance_x_meters).to_crs(crs)
将函数应用于站点的位置后，我们可以绘制每个自行车站点周围的影响区域:
自行车站影响区域
细川玉子区——骑自行车的最佳区域
在这种情况下，所采用的方法是，地图中的每个区域都可能受到许多站点的影响，因此该区域中可用的自行车数量将是与该区域接触的影响区域中可用自行车数量的总和。
自行车站点区域和统计数据
我们需要计算所有影响区域之间所有交集的并集。这个操作可以使用库shapely来解决(感谢 stackoverflow ):
shapes = areas.to_list()
areas_overlay = list(polygonize(unary_union(list(x.exterior for x in shapes))))
该操作为每个叠加交点创建一个形状(区域)，因此总共生成 68 个形状。如前所述，我们可以绘制图形来更好地理解操作本身:
自行车站区
免费自行车的数量可通过在区域和车站影响区域之间执行空间连接，并将同一区域中的所有行相加来获得:
df_gracia_area_stats = gpd.sjoin(df_areas_overlay,
df_gracia_stats.reset_index(),
how="left", op='intersects')\
.reset_index()\
.rename(columns={'index':'index_left'})
df_gracia_area_stats['free_bikes_area'] = df_gracia_area_stats.groupby("index_left")['free_bikes'].transform(sum)
同样，实例df_gracia_area_stats是一个 geopandas 数据帧，可以投影到地图上，直观地检查 Gràcia 区的自行车可用性:
free_bikes_map = KeplerGl(height=600, data={"areas": df_gracia_free_bikes})
格拉西亚区
与之前相同，可以用表格的方式提取最佳区域:
df_gracia_free_bikes.loc[df_gracia_free_bikes.free_bikes_median.idxmax()]
总之，geopandas和kepler.gl是您数据分析工具箱的强大补充。geopandas基于 Pandas dataframe API，而kepler.gl在其核心支持 Pandas dataframe，因此这两个库的学习曲线都很快，同时好处也很大。
这篇文章的完整代码可以在这里找到——https://github.com/franperezlopez/blog_spatial
安装kepler.gl
安装这个库有点棘手，这方面缺少文档，所以我会尽量简化这一步。首先，这个库似乎不能在 Windows 中工作，所以你需要使用 MacOS、Linux 或 WSL(2)。
其次，这个库需要安装一个 Jupyter 扩展。按照我自己关于如何安装 Jupyter Notebook 的建议，您应该在您的 base 环境中安装这个扩展，执行以下命令:
conda activate base
pip install --user --use-feature=2020-resolver keplergljupyter nbextension install --system --py keplergl
jupyter nbextension enable --system --py keplergl
参考
地理空间聚类:种类和用途
对所有不同种类的集群及其用例的详细回顾
丹尼斯·库默在 Unsplash 上拍摄的照片
地理空间聚类
地理空间聚类是将一组空间对象分组为称为“簇”的组的方法。一个聚类中的对象表现出高度的相似性，而这些聚类尽可能地不相似。聚类的目标是进行概化并揭示空间和非空间属性之间的关系。
让我们用一个小例子来理解空间聚类。
假设你是一家大都市食品配送连锁店的负责人，你想了解顾客的偏好，以便扩大业务规模。因为查看每个客户的详细信息是不可行的，所以你将他们分成不同的组，并为每个组/群制定商业计划。
空间聚类可分为以下五大类型:
1.划分聚类
2.分层聚类
3.模糊聚类
4.基于密度的聚类
5.基于模型的聚类
通过 Locale ，我们致力于让每个在地面上移动资产的企业都可以访问位置数据。让我们深入研究集群的类型，并详细了解每一种类型！
划分聚类
分区聚类是将数据点分离成不重叠的子集(聚类),使得每个数据点恰好在一个子集中。基本上，它通过满足这两个要求 : 将数据分类到 组
1.每个数据点仅属于一个聚类。
2.每个聚类至少有一个数据点。
划分聚类有三种类型:K-means 聚类、K-medoids 聚类/PAM 和 CLARA(分类大型应用)
1.k 均值聚类
K-means 聚类是一种分区方法，这种方法根据属性将数据集分解成一组 K-分区。你可以在这里阅读更多关于 k 的意思。
2.K-medoids 聚类/PAM
K-medoids 聚类是一种类似 K-means 聚类的划分方法*。*med oid 是聚类中与该聚类中所有其他点具有最小相异度的点。查看这个来了解更多关于帕姆的信息。
3.克拉拉
这是 PAM 方法对大型数据集的扩展。CLARA 不是为整个数据集寻找几何图形，而是考虑一个固定大小的小样本数据，并应用该算法为该样本生成一组最佳几何图形。
使用案例:
通常，基于分区的聚类用于查找数据中没有明确标记的组。这有助于将任何新数据点分配到正确的聚类。企业使用基于分区的聚类来划分购买历史、按销售活动分组库存、在健康监控中识别组等。
优步使用聚类进行意图检测，他们将其用于一键式聊天/建议回复系统【来源】
分层聚类
分层聚类是一种类似于基于分区的聚类的聚类方法，但它对数据点进行分类的方式不同。它首先将每个数据点视为一个单独的聚类。然后合并彼此最接近的相似聚类。它不断迭代，直到所有的集群被合并。
两种类型的层次聚类如下:聚集和分裂
凝聚层次聚类
这与计算点/簇的邻近矩阵的简单算法一起工作。在每次迭代中，它合并最近的点/簇，并更新邻近矩阵。这一直持续到形成一个簇或 k 个簇。这可以被认为是一种*“自下而上”*的方法。
分裂层次聚类
这与聚集聚类正好相反。在该方法中，最初，所有数据点被认为属于单个聚类。在每次迭代中，将不相似的数据点从聚类中分离出来。每个分离的数据点被认为是一个单独的聚类。这个过程一直持续到我们有了 K-clusters。这可以被认为是一种【自上而下】的方法。
使用案例:
虽然分层聚类在计算上可能是昂贵的，但是产生直观的结果。由于层次聚类不需要太多的假设，所以当手头的数据不太为人所知时，这是很有用的。分层聚类广泛应用的一个现实场景是在流行病传播期间绘制病毒图，以及在银行业或零售业中进行客户细分。
使用层次聚类将美国参议员分组到他们各自的党派中。红色代表共和党，蓝色代表民主党，黑色代表独立【来源】
模糊聚类
在模糊 c 均值聚类中，我们找出数据点的质心，然后计算每个数据点到给定质心的距离，直到形成的聚类变得恒定。它不同于基于分区的聚类，它允许将数据点部分分类到多个聚类中。
卫星图像的模糊聚类[ 来源
理论上，每个数据点可以属于隶属函数范围在 0 和 1 之间的所有组。 0 是数据点距离聚类中心可能最远的点， 1 是数据点距离聚类中心最近的点**。**
使用案例:
当您需要进行图像分割或者您的目标是分割卫星图像中的水域、植被和岩石区时，模糊聚类非常有用。这在不能先验地确定聚类数的情况下是有用的。在这种情况下，可以合并具有弱边界的聚类。与 K-means 相比，模糊聚类在计算上是昂贵的，因为对于每个点，is 计算它属于每个聚类的概率。
基于密度的聚类
基于密度的聚类通过将高密度的区域分组并将它们与低密度的区域分开来工作。最著名的基于密度的聚类算法是 DBSCAN 算法(应用噪声的基于密度的空间聚类)。
使用以下两个参数计算密度
-
EPS: 这定义了数据点周围的邻域，即如果两点之间的距离小于或等于 EPS，则称它们为邻居
-
MinPts :定义形成一个邻域的数据点的最小数量。数据集的大小和 MinPts 的值成正比。
DBSCAN 算法访问每个点，并且如果它包含 eps 内的 MinPts，则聚类形成开始。任何其他点都被定义为噪声。这一过程一直持续到形成一个密度连通的集群，然后从一个新的点重新开始。
左图:传统聚类；右图:DBSCAN 聚类(DBSCAN 允许点采用任何形状或维度来形成聚类)
使用案例:
DBSCAN 主要用于平面空间的聚类。如果将它用于绘制自然灾害的影响图或绘制城市中气象站的位置图，可以获得良好的结果。当数据由非离散点组成时，也可以使用这种方法，这有利于处理异常值。推荐引擎/系统利用 DBSCAN 向他们的客户推荐产品/节目。
基于模型的聚类:
这种分类方法使用特定的分类模型，并尝试优化数据和模型之间的匹配。在基于模型的聚类方法中，数据被视为来自概率分布的混合，每个概率分布代表一个不同的聚类。换句话说，在基于模型的聚类中，假设数据是由混合概率分布生成的，其中每个分量代表一个不同的聚类。每个分量(即，集群)由正态分布或高斯分布建模。
期望最大化是一种众所周知的基于模型的聚类算法。当数据符合模型时，特定的聚类算法被认为工作良好。
使用案例:
这在聚类具有弱边界并且数据点在聚类之间具有混合成员的情况下非常有用。它在聚类协方差方面也更加灵活。在这种情况下，集群分配更加灵活，集群可以根据分布采用任何形状。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
