本文为《通过深度学习了解建筑年代和风格》论文复现的第三部分——获取阿姆斯特丹高质量街景图像的上篇，主要讲了如何获取利用谷歌街景地图自动化获取用于深度学习的阿姆斯特丹的高质量街景图像，此数据集将用于进行建筑年代的模型训练。

我们会从上文 Part2.下载和预处理建筑足迹数据集 获取到的阿姆斯特丹的 163210 条建筑足迹数据开始，获取用于下载街景图像的谷歌街景网址 128876 条 url，在下一篇文章我们会通过获取的 url 通过 selenium 进行街景图像的采集并分享我下载的完整的街景图像。

---

街景图像就是把建筑环境虚拟化展现出来，已经成为了现实世界体验和感知的一个替代品。

阅读前必看知识点

在阅读本文前需要了解的知识点，大部分都能在菜鸟教程找到，你也可以去相应的官方网站查找更多信息、

• Python 的基本知识，函数isinstance的使用
• 向量的点积（dot product）的概念
• numpy 的基本知识，dot（向量点积）、norm（范数）、arccos（反余弦）、degrees（弧度转换为度数）等部分函数的使用
• pandas 的基本知识
• geopandas 包如何处理 gdb、gpkg 等文件地理数据库
• Geodataframe 和 GeoSeries 的属性查看，切片和索引，apply 函数的使用，
• shapely 包中的几何对象：Point、MultiPoint、Polygon、MultiPolygon
• matplotlib 基本绘图方法
• Arcpy 的基本知识，要素属性表的查询和更新游标，Arcpy 中的几何对象 MultiPoint
• OSMnx包的介绍，如何下载 openstreet map 上的矢量数据

一、从街景图像的获取开始思考

1 方法一，超额收费：通过谷歌街景 API 获取街景图像

论文作者使用谷歌 API 去获取的街景图像，这是收费的，大约 7 美元一千张照片。如果我们去谷歌的开发平台绑定信用卡，可以获得每月 200 美元的额度，但是也只能能获取大约 3 万张照片，而我们的项目预估仅仅阿姆斯特丹一个城市就需要 7 万张照片。无奈我“囊中羞涩”，只能去使用一种免费的方法——selenium库实现浏览器自动化截图，去获取建筑立面的图像。

如果你恰好有经费，更推荐使用此方法下载街景图片，速度会比方法二快很多，它的代码主要使用streetview库获取pano_id变量后，使用get_streetview()函数下载：

from streetview import get_streetview
import os

# 下载街景图片
image = get_streetview(
    pano_id=pano_id, # 使用streetview包查询到的，后文会说
    api_key=GOOGLE_MAPS_API_KEY, # 在https://console.cloud.google.com/获取你的api
)

image.save("image.jpg", "jpeg")

2 方法二，完全免费：通过 selenium 实现批量街景图像的采集

Selenium 是一个强大的工具，用于控制 web 浏览器通过程序进行自动化操作。它的 Python 版本（Selenium WebDriver）允许你使用 Python 编写程序来做例如打开网页、填充表单、点击按钮等常见浏览器操作。

对于数据采集，使用 selenium 意味着只要能用浏览器打开的网页，使用 selenium 都能抓取到数据，即使涉及到了验证码都能处理。而且，你可以使用 selenium 去大麦抢演唱会票（ClassmateLin / dm-ticket）。我们这次用到的就是它的截屏功能：

为了获取不同建筑的照片，我们通过更改街景网页地址中的相关参数，来改变位置、转换角度等功能，就像你自己拿着一个摄像机 📸，最终针对每一个建筑调整到合适的角度，实现一个正对建筑立面的效果。

所以我们先分析网页地址中的相关参数，网页为https://www.google.com/maps/@52.3663507,4.8861953,3a,75y,42.63h,92.9t/data=!3m6!1e1!3m4!1sZwqc-JnbpT03nRj8Ublqjw!2e0!7i16384!8i8192?entry=ttu，其中我们要关注的参数如下表：

	Parameter	Description
1	@52.3663507,4.8861953	纬度和经度坐标的位置
2	3a	表示街景模式
3	75y	视野（缩放级别）
4	42.63h	相机指向的方向或方向（以度为单位），0-360°，正北方为起点。
5	92.9t	相机倾斜或角度（以度为单位）
9	1sZwqc-JnbpT03nRj8Ublqjw	唯一的全景 ID

• 1 纬度和经度
• 4 相机指向的方向或方向（以度为单位）
• 9 唯一的全景 ID

以上三个参数是我们着重需要调整的，你可以手动转动街景地图，或者是向前后移动，初步观察网址中相应参数随网页街景图变化的规律：

• 相机指向的方向在 0-360 度中变化，最北的为 0 度，角度随顺时针方向增加。
• 全景 ID：似乎没有规律。
• 纬度和经度：需要查阅资料。
• 视野（缩放级别）：数值最大为 90，值越大距离建筑物越远。

3 详解谷歌街景网页 URL 中的三个重要参数：

1）纬度和经度 lat, lng

查阅资料后发现，谷歌街景（Google Street View, 简称 GSV）使用的坐标系是 WGS 84（World Geodetic System 1984），因此，需要将之前提取的阿姆斯特丹的建筑矢量数据（要素Amsterdam_buildings）转为 WGS 84 坐标系。这很好办，我们使用 ArcGIS Pro 中的投影工具进行转换

2）相机指向的方向或方向 heading θ

既然是需要正面的建筑照片，我们肯定是在道路上拍摄的，网页中的角度也是以道路上的点拍摄的，并且这个点最好要满足距离建筑足够近，这样才能获取比较清晰的建筑立面照片，这个最佳的街景拍摄点我们称之为Point S，也就是 Street Point。此时Point S的方位角（以北为起点，顺时针旋转的角度）叫做 θ，就是网页中需要填入的角度。那如何找到此点，论文提出了一种方法找到此点。并计算 θ：

1. 找到Point S：为了使拍摄的照片包含建筑的大部分里面信息，我们的拍摄方向一定从街景拍摄点Point S(xs, ys)的朝着建筑物最靠近街道的底面边的中点拍摄的，假设存在多个建筑底面边的中点，谁距离它最近的道路的距离最短，那一定就是我们找的建筑上的点，我们称之为Point C(xc, yc)。

2. 计算 θ：我们上一步获取了两个点的坐标，然后通过向量之间的点积（dot product）就可以求出*Point S(xs, ys)*和 Point C(xc, yc)的夹角，具体是如**方程式(1)**所示：

​ 其中矢量 Vn 是北向矢量，矢量 Vsc 是从 S 点到 C 点的矢量。

北向量 Vn 被定义为从点 S(xs, ys) 指向正北方向的向量。正北方向通常与垂直轴（即 y 轴）对齐。因此，为了得到一个指向正北方向的单位向量，我们可以将 y 坐标增加 1，而 x 坐标保持不变。

更详细点说：

因为两个向量的点积和它们的模（magnitude）的乘积之间的比和这两个向量之间的 cosine 值存在关系，例如给定两个向量 A 和 B：

它们之间的点积定义为：

其中，θ 是向量 (A) 和 (B) 之间的夹角。

从上述公式中，我们可以得到:

这就是为什么点积和两个向量的模的乘积之间的比值可以得到这两个向量之间的 cosine 值。这个关系在计算向量之间的角度时非常有用，因为我们可以使用 arccos 函数来从 cosine 值得到实际的角度。

我们用 Python 中的numpy，按照上述**公式(1)**实现计算向量 A 和向量 B 的角度 θ ：

import numpy as np

def calculate_angle(xs, ys, xc, yc):
    """用于计算两个向量之间的夹角。这两个向量是从点S到点C的向量Vsc和北向量Vn。返回消毒"""
    # 定义北向量 Vn
    Vn = np.array([xs, ys + 1])

    # 定义从点S到点C的向量 Vsc
    Vsc = np.array([xc - xs, yc - ys])

    # 计算Vn和Vsc的点积
    dot_product = np.dot(Vn, Vsc)

    # 计算Vn和Vsc的模（或长度）
    norm_Vn = np.linalg.norm(Vn)
    norm_Vsc = np.linalg.norm(Vsc)

    # 计算夹角的余弦值
    cos_theta = dot_product / (norm_Vn * norm_Vsc)

    # 根据余弦值计算夹角theta（以度为单位）
    if (xc - xs) > 0:
        theta = np.degrees(np.arccos(cos_theta))
    else:
        theta = 360 - np.degrees(np.arccos(cos_theta))

    return round(theta, 2) # 保留小数点两位

# 测试函数
xs, ys = 2, 3  # 点S的示例坐标
xc, yc = 4, 1  # 点C的示例坐标
calculate_angle(xs, ys, xc, yc)

out:

108.43

😵😵😵 线代的知识看不懂也没关系，我们可以直接使用 ArcGIS Pro 中的邻近分析来找到最近点并且计算好角度 θ，后文会说。

3）全景 ID

搜索之后发现，全景 ID 是一个唯一值，需要使用streetview获取它，首先你要确保能科学上网，然后需要安装streetview：

pip install streetview==0.0.6

from streetview import search_panoramas

panos = search_panoramas(lat=41.8982208, lon=12.4764804)
first = panos[0]

print(first)

pano_id='_R1mwpMkiqa2p0zp48EBJg' lat=41.89820676786453 lon=12.47644220919742 heading=0.8815613985061646 pitch=89.001953125 roll=0.1744659692049026 date='2019-08'

其中 pano_id、lat 和 lon 即全景 id，纬度和经度我们都需要进行储存。

4 获取街景采集点 Point S 和 heading θ 的思路

1. 定义建筑物的底面边界：读取阿姆斯特丹的建筑足迹数据，获取每个建筑的各个边。

2. 计算建筑物各边的中心点：遍历建筑物的所有边，计算每条边的中心点。如果建筑物的每条边由点 A 和点 B 定义，那么中心点 C 的坐标是 ((A.x + B.x) / 2, (A.y + B.y) / 2)。

3. 获取道路矢量数据：下载 Open Street Map 的道路矢量数据，我们可以通过OSMnx包去下载并进行简化。

4. 找到最近的点 Point C：对于建筑物的每个边的中心点，计算它到道路的每个段的最近距离。

   • 计算点到线段的垂直距离，可以通过向量数学或使用一些专用的几何算法来完成。也可以使用 Shapely 库计算最短距离。
   • 对于每个中心点，您将遍历道路上的所有线段，找到点到线段的最近距离。保存这个距离和对应的线段。

5. 比较距离：一旦您有了从各个中心点到道路的距离，您就需要找出哪个距离最短。通过比较所有计算出的距离来完成这一点。

6. 确定最短距离的坐标 Point S：找出最短距离后，返回对应的中心点的坐标以及该点到最近道路边界的距离。

7. 计算 heading θ：利用两个点的坐标计算角度

很明显，Point S在道路上，并且是车行的公共道路。并且，Point S必须要是距离Point C较近的点。在阿姆斯特丹，OSM(openstreet map)的数据非常全面，所以我们通过 Python 的OSMnx的库来下载和处理 OSM 的道路数据。

5 获取阿姆斯特丹的道路矢量数据

OSMnx是一个 Python 库，用于从 OpenStreetMap 下载、建模、分析和可视化街道网络和其他地理空间功能。您可以下载和建模步行、驾驶或骑自行车的网络，只需一行代码，然后轻松地分析和可视化它们。您可以轻松地处理城市设施/兴趣点、建筑物占地面积、公交站点、高程数据、街道方向、速度/行驶时间和路线。

# 下载道路数据
G = ox.graph_from_place("Amsterdam, North Holland, Netherlands", network_type="drive")

# 简化道路数据
G2 = ox.consolidate_intersections(G_proj, rebuild_graph=True, tolerance=15, dead_ends=False)

# 绘制道路
fig, ax = ox.plot_graph(G2, node_size=0) # 不绘制node节点

• ‘drive’ - 获得可驾驶的公共街道
• ‘drive_service’ - 获得可驾驶的公共街道，包括服务道路
• ‘walk’ - 获取行人可以使用的所有街道和路径（这种网络类型忽略单向方向性）
• ‘bike’ - 获取骑自行车者可以使用的所有街道和路径
• ‘all’ - 下载所有（非私有）OSM 街道和路径
• ‘all_private’ - 下载所有 OSM 街道和路径，包括私人访问的

我们直接导出并用 Arcgis Pro 打开看看，如果你想用 geopandas 继续处理，可以转为 geodataframe 对象然后处理：

ox.save_graph_geopackage(G2, filepath="../data/UNZIP/osm/Amsterdam_road.gpkg")

放大之后我们发现部分道路（下左）还是没有被简化，但是对比谷歌街景地图的轨迹（下右）之后，复杂度差不多，因为最终的经纬度要根据streetview包中返回的经纬度决定，这一步的道路用于找到大致的经纬度，所以也不需要简化了。

二、获取建筑物的各边中点

简化建筑物指的是在保持建筑物基本形状和大小不变的前提下简化建筑物面的边界或轮廓。

在本项目中，未简化的建的边会过多（如下图 single_building），而过短的边会增加后续寻找建筑的街景点的计算量，从而影响计算效率。

在深度学习提取建筑轮廓的时候，我们也需要对识别的建筑物进行简化，也称为 Polygonization，多边形化。目的是为了让建筑物的各边更规则：

在 ArcGIS Pro 中有一个工具叫规则化建筑物覆盖区 (3D Analyst)也实现了同样的功能。

---

接下来我会演示两种方法来获取建筑物的各边中点，分别是方法 1：用 geopandas 和shapely处理建筑并获取中心点，方法 2：使用 ArcGIS Pro 的游标获取中心点，会使用到 Arcpy。

2.1 方法一：用 geopandas 和 shapely 处理建筑并获取中心点

1) 简化建筑物

我们先用 geopandas 读取建筑足迹数据，注意需要安装高版本的 geopandas 才能读取文件地理数据库(gdb)。

import geopandas as gpd
print(gpd.__version__) # 0.12.0 低版本的不能读取gdb数据库

gdb = "../../5-ArcgisPro工程.建筑风格和年代深度学习.gdb"
gdf = gpd.read_file(gdb, layer='Amsterdam_buildings_Project', rows=5000)  # layer就是数据库中的一个要素类 rows代表只读取前多少条数据

读取的数据的类型为 GeoDataframe，简称为 gdf。

可以使用 shapely 的object.simplify(tolerance, preserve_topology)方法对单个集合体进行简化。

object.simplify(tolerance, preserve_topology) 方法来自 Shapely 库，用于简化几何对象。该方法基于 Douglas-Peucker 算法返回输入几何的简化版本。参数包括：geometry（几何对象或类似数组），tolerance（float 或类似数组，表示允许的最大几何位移），以及一个可选的 preserve_topology（默认为 True，用于保留几何拓扑结构），详见[官方文档](https://shapely.readthedocs.io/en/stable/reference/shapely.simplify.html#:~:text=shapely.simplify ,The maximum allowed geometry displacement)。

简化前的建筑

# 选一个样本
sample = 3255

# 简化前
gdf.loc[sample, "geometry"]

gdf.loc[sample, "geometry"].geom_type, type(gdf.loc[sample, "geometry"])

out：

('MultiPolygon', shapely.geometry.multipolygon.MultiPolygon)

可以看到简化是一个 MultiPolygon，存在一些不需要的边。

简化后的建筑

我们对这个多边形进行简化，并查看结果：

# 简化
tolerance_m = 1 # 容差通常以地图单位为单位（例如，米、英尺等），这取决于您的地图或空间数据的坐标系统。此处是米。
gdf.loc[sample, "geometry"].simplify(tolerance=tolerance_m, preserve_topology=True)

不错 👍！短边被去除了！！

对所有建筑物都进行简化

那么如何对所有的建筑物几何体都进行简化？

方式一是对所有的建筑物逐个应用 shapely 的object.simplify(tolerance, preserve_topology)方法，即使用 apply 的函数对读取的 gdf，

然后我们使用geopandas.GeoSeries.simplify对gdf.geometry进行简化，而不需要写匿名 lambda 函数：

geopandas.GeoSeries.simplify方法用于返回一个包含每个几何体简化表示的GeoSeries。这个方法基于 Douglas-Peucker 算法，该算法递归地将原始线分割成较小的部分，并通过直线连接这些部分的端点。然后，它会移除所有到直线距离小于tolerance的点。该方法不会移动任何点，并且总是保留原始线或多边形的端点，详见官方文档。

• 参数:
  • tolerance (float): 简化几何体的所有部分将与原始几何体的距离不超过tolerance。它与GeoSeries的坐标参考系统的单位相同。例如，在投影的坐标参考系统中，如果单位是米，那么tolerance=100意味着现实中的 100 米的距离。
  • preserve_topology (bool, 默认值为 True): 如果为False，则使用更快的算法，但可能会产生自相交或其他无效的几何体。
• 注意:
  • 如果不保留拓扑结构，简化可能会导致无效的几何对象，而且简化可能对坐标的顺序敏感：仅在坐标顺序上不同的两个几何体可能会被不同地简化.

我们对通过 gdf.geometry 将其转化为 GeoSeries 对象，然后使用 simplify 方法：

s = gdf.geometry.simplify(tolerance=tolerance_m, preserve_topology=False)

2）使用 Shapely 获取建筑各边的中心点

要获取 GeoPandas 集合体（例如 GeoSeries 或 GeoDataFrame）中每个多边形的外边界上所有中点，你可以使用 Shapely 库的几何对象方法和属性。每个多边形的外边界可以通过polygon.exterior获得，该属性返回一个LinearRing对象。然后，你可以使用LinearRing对象的coords属性来访问边界上的坐标点，并计算相邻点之间的中点。

详解 polygon.exterior.coords

在 Python 的 Shapely 库中，polygon.exterior.coords 不是一个函数，而是一个属性。polygon.exterior 返回一个 LinearRing 对象，该对象代表多边形的外部轮廓，而 .coords 属性则提供该外部轮廓的坐标点。

使用len(polygon.exterior.coords)找到该对象的长度并可以像列表一样索引对象。 例如，要获取第一个顶点，使用polygon.interiors[0].coords。❗ 注意，第一个点和最后一个点是相同的。 所以要包含没有重复点的顶点的列表，需要使用polygon.interiors[0].coords[:-1]。

获取多边形外轮廓的坐标列表

# 转换多边形的边界为线性环
linear_ring = gdf_simplify.loc[sample, "geometry"].exterior

# 将坐标转换为一个列表，每个元素都是一个坐标对 (x, y)
coords_list = list(linear_ring.coords)
print(len(coords_list), coords_list)

out:

5 [(632140.1681000004, 5803521.2783), (632140.5662000002, 5803509.002), (632125.7775999997, 5803508.522), (632125.3492999999, 5803520.7722), (632140.1681000004, 5803521.2783)]

可以看到第一个点和最后一个点是相同的。

获取每相邻两个坐标点的中点坐标

因为我们是在投影坐标系下，获取中点只需要分别计算经度和维度的两个点的坐标的平均值：

from shapely.geometry import Point
# 初始化列表，用于存储每条边的端点
mid_points = []

# 遍历坐标列表中的点，每两个点组成线段的两个端点
for i in range(len(coords_list) - 1):  # 减1，因为我们回到了起始点，所以最后一个点不需要形成新的线段
    # 获取线段的端点
    point1 = coords_list[i]
    point2 = coords_list[i + 1]

    # 计算中点的坐标
    mid_x = (point1[0] + point2[0]) / 2.0
    mid_y = (point1[1] + point2[1]) / 2.0

    # 返回中点
    mid_point = Point(mid_x, mid_y)

    mid_points.append(mid_point)


print(len(mid_points), mid_points)

out：

4 [<POINT (632140.367 5803515.14)>, <POINT (632133.172 5803508.762)>, <POINT (632125.563 5803514.647)>, <POINT (632132.759 5803521.025)>]

我们从一个四边形中获取到了他的四个边的中心点。

将中心点转为 GeoDataFrame 并绘制出来

# 将中心点构造成dataframe
points = gpd.GeoDataFrame(geometry=mid_points)

# 为 GeoDataFrame 设置坐标参考系统（CRS）
points.crs =  gdf_simplify.crs
points

绘制：

# 绘制原始多边形和简化后的多边形
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(7, 7))

# 绘制原始多边形。我们通过将 'loc' 结果转换为一个新的 GeoDataFrame 来确保 'plot' 方法的可用性。
# 'loc' 通常返回一个 Series 或单个值，但我们需要一个 GeoDataFrame 来使用 'plot' 方法。
gdf_simplify.loc[[sample]].plot(ax=ax, facecolor='gray') #注意，这里用双括号来确保返回一个 DataFrame

# 绘制中点
points.plot(ax=ax, color='red', markersize=100)

利用 apply 批量处理所有建筑物

我们为 GeoDataFrame gdf_simplify 中的每个多边形和多多边形计算边的中点，并将这些中点存储在新的 midpoints 列中。每个 midpoints 值都是一个 MultiPoint 对象，包含一个多边形或多多边形的所有边的中点。

import geopandas as gpd
from shapely.geometry import Point, MultiPoint, Polygon, MultiPolygon

def calculate_midpoints(geometry):
    """
    计算Polygon的每条边的中点或MultiPolygon中每个Polygon的中点。
    接受一个几何对象（geometry）作为输入，并返回一个包含所有中点的 MultiPoint 对象。
    """
    midpoints = []

    def midpoint_for_polygon(polygon):
        # 为单个多边形计算中点
        coords = polygon.exterior.coords
        for i in range(len(coords)-1):
            mid_x = (coords[i][0] + coords[i + 1][0]) / 2.0
            mid_y = (coords[i][1] + coords[i + 1][1]) / 2.0
            midpoints.append(Point(mid_x, mid_y))

    if isinstance(geometry, Polygon):
        midpoint_for_polygon(geometry)
    elif isinstance(geometry, MultiPolygon):
        # 对于MultiPolygons，我们需要遍历每个多边形。
        for poly in geometry.geoms:  # 在这里，我们使用.geoms属性进行迭代。geoms属性返回它包含的Polygon实例。
            midpoint_for_polygon(poly)

    return MultiPoint(midpoints)

# 应用函数并创建一个包含中点的新GeoDataFrame
gdf_simplify['midpoints'] = gdf_simplify.geometry.apply(calculate_midpoints)

❗ 如何处理两列带有 geometry 属性的 GeoDataFrame：

我们的gdf_simplify现在有两列带有 geometry 属性，他们可以同时存在，而且后续会继续使用到这两列（上图左侧建筑面，上图右侧建筑各边中点）数据。需要说明的是，gdf_simplify.geometry的调用的属性返回的是左侧的那一列，我们可以使用gdf_simplify.set_geometry('midpoints')来改变gdf_simplify.geometry的属性，从而让他返回右侧（midpoints）列。此时要想调用左侧的 geometry，需要使用gdf_simplify.loc[:, "geometry"]。

---

我们此时的列很多，不需要这么多列，我们只保留[“identificatie”, “bouwjaar”, “midpoints”, “geometry”]这四列，同时使用gdf_simplify.set_geometry('midpoints')来改变gdf_simplify.geometry的属性：

gdf2 = gdf_simplify.copy()
gdf2 = gdf2.set_geometry('midpoints')
gdf2 = gdf2[["identificatie", "bouwjaar", "midpoints", "geometry"]]

2.2 方法二：用 ArcGIS Pro 和 ArcPy 处理建筑并获取中心点

1）简化建筑物

为了不受建筑短边的影像，我们需要对建筑进行简化，使用 ArcGIS Pro 的简化建筑物工具能达到 shapely 的效果，而且还能输入障碍图层，障碍图层可以是挡住视线的其他建筑物，我们先不考虑障碍图层。

简化前后对比：

2） 获取建筑各边中心点

接下来我们在 ArcGIS 软件中的 notebook 中进行获取建筑中心点的操作：

首先我选择了单个建筑input_polygon_feature_class，通过查询游标获取此建筑多边形的各个边的中点(mid_x, mid_y)，同时定义一个Arcpy的多点几何对象(Multipoint)，最后确定坐标系为 UTM 的投影坐标系（WGS 84 / UTM 北 N ZONE 31），EPSG 编号是 32631，我们通过编号定义一个空间参考sr。

import arcpy

# 输入的已简化多边形要素类
input_polygon_feature_class = "single_buildin_SimplifyPolyg2"

# 定义SpatialReference对象，这里以WGS 1984坐标系为例
sr = arcpy.SpatialReference(32631)

# 使用SearchCursor遍历多边形
with arcpy.da.SearchCursor(input_polygon_feature_class, ["identificatie", "SHAPE@"]) as cursor:
    for row in cursor:
        polygon_id = row[0] # 建筑的id
        polygon = row[1]

        # 创建一个Array
        point_array = arcpy.Array()

        # 遍历多边形的每一段，计算并存储中心点
        # 遍历poly中所有的多边形，part为每个单独的多边形
        for part in polygon:
            # 遍历每个边
            for i in range(len(part) - 1):
                start_point = part[i]
                end_point = part[i + 1]
                mid_x = (start_point.X + end_point.X) / 2
                mid_y = (start_point.Y + end_point.Y) / 2
                print(mid_x, mid_y)

                # 创建arcpy的点
                midpoint = arcpy.Point(mid_x, mid_y)

                # 添加到Array
                point_array.add(midpoint)

            # 使用Array创建Multipoint对象
            multipoint = arcpy.Multipoint(point_array, sr)

OUT：

627016.5058999998 5805486.551
626992.1195999999 5805457.27365
626960.1142500001 5805483.0801
626984.50055 5805512.35745

有关 ArcPy 中的几何对象和游标的使用可以查看六、处理几何数据【ArcGIS Python 系列】进行更多了解。

# notebook内直接查看multipoint的形状
multipoint

multipoint 是由四个点组成的对象。

---

接下来将其保存到数据库中，以便后续操作，几何对象可以作为要素直接进行运算，但是邻近分析需要对几何对象新增字段和更新字段，所以得先保存为要素。

arcpy.CopyFeatures_management(multipoint, "multipoint")

三、找到街景采集点和对应的建筑物中点并 heading 角度

在前文已经讲解了如何获取街景采集点的位置和 heading 角度，我们回顾一下：

我们现在得到了路网，为了实现街景图形和建筑一一对应，论文用了一种巧妙的方法获取建筑立面所对应的 GSV 的位置：

图 3. 从 GSV 检索建筑立面的过程。在步骤 1 中，建筑物外墙的中点（红点）投影到最近的街道，该点用作请求 GSV 的位置，在步骤二中，计算向量北与从请求点到外墙中点的向量之间的角度并将其输入 Google 地图 API 作为相机角度。

3.1 使用 geopandas 找到街景点（方法 1）

建议用方法一，因为速度更快。如果你想学如果使用 ArcGIS Python 也就是 Arcpy 如何处理空间数据，也推荐看看第二种方法。

我们要利用 shapely 计算距离，但是我们的 road 要素是整个地区的线要素，我们需要裁剪到到我们的建筑物的外围，方便计算。我们先对读取 s 行问获取的阿姆斯特丹的道路，进行简单修复，然后对对建筑物做缓冲区，用于提取建筑物周围的道路，减少计算量。

1）读取阿姆斯特丹矢量道路数据

# 读取道路
road = gpd.read_file(r"..\5-ArcgisPro工程\Amsterdam_road.gpkg", layer="edges")

# 进行简单处理
road = road.to_crs(gdf2.crs)[["geometry"]].dropna()

# 修复无效的几何形状
road['geometry'] = road['geometry'].apply(lambda geom: geom if geom.is_valid else geom.buffer(0))

road

out：

2）对建筑做缓冲区

sample = random.randint(0, len(gdf2))
gdf2.loc[sample, "geometry"]

# 获取建筑物外轮廓并扩展一定距离
gdf2.loc[sample, "geometry"].buffer(50)

3）裁剪道路数据

#裁剪道路
road_clip = road.geometry.intersection(gdf2.loc[sample, "geometry"].buffer(50))

# 先扔掉空值的行
road_clip.dropna()

# 过滤掉空的几何形状然后将合并分散的线
road_nearby = GeometryCollection([geom for geom in road_clip if not geom.is_empty])
road_nearby

4）使用 shapely 的 nearest_point 找出最近的两个点

from shapely.geometry import Point
from shapely.ops import nearest_points

# 初始化变量
shortest_distance = float('inf')  # 初始化为无穷大
nearest_center_point = None  # 最近的点将存储在这里

center_points = gdf2.loc[sample, "midpoints"].geoms # 从MultiPoint中获取所有的点

for point in center_points:
    # 求当前点到道路的最近点（道路上的点）
    nearest_road_point = nearest_points(point, road_line)[1]

    # 计算距离
    distance = point.distance(nearest_road_point)

    # 如果当前距离比记录的距离更短，则更新最短距离和最近的中心点
    if distance < shortest_distance:
        shortest_distance = distance
        nearest_center_point = point

# 输出结果
print(f"最近的中心点是: {nearest_center_point}")
print(f"最短距离是: {shortest_distance}")
print(f'最近的道路点是: {nearest_road_point}')

out:

最近的中心点是: POINT (632175.4186 5806202.752800001)
最短距离是: 16.51552203497222
最近的道路点是: POINT (632184.9900678364 5806216.517023285)

5）使用向量相乘的原理计算两个点间的角度 heading θ

import numpy as np


def calculate_angle(xs, ys, xc, yc):
    # 定义北向量Vn
    Vn = np.array([xs, ys + 1])

    # 定义从点S到点C的向量Vsc
    Vsc = np.array([xc - xs, yc - ys])

    # 计算Vn和Vsc的点积
    dot_product = np.dot(Vn, Vsc)

    # 计算Vn和Vsc的大小（范数）
    norm_Vn = np.linalg.norm(Vn)
    norm_Vsc = np.linalg.norm(Vsc)

    # 计算角度的余弦
    cos_theta = dot_product / (norm_Vn * norm_Vsc)

    # 以度为单位计算角theta的角度
    if (xc - xs) > 0:
        theta = np.degrees(np.arccos(cos_theta))  # 如果角在x轴的正半轴
    else:
        theta = 360 - np.degrees(np.arccos(cos_theta))  # 如果角在x轴的负半轴

    return theta

# 测试一下
xs, ys = points.x, points.y
xc, yc =  pointc.x, pointc.y
print(xs, ys, xc, yc)
calculate_angle(xs, ys, xc, yc)

OUT：

(632184.9900678364 5806216.517023285 632175.4186 5806202.752800001)
208.60033988006262

6）我们绘制出相应的点来验证是否计算正确

# 和实际结果验证一下绘制结果
import matplotlib.pyplot as plt

# 创建一个图来绘制
fig, ax = plt.subplots()

# 绘制四个中心点
gdf2.loc[[sample], "midpoints"].plot(ax=ax)

# 绘制建筑物
gdf2.loc[[sample], "geometry"].plot(ax=ax, facecolor='gray', alpha=0.5)

# 绘制箭头
ax.quiver(xs, ys, xc - xs, yc - ys, angles='xy', scale_units='xy', scale=1, color='red')
# 绘制中心点
ax.plot(pointc.x, pointc.y, marker='o', color='red', markersize=5)
# 绘制道路上的点
ax.plot(points.x, points.y, marker='o', color='green', markersize=5)


# 设置额外的图形属性，如标题和标签（如果需要）
ax.set_title(f'id:{gdf2.loc[sample, "identificatie"]}  Location')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')

# 显示图形
plt.show()

角度 heading θ 是最容易出错的地方，我验证了多次，并且在后续或者街景地图之后也进行了验证，确保不会出错，否则后续获取街景图的工作就废了。

7）整合并简化用 geopandas 寻找街景点的代码

上述只用了一个建筑物去测试，我们要使用 apply 函数对整个建筑数据集进行处理，同时使用concurrent.futures多线程优化代码处理速度，并且创建calculate_midpoints、calculate_angle、read_road、process_geometry_and_calculate_angle、parallel_processing、save_results几个函数。

《使用 geopandas 寻找街景点》的完整代码文件获取方式：

关注本公众号renhailab，选择点赞、在看或者转发本文之后，私信20231027获取代码。码字不易、多多点赞。

3.2 使用使用 ArcGIS Pro 进行邻近分析，找到街景点和角度（方法 2）

1）投影到同一个 UTM 坐标

我们要找到Amsterdam_buildings靠近其最近道路nodes的立面的重点Point C，然后确定Point S。

首先我们统一坐标系到 WGS1984 地理坐标系和 WGS_1984_UTM_Zone_31N 投影坐标系：

2）找到建筑立面的中心点并进行邻近分析

进行建筑简化之后生成中心点，最后使用 ArcGISPro 的邻近分析 (分析)工具找到建筑物外墙的中点（Point C）投影到最近的街道，该点用作请求 GSV 的位置(Point S)。

之后进行邻近分析，从每个点找到距离周围道路最近的点 Point C，以及 Point S，heading 值：

# 邻近分析 (分析)
in_features = "multipoint"
near_features = "main.edges"
search_radius = "25 Meters"
arcpy.analysis.Near(in_features,
                    near_features,
                    "25 Meters", "LOCATION", "ANGLE", "PLANAR", "NEAR_FID NEAR_FID;NEAR_DIST NEAR_DIST;NEAR_X NEAR_X;NEAR_Y NEAR_Y;NEAR_ANGLE NEAR_ANGLE")

在地图中查看：

同时可以看到 multipoint 对象的属性表也被更新了：

3）查询 multipoint 中的角度并进行角度转换

NEAR_DIST是街景点与最相邻建筑边终点的距离，NEAR_X和NEAR_Y则为要找的街景点。我们可以用XY 表转点将两个字段转为一个点：

可以看到蓝点就是我们要找到街景采集点 Point S。

---

角度NEAR_ANGLE有特殊的规定：生成的角度是输入要素（建筑物上的点）对于邻近要素（街道上的点）的角度，并且转换前的角度的表示方式是：在方法参数中使用平面方法时，角度在 -180° 到 180° 的范围内，0° 代表东，90° 代表北，180°（或 -180°）代表西，-90° 代表南。
转换后，因为要获取街道街景，要以（街道上的点）为原点，朝向建筑物上的点的角度，并且角度表示为：角度范围在 0-360 度，0° 代表北，90° 代表东，180° 代表南，270° 代表西。，即：

def transform_angle(original_angle):
    """
    将角度从一个坐标系转换为另一个，并更改方向表示。

    :param original_angle: 初始的角度（基于东为0°的系统）
    :return: 转换后的角度（基于北为0°的系统）
    """
    # 从建筑物到街道的角度需要将角度旋转180度以“反转”方向
    reversed_angle = original_angle + 180

    # 规范化角度在0到360之间
    if reversed_angle >= 360:
        reversed_angle -= 360
    elif reversed_angle < 0:
        reversed_angle += 360

    # 现在，我们需要将“东为0度”转变为“北为0度”，这需要一个90度的逆时针旋转
    north_based_angle = reversed_angle + 90

    if north_based_angle >= 360:
        north_based_angle -= 360

    return north_based_angle

print(transform_angle(-178.9227630675157))

接下来把得到的数据保存下来，首先进行查询：

# 定义要查询的字段
field_name_list = ["NEAR_DIST","NEAR_X","NEAR_Y", "NEAR_ANGLE"]

# 使用SearchCursor迭代访问每个记录
with arcpy.da.SearchCursor(in_features, field_name_list) as cursor:
    for row in cursor:
        # 获取NEAR_ANGLE字段的值并添加到列表中
        print([i for i in row])
        angle = row[3]
        print("转换前的角度：", angle)
        angle2 = transform_angle(angle)
        print("转换后的角度",angle2)


# 打印所有的NEAR_ANGLE值
print("当前建筑的id", polygon_id)

OUT:

[15.729911607404052, 626640.935368494, 5799565.214474143, -178.9227630675157]
转换后的角度 91.07723693248431
当前建筑的id 0363100012061237

然后进行保存：

# 追加保存到json文件，然后通过streetview进行构建url的操作
import pandas as pd

# 从字典创建 DataFrame
df = pd.DataFrame({
    "polygon_id": [polygon_id],
    "NEAR_DIST": [row[0]],
    "lat": [row[2]],
    "lng": [row[1]],
    "heading": [angle]
})

# 追加数据到 CSV
df.to_csv("./temp.csv", mode='a', header=False, index=False)

3）简化和整合代码

我们进行代码整合，对全市的建筑进行处理，完整代码见：。

本来想进行多线程处理，但是在 ArcGIS 的采用多线程会遇到数据锁，不能同时被占用，真实头疼。为了加快速度，我手动进行了"多线程"——选择了克隆项目然后运行多个处理不同数据子集的 Python 程序，我按照 5000 一组分为了 31 个子数据集，5000 个建筑处理大约为 40 分钟，时间很长，如果有更快的方法，我会在仓库更新。

《使用 Arcpy 寻找街景点》的完整代码文件获取方式：

关注本公众号renhailab，选择点赞、在看或者转发本文之后，发送20231027获取代码。码字不易、多多点赞。

---

四、构建谷歌街景图片 url

1）通过 streetview 获取经纬度、朝向

上文已经将谷歌街景网页中的三个重要参数获取了到了，接下来我们构建谷歌街景图片的网址(url)。

首先用 pandas 读取保存有经纬度、朝向和建筑 id 的表格数据（此时读取的测试数据，不代表最终结果）：

import pandas as pd
df = pd.read_csv(r'../5-ArcgisPro工程/temp3.csv', header=0, encoding='utf-8')
df.head()

我们通过streetview库获取最新的一个的全景 id（panos），然后和其他参数组合成最终的谷歌街景网页地址，用于后续 selenium 获取街景：

## 构建url
from streetview import search_panoramas
from datetime import datetime

print(df.loc[1,:])
lat = df.loc[1,'lat']
lng = df.loc[1,'lng']

heading = df.loc[1,'heading']
panos = search_panoramas(lat, lng)

# 获取当前日期和时间
current_date = datetime.now()

# 初始化最近日期和对应的Panorama对象
closest_date = None
closest_pano = None

# 遍历panos列表
for pano in panos:
    if pano.date:
        # 将日期字符串转换为datetime对象
        pano_date = datetime.strptime(pano.date, '%Y-%m')
        # 计算日期之间的差值
        delta = current_date - pano_date
        # 如果closest_date为None或者当前日期更接近于closest_date，则更新closest_date和closest_pano
        if closest_date is None or delta < closest_date:
            closest_date = delta
            closest_pano = pano

# 输出最近日期的
closest_date, closest_pano

OUT:

(datetime.timedelta(days=955, seconds=55312, microseconds=899666),
 Panorama(pano_id='a8USRXjhCtTGXrjwHS5HJA', lat=52.35513228204291, lon=4.992971208322429, heading=130.8642120361328, pitch=88.9742431640625, roll=2.590848445892334, date='2021-03'))

2）组合 url

我们获取到了日期date为 2021-03 的一张照片，pano_id=‘a8USRXjhCtTGXrjwHS5HJA’，注意此时角度heading不是我们要的角度，但是lat和lon是需要的，我们将他们为变量：

lat = closest_pano.lat
lng = closest_pano.lon
pano_id = closest_pano.pano_id

url = f"https://www.google.com/maps/@{lat},{lng},3a,60y,{heading}h,95t/data=!3m6!1e1!3m4!1s{pano_id}!2e0!7i16384!8i8192"
print(url)

OUT：

https://www.google.com/maps/@52.35513228204291,4.992971208322429,3a,60y,40.83320085243709h,95t/data=!3m6!1e1!3m4!1sa8USRXjhCtTGXrjwHS5HJA!2e0!7i16384!8i8192

我们打开此链接：

3）整合并简化代码

我们将上述代码合成一个整体，并简化代码。通过使用pandas的apply方法更高效地遍历 df 中的每一行。
通过使用列表推导式和min函数，可以更高效地找到日期最近的pano。

---

上一篇：Part2.下载和预处理建筑足迹数据集——《通过深度学习了解建筑年代和风格》

下一篇：Part3.获取高质量的阿姆斯特丹建筑立面图像（下）——《通过深度学习了解建筑年代和风格》

---

因为其他平台不能同步修改，论文解读文章将最先在我的博客发布，你可以点击阅读原文获得查看博客上的原文。

好了，码字不易，如果如果你觉得本系列文章有用，欢迎关注博客，点赞 👍 和收藏，也欢迎在评论区讨论，有任何问题都可以私信我。涉及到代码报错问题可以在 Github 的本仓库的 issue 页面提问，Gitee 仅用于同步：

• 如果你觉得本系列文章有用，欢迎关注博客，点赞 👍 和收藏，也欢迎在评论区讨论，也欢迎访问我的爱发电支持我，或者对此文章进行赞赏。

  

  其他平台账号：
  

  • 我的知乎
  • 我的博客
  • 我的 GITHUB
  • 我的 GITEE
  • RSS

写在最后

论文引用：

Maoran Sun, Fan Zhang, Fabio Duarte, Carlo Ratti,
Understanding architecture age and style through deep learning,
Cities,
Volume 128,
2022,
103787,
ISSN 0264-2751,
https://doi.org/10.1016/j.cities.2022.103787.
(https://www.sciencedirect.com/science/article/pii/S0264275122002268)