本文为《通过深度学习了解建筑年代和风格》论文复现的第六篇——对建筑年代深度学习模型的进行评价，我们首先会通过对测试数据集的预测来展示模型的预测能力，其中，我们会介绍对模型进行评估的几种方法，包括混淆矩阵、召回率 (Recall)、精确度 (Precision)、F1 分数 (F1 Score)，然后，我们会利用类激活映射（Class Activation Mapping，简称 CAM）查看模型关注哪些方面，最后从空间上观察建筑年代的预测结果在空间上的表现。

《通过深度学习了解建筑年代和风格》论文复现代码已上传到Github和Gitee，还在更新中，请访问对应的主页查看。

本文主要对应的是论文5. Results 部分，会复现以下几张图：

⬇️ 模型预测可视化结果

⬇️ 表 4 混淆矩阵（百分比）

⬇️ 图 10 CAM 去识别不同年代模型的关注点

1. 左侧小图是将 CAM 叠加在原始图像上。图像的红色区域主要覆盖一楼和二楼之间的窗户或门。
2. 右侧小图：根据 CAM 裁剪的图像显示了窗户的演变。早期的窗户通常框架较宽，装饰较多，而且较窄。最近的窗户样式以方形和水平形状为特点，框架更薄，装饰更少，深度更小。

⬇️ 图 7 阿姆斯特丹市中心建筑年代预测结果空间分布

建筑年代预测结果的空间分布
蓝色表示旧建筑被预测为新建筑，而粉色表示模型将新建筑预测为旧建筑。灰色表示预测正确。

⬇️ 图 8 ：建筑年代预测结果在 150 米网格范围的准确度

长文预警，大约需要 35 分钟。

🛩️🛩️🛩️🛩️🛩️🛩️🛩️🛩️🛩️🛩️🛩️🛩️🛩️🛩️🛩️🛩️

---

一、加载测试数据集

1.1 读取阿姆斯特丹的街景数据并选出测试集

论文中选择了 20%的图像来进行验证，也就是我们代码中的测试集。

由于我们固定了随机种子torch.manual_seed(8)，所以我们现在的测试集test_data_raw是没有被模型训练过的，也就是说，我们的模型还没有见过测试集的数据。

# 重新加载
import torch
from torch import manual_seed
from torchvision import datasets, transforms
from torch.utils.data import DataLoader, random_split, Dataset

img_root = '/root/autodl-tmp/GSV/clip' # r"../../data/GSV/clip"
all_data = datasets.ImageFolder(root=img_root)  # 不要应用tranform

# 拆分数据
train_size = int(0.8 * len(all_data))
test_size = len(all_data) - train_size

# 固定随机种子
torch.manual_seed(8)
train_data_raw, test_data_raw = random_split(all_data, [train_size, test_size])

我们看一看测试集的总数：

len(test_data_raw)
>>> 15911

15911 条测试集。

1.2 获取建筑年代类别名称和其映射关系字典

# 数据集的类别名称
class_names = all_data.classes

# 数据集的类别的字典形式
class_dict = all_data.class_to_idx
print(class_dict)

{'1653–1705': 0, '1706–1764': 1, '1765–1845': 2, '1846–1910': 3, '1911–1943': 4, '1944–1977': 5, '1978–1994': 6, '1995–2023': 7, 'pre-1652': 8}

要注意pre-1652被排到了末尾。

1.3 自定义 Dataset

为了能够进行后续的空间分析，我们需要建筑的 id 来进行定位，所以我们进一步修改CustomDataset类中的__getitem__方法，用来从Dataloader中获取数据时，不仅能返回 image_tensor 和 Label，还能返回图像文件名中的建筑 id。

图像的文件名比如“subset_1–11739–363100012571333–2023-03”使用“–“分割的字符串，建筑 id 我们只需要使用 split 从图像文件名中提取。

在预测过程中，我们会在预测中收集对应建筑 id，并在所有预测完成后将它们预测结果、真实标签一起保存到 CSV 表格文件中。我们来修改CustomDataset类：

class CustomDataset(Dataset):
    """包装PyTorch数据集以应用转换。"""
    def __init__(self, subset, transform=None):
        self.subset = subset
        self.transform = transform
        self.imgs = subset.dataset.imgs

    def __getitem__(self, index):
        img, y = self.subset[index] # 这里的y是类别的索引

        if self.transform:
            img = self.transform(img)
#####################仅修改下列代码#####################
        # 获取文件名
        file_name = self.imgs[self.subset.indices[index]][0]  # 通过传入的index来定位到图像的文件名
        # 通过split分割字符串获取文件名中的id
        id = file_name.split('--')[-2]
        return img, y, id
#####################仅修改以上代码#####################

    def __len__(self):
        return len(self.subset)

1.4 定义 transform 并加载测试集

# 只需要调整尺寸和转换为张量
test_transform = transforms.Compose([
        transforms.Resize(size=(400, 400), antialias=True),
        transforms.ToTensor()
    ])

# 加载数据集
test_data = CustomDataset(test_data_raw, transform=test_transform)

我们测试一下能不能获取到正确文件名：

# 获取数据集中的前几个项
for i in range(5):  # 检查前5项
    img, y, id = test_data[i]
    print(f"Item {i}:")
    print(f"    ID: {id}")
    print(f"    Label: {y}")
    # 如果图片是一个张量，您可以打印其形状
    print(f"    Image shape: {img.shape if hasattr(img, 'shape') else 'not a tensor'}")
    print("\n")

OUT：

Item 0:
    ID: 363100012242337
    Label: 7
    Image shape: torch.Size([3, 400, 400])
Item 1:
    ID: 363100012100709
    Label: 7
    Image shape: torch.Size([3, 400, 400])
Item 2:
    ID: 363100012069617
    Label: 5
    Image shape: torch.Size([3, 400, 400])
Item 3:
    ID: 363100012115277
    Label: 4
    Image shape: torch.Size([3, 400, 400])
Item 4:
    ID: 363100012074646
    Label: 4
    Image shape: torch.Size([3, 400, 400])

1.5 平衡数据集

虽然测试集数据也不平衡，但是测试集反映的是真实世界的情况，我认为不需要进行数据平衡，在代码中就没必要应用随机采样(WeightedRandomSampler)去平衡数据。

二、加载模型

2.1 使用 load_state_dict 加载模型

上文我们将模型结果保存为了字典：‘model.pth’。在Pytorch中，我们重新使用模型需要定义相同的模型架构，并且加载模型的字典数据。所以，我们会重新加载 densenet121 模型构架，然后将模型的最后一层分类器调整为 9 类，最后加载字典：

from torchvision.models import densenet121
from torchvision.models.densenet import DenseNet121_Weights
import torch
import torch.nn as nn

# 定义使用gpu还是cpu
device = "cuda" if torch.cuda.is_available() else "cpu"

# 加载预训练的DenseNet121模型
model = densenet121(weights=DenseNet121_Weights.DEFAULT)

# 修改最后一层的输出特征数
num_features = model.classifier.in_features
# 将原始的densenet121的1000个类别修改为9个类别，保证模型网络结构一致
model.classifier = nn.Linear(num_features, 9)

# 加载建筑年代的模型
model_path = '../models/weights_6/model_epoch_32.pth' # 修改为你保存的模型字典地址
model.load_state_dict(torch.load(model_path, map_location=device)) # map_location字段避免无GPU的电脑出错，因为此模型默认加载在cuda中。

接下来需要调整到评估模型，这样不会对模型的参数进行更新：

# 调整到eval评估模式
model.eval()
# 将模型发送到指定设备（gpu）
model.to(device)

2.2 创建 DataLoader

我在云端的 4090 显卡上运行的，你可以根据情况调整BATCH_SIZE和num_workers。

BATCH_SIZE = 128

test_loader = DataLoader(test_data, batch_size=BATCH_SIZE, shuffle=False, num_workers=10)
len(test_loader) # 等于数据集总长度除以每批次的大小（BATCH_SIZE）

OUT:

166

BATCH_SIZE大写代表常量，即不需要修改的变量。在 Python 的官方风格指南 PEP8 建议使用全部大写的方式来命名常量。如果你对 PEP8 感兴趣可以阅读PEP8 官方文档。

三、开始预测

3.1 对整个测试集进行预测

我们预测图像的最终目标是获取每个图像的预测标签，用来对比是否和真实标签相等，从而进行接下来的分析。

我们先定义三个空列表，储存真实标签、预测标签和建筑 id：

true_labels = []
pred_labels = []
ids_list = []

然后我们会用 for 循环遍历 166 个DataLoder，每个DataLoder中有 128 个图像（BATCH_SIZE 的大小）：

from tqdm import tqdm # 用来显示进度条

with torch.inference_mode():
    for images, labels, id in tqdm(test_loader, desc="Predicting"): # 遍历test_loader会使用__getitem__方法,返回img, y, id

        # 将数据移动到GPU上（如果可用）
        images, labels = images.to(device), labels.to(device)

        # 运行模型以获取预测（向前传递）
        outputs = model(images)

        # 使用argmax获取最大值的索引
        test_pred_labels = outputs.argmax(dim=1)

        # 将预测标签、真实标签和建筑id添加到列表中
        true_labels.extend(labels.cpu().numpy())
        pred_labels.extend(test_pred_labels.cpu().numpy())
        ids_list.extend(id)

    # 如果您想查看这一批的结果，可以打印或处理这些列表
    # print("真实标签", true_labels)
    # print("预测标签:", pred_labels)
    # print("建筑id", ids_list)

在我们的代码中，true_labels 和 pred_labels 是一维数组，ids_list 是列表。

1）logits > pred_labels

重点说一下如何通过模型的预测结果(output，称为logits，原始输出)得到它的预测标签(test_pred_labels)：

在我们的多类分类问题中，模型的输出是一个的概率分布，表示 9 个类别的预测概率。变量 outputs 是一个二维张量，其中包含了批次中每个样本对应每个类别的预测分数或概率。第一维（dim=0）表示批次中的样本索引。第二维（dim=1）表示每个类别的预测分数。

例如，我们有一个批次大小为 32 的数据，且分类问题有 9 个类别，那么 outputs 的大小是 [32, 9]使用 argmax 函数: argmax(dim=1) 在类别的维度上找到最大值索引。然后，在这种情况下，在它每一行（对应一个样本的所有类别预测）上找到最大值的索引。这个索引实际上是模型预测的类别标签（0-8）。

所以最终， test_pred_labels 包含每个输入样本的预测标签。这些标签是根据模型给出的最高分数（概率）选择的类别。如果你有一个 [32, 9] 的输出，那么 test_pred_labels 将是一个长度为 32 的一维数组，每个元素都是 0 到 8 之间的一个整数（对应 9 个类别）。

2）将测试集的数据保存为表格

可以将预测结果保存为表格，方便后续加载。

import pandas as pd

# 创建一个数据框来保存文件名和预测
df_predictions = pd.DataFrame({
    'id': ids_list,
    'prediction': pred_labels,  # 这是之前收集的预测列表
    'true_label': true_labels  # 这是之前收集的真实标签列表
})

# 将数据框写入CSV文件
df_predictions.to_csv('predictions_with_building_age_model_6_on_test_data.csv', index=False)

3.2 可视化某一批次图像的预测结果

我们直接 matplotlib 用绘制结果，但是，数据集太大了，我们只想绘制某一批次的数据。所有我们先从DataLoader取出一些数据：

1） 使用迭代器

我们使用从DataLoader中抽取第一批数据来进行绘制。但是DataLoader并不是列表，也不是迭代器，是一个 Pytorch 的DataLoader对象，为了能够从中取出数据，需要先使用iter()将DataLoader转换为迭代器（也称为生成器，它的特性是不会将数据全部加载到内存，调用它的时候才会进入内存），然后进行for循环遍历，或者直接使用next()获取迭代器的下一个批次的数据，第一次调用next()则获取第一批数据。

继续使用next()会获取第二批数据，以此类推。

看看我们的代码实现，在下列代码中，如果你是在 jupyter notebook 中运行，我们先将 num_workers 设为 0 以避免多线程 bug：

自定义数据集时并且自定义数据集的函数不在当前单元格、同时 num_workers 大于 0 就会出现此 bug：
当您在使用 DataLoader 时设置 num_workers 大于 0 以使用多个子进程加载数据时，PyTorch 使用 multiprocessing 来创建这些子进程。但是，multiprocessing 需要能够从主进程中找到并加载任何自定义函数或类，这在 Jupyter Notebook 或其他交互式环境中可能会出问题。
两种解决方法： 1. 设置num_workers=0 2. 将自定义数据集CustomDatasets()放入__name__ == '__main__'中。

BATCH_SIZE = 8 # 此处代表你要绘制的多少图像的预测结果

test_loader = DataLoader(test_data, batch_size=BATCH_SIZE, shuffle=False, num_workers=0)

if __name__ == '__main__': # 以尝试将启动训练过程的代码放入此保护块中。这有助于防止 multiprocessing 在它不应该这样做的时候启动新进程。
    test_data_iter = iter(test_loader)
    test_samples, test_labels, ids = next(test_data_iter) # next() 函数是用来获取迭代器的下一个批次的数据
    print(test_samples.shape, test_labels.shape, len(ids) )
    # print(test_samples, test_labels, ids) # 可以打印具体结果

OUT:

torch.Size([8, 3, 512, 512]) torch.Size([8]) 8

test_labels是 8 个一维数组代表真实标签，ids_list 是自定义 dataset 返回的列表，此时返回包含 8 个建筑 id 的列表，如果想保持他们的一致性，我们也可在自定义数据集中将 ids_list 定义为一维数组。

其中，test_samples是一个有四个维度的张量，每个维度的大小分别为 8、3、512 和 512，[BATCH_SIZE, C, Height, Width]，这种维度设置通常是在深度学习框架中使用的“NCHW”（或“BHWC”）数据格式，其中 N/B 表示批量大小、C 表示通道数、H 表示高度、W 表示宽度。

2） 将预测的标签从索引转到其真实名称

我们要在图片上显示出建筑 id、预测和真实类别，但是现在的 test_labels 还是索引值，我们要从 class_dict 获取真实年代标签进行替换，方便阅读：

class_dict ={'1653–1705': 0,
 '1706–1764': 1,
 '1765–1845': 2,
 '1846–1910': 3,
 '1911–1943': 4,
 '1944–1977': 5,
 '1978–1994': 6,
 '1995–2023': 7,
 'pre-1652': 8}
# 创建一个值到键的反向映射 方便取值
reverse_dict = {value: key for key, value in class_dict.items()}
reverse_dict

out：

{0: '1653–1705',
 1: '1706–1764',
 2: '1765–1845',
 3: '1846–1910',
 4: '1911–1943',
 5: '1944–1977',
 6: '1978–1994',
 7: '1995–2023',
 8: 'pre-1652'}

# 通过反向映射，我们可以直接用值获取键，现在获取字典中键为整数0的值
key_with_value = reverse_dict[0]
key_with_value

OUT:

'1653–1705'

3） 加载中文字体

matplotlib 默认不支持中文字体，需要在代码中定义：

from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"C:\Windows\Fonts\simhei.ttf", size=10) # 选择你系统中的中文字体的路径，并且通过size定义大小

4） 预测并绘制

from tqdm import tqdm
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
import numpy as np
# 存储真实标签 预测标签 文件名
true_labels = []
pred_labels = []
ids_list = []

images_so_far = 0
fig = plt.figure(figsize=(10, 20))
font = FontProperties(fname=r"C:\Windows\Fonts\simhei.ttf", size=10)
num_images = 8

# # 可以继续调用获取第二批数据
next_samples, next_labels, ids = next(test_data_iter)

with torch.inference_mode():
        # 将数据移动到GPU上（如果可用）
        images, labels = next_samples.to(device), next_labels.to(device)
        # 运行模型以获取预测（向前传递）
        outputs = model(images)
        # 使用argmax获取最大值的索引
        test_pred_labels = outputs.argmax(dim=1)

        # 选择要显示的图片
        for j in range(images.size()[0]):
            images_so_far += 1
            ax = plt.subplot(num_images//2, 2, images_so_far)
            ax.axis('off')
            pred_label = reverse_dict.get(int(test_pred_labels[j]))
            true_label = reverse_dict.get(int(labels[j]))
            ax.set_title(f'建筑id:{ids[j]}\n预测类别: {pred_label}，    真实类别: {true_label}', fontproperties=font)

            # 将图形转移到cpu 并且更改通道顺序 从[C, Height, Width]更改为[Height, Width, C]
            image = images.cpu().data[j].numpy().transpose((1, 2, 0))

            plt.imshow(image)

四、混淆矩阵、召回率、精确度、F1 分数

4.1 概念解释

1）混淆矩阵

混淆矩阵（Confusion Matrix）是在分类问题中用于评估模型性能的一种表格形式。它以实际类别（真实标签）和预测类别为基础，将样本的分类结果进行统计和总结。混淆矩阵的每一行代表了真实类别，每一列代表了预测类别。

混淆矩阵的常见形式如下，我写成英文更容易理解：

用一个例子理解：

混淆矩阵中的四个关键术语是：

• True Positive (TP)： 即实际为正且被预测也为正的样本数。图中 True Positives (TP) = 86。
• False Positive (FP)： 即实际为负但被错误地预测为正的样本数。图中 True Negatives (TN) = 79。
• False Negative (FN)： 实际为正但被错误地预测为负的样本数。图中 False Positives (FP) = 12。
• True Negative (TN)： 即实际为负且被预测为负的样本数。图中 False Negatives (FN) = 10。

基于上述情况，我们可以定义（召回率、精确度和 F1 分数）：

2）召回率 (Recall):

• **概念：**召回率衡量了所有真实为正的样本中，被模型正确预测为正的比例。

• 公式：

$$
Recall = True Positives (TP) / (True Positives (TP) + False Negatives (FN))
$$

• **作用：**召回率特别适用于那些错过真实正样本的代价很高的情境，例如疾病诊断。在这种情况下，我们更希望模型能够捕获所有可能的正样本，即使这意味着会有一些误报。
• **在我们的例子中：**Recall = 86 / (86 + 10) = 0.8983 = 89.83%

3）精确度 (Precision):

• **概念：**精确度衡量了被模型预测为正的样本中，真实为正的比例。

• 公式：

$$
Precision = True Positives (TP) / (True Positives (TP) + False Positives (FP))
$$

• **作用：**精确度特别适用于那些误报代价很高的情境，例如垃圾邮件检测。在这种情况下，我们不希望误将正常的邮件标记为垃圾邮件。
• **在我们的例子中：**Precision = 86 / (86 + 12) = 0.8775 = 87.75%

4）F1 分数 (F1 Score)

• **概念：**F1 分数是召回率和精确度的调和平均值，它试图在召回率和精确度之间找到一个平衡。

• 公式：

$$
F1 Score = 2 * Precision * Recall / (Precision + Recall)
$$

• 作用：当我们需要同时考虑召回率和精确度时，F1 分数是一个很好的指标。它特别适用于那些正负样本不均衡的数据集。
• **在我们的例子中：*F1-Score = (2 0.8775 * 0.8983) / (0.8775 + 0.8983) = 0.8877 = 88.77%

总之，选择哪个指标取决于具体的应用场景和业务需求。

在某些情况下，我们可能更关心召回率（例如，医院确保所有患者都被正确诊断），而在其他情况下，我们可能更关心精确度（例如，确保只有真正的垃圾邮件被标记）。

当我们需要同时考虑召回率和精确度时，F1 分数提供了一个综合的评估指标。

4.2 读取预测结果

我们将使用sklearn提供的工具来计算混淆矩阵、召回率、精确度和 F1 分数。

import pandas as pd
df = pd.read_csv('predictions_with_building_age_model_6_on_test_data.csv')
df.head()

true_labels = df['true_label'].tolist()
pred_labels = df['prediction'].tolist()

因为我们的原始 class_dict 的pre-1652没有按照时间顺序排在最前，所以需要更改一下顺序，我们将“pre-1652”的索引从 8 更改为 0，其他类别相应地向后移动：

updated_class_dict = {
    0: '–1652',  # 这个现在是第一个
    1: '-1706',
    2: '-1765',
    3: '-1846',
    4: '-1911',
    5: '-1944',
    6: '-1978',
    7: '–1995',
    8: '–2023'  # 这个现在是最后一个
}

updated_true_labels = [(label + 1 if label < 8 else 0) for label in true_labels]
updated_pred_labels = [(label + 1 if label < 8 else 0) for label in pred_labels]

4.3 使用 sklearn 创建混淆矩阵

from sklearn.metrics import confusion_matrix

conf_matrix = confusion_matrix(updated_true_labels, updated_pred_labels)
print("Confusion Matrix:")
print(conf_matrix)

4.4 使用 seaborn 进行可视化

import seaborn as sns
import matplotlib.pyplot as plt

# 设置可视化效果的样式
sns.set(style='whitegrid', palette='muted')

# 将混淆矩阵转换为DataFrame
conf_matrix_df = pd.DataFrame(conf_matrix, index=class_names, columns=class_names)

# 创建热图
plt.figure(figsize=(8, 6))
sns.heatmap(conf_matrix_df, annot=True, fmt='d', cmap='Blues')
plt.ylabel('Actual')
plt.xlabel('Predicted')
plt.title('Confusion Matrix')
plt.show()

4.5 通过混淆矩阵分析模型预测结果

我们来分析一下我们的混淆矩阵：

1. 主对角线：从左上角到右下角的数字表示模型正确预测的数量。
2. 横轴和纵轴：横轴（Predicted）代表模型的预测类别，纵轴（Actual）代表实际的类别。
3. 矩阵中的其他元素：表示模型的误判数量。
4. 颜色：颜色深浅代表了数值的大小。深色代表数字大，浅色代表数字小。从图中可以看出，对角线上的颜色比较深，说明模型在这些类别上的预测较为准确。而其他位置的颜色较浅，表示误判的数量相对较少。

基于这个混淆矩阵，我们可以得出一些结论：

• 主对角线表现：大部分的样本被正确地分类，这可以从对角线上的深蓝色区域看出。这说明模型在许多类别上的预测都是准确的。模型在“-1911”、“-1944”、“-1978”和”-2023“类别模型都有较高的准确率。
• 误分类情况：尽管某些类别的预测准确性较高，但仍然存在一些误分类的情况。“-1911”、“-1944”类别，除了主对角线上的 5781 外，还有其他非对角线上的值如 373、180、49 和 47，这表明这些实例被误分类到其他类别。
• 难以分类的类别：某些类别如"-1652"、“-1706”、“-1765”、"-1846"在主对角线上的值相对较小，并且有多个非对角线的值也相对较大，这意味着模型在这些类别上的预测存在较大的困难。

4.6 使用 sklearn 生成各种分类指标

分类报告（classification report）为我们提供了每个类别的主要分类指标的细分，这有助于我们理解模型在预测每个特定类别时的性能：

# 借助混淆矩阵计算各种分类指标(召回率、精确度和F1分数)
class_names = list(updated_class_dict.values())
report = classification_report(updated_true_labels, updated_pred_labels, target_names=class_names)

print("\nClassification Report:")
print(report)

4.7 使用分类报告分析模型预测结果

从这个分类报告中，我们可以看出：

1. 准确率 (Precision): 是模型预测为正例中实际也为正例的比例。该模型在不同的类别上的准确率有很大的差异。例如，“-1944” 类别的准确率高达 0.93，表示对于这个类别的预测很准确。然而，“-1652” 类别的准确率仅为 0.24，表示在预测为这个类别的结果中，有相当一部分是误判的。
2. 召回率 (Recall): 是模型正确预测的正例占所有实际正例的比例。对于 “-1911” 类别，召回率达到了 0.83，这意味着模型能够捕获大部分真实的 “-1911” 样本。但对于 “-1706” 类别，召回率仅为 0.25，表示很多真实的 “-1706” 样本都没有被正确预测。
3. F1 得分 (F1-Score): 是准确率和召回率的调和平均值，用于考虑准确率和召回率之间的平衡。例如，“-1944” 类别的 F1 得分为 0.91，表现很好。然而，“-1652” 类别的 F1 得分为 0.22，这意味着这个类别的预测性能较差。
4. 支持 (Support): 是每个类别的样本数量。从支持（support）列可以看出，类别 “-1944” 的样本数量最多，达到 6488，而类别 “-1706” 的样本数量最少，仅为 169。这种数据不平衡可能影响到模型的预测性能，特别是对于样本数量较少的类别。
5. 总体准确率 (Accuracy): 是模型正确预测的样本占总样本的比例。模型的整体准确率为 0.82，表示模型在所有的预测中有 82% 是正确的。
6. 宏平均 (Macro Avg): 是所有类别的平均准确率、召回率和 F1 得分。此模型的宏平均精确度、召回率和 F1 得分都为 0.59。这意味着在所有类别上，模型的平均性能是相对一致的。
7. 加权平均 (Weighted Avg): 考虑到每个类别的样本数量，模型的加权平均精确度、召回率和 F1 得分都为 0.82。这与总体准确率相符，说明模型在样本数量较多的类别上的性能较好。

综上分类报告，我们可以进行总结模型的预测结果：

1. 整体性能：模型的总体精确度为 82%，这意味着约 82% 的预测是正确的，这是一个相对较高的准确度。
2. 最佳性能类别：模型在 “-1944”, “-1978”, 和 “-2023” 这三个类别上的预测性能较好，其中 “-1944” 的 F1-score 为 0.91，是所有类别中最高的。
3. 需要改进的类别：模型在 “-1652”, “-1706”, 和 “-1846” 这三个类别上的预测性能较差，特别是在 “-1652” 上，其 F1-score 为 0.22，这意味着该类别的预测性能有待提高。
4. 不平衡的性能：虽然 “-1911” 类别的 recall 达到了 0.83，但其 precision 为 0.74，这意味着模型更倾向于将样本分类为这个类别，但同时这也可能导致其他类别的误分类增加。

总的来说，虽然模型在某些类别上的预测性能较好，但在其他类别上仍存在改进的空间。但在某些类别上可能需要进一步优化。

对比论文中的模型评估结果(下图)，我们的模型不够完美，差距还比较大：

将我们的混淆矩阵转化为百分数：

虽然我们和作者的数据集不一样，但是我的研究方法是没错的，如果后期学到更多处理技巧，或者发现错误，我会在我的博客更新，也欢迎看到这里的大佬能给出一些优化意见！！！

4.8 可以进一步优化的地方

1. 街景数据集：获取更多的-1652, -1706, -1765 和-1846 四个类别街景图。
2. 进一步达到数据平衡: 在不平衡的分类问题中，可以使用过采样、欠采样或合成数据技术，如 SMOTE，来平衡数据。
3. 建筑足迹数据有待优化：正如论文中所述：“BAG 数据集中建筑物的施工年份被定义为“建筑物最初或将在施工准备就绪时交付的年份”。建筑物的改建、扩建和增建不会改变原来的建造年份。这种限制反映在图 8 中，我们的年龄纪元估计模型对阿姆斯特丹市中心的预测不太准确，因为市中心的建筑物有更高的可能性进行翻新。”——更老的建筑也很有可能被翻新。

五、类激活映射

观察模型观察的是建筑的哪一个部分，有助于了解不同建筑年代的建筑局部的差异点在哪？是建筑外墙的材料？门或窗的形式？我们可以使用一种叫计算机视觉中的类激活映射（Class Activation Mapping，简称 CAM）技术的方法查看模型关注以上哪些方面。

类激活映射（Class Activation Mapping，简称 CAM）是一种在计算机视觉中广泛使用的技术，特别是在深度学习和卷积神经网络（CNN）的上下文中。它用于可视化输入图像的哪些部分被模型用来识别特定的类别。换句话说，CAM 帮助我们理解模型的决策过程，特别是模型是如何从视觉信息中“学习”并做出分类决策的。

我们可以自己“手搓”一个 CAM，也可以直接用别人的开源项目，我在 Github 发现了两个高星 CAM 开源项目：frgfm/orch-cam和zhoubolei/CAM。

这两个库该选择哪一个？

“frgfm/torch-cam” 主要是一个为PyTorch用户设计的库，提供了一个模块化和易于集成的 CAM 解决方案。这个库是为PyTorch框架设计的，它提供了一个高度模块化和可定制的系统，可以轻松地与现有的PyTorch模型集成。如果你已经在使用PyTorch并且需要一个能够快速集成并且具有良好文档支持的解决方案，那么这个库可能是最佳选择。

而 “zhoubolei/CAM” 则更像是一个完整的研究项目，包含了从理论解释到实际应用的所有内容，使用的是 MATLAB 和 Caffe 框架，并提供了多种预训练模型。这个库提供了一个更全面的解决方案，包括预训练模型、详细的使用说明和评估脚本。它更多地侧重于教育和研究，提供了对 CAM 理论的深入理解。但是，它主要基于MATLAB和Caffe，我不熟悉这些工具。

最终我选择了"frgfm/torch-cam" 库，它使用了面向对象的方法，定义了一个基础类 _CAM，用于实现类激活映射（CAM）的核心功能。这种设计允许扩展不同类型的 CAM 方法。而且用户可以指定观察的目标层，或者让系统自动选择。

5.1 使用"frgfm/torch-cam" 库对单个图像进行测试

1）安装

可以使用 pypi 安装软件包的最新稳定版本：

pip install torchcam
# 或使用conda：
conda install -c frgfm torchcam

2）定义 CAM 的目标层

# 设置CAM
# step:1使用上文定义的模型并且设置为评估模式
from torchinfo import summary
summary(model=model,
        input_size=(32, 3, 300, 300), # make sure this is "input_size", not "input_shape"
        # col_names=["input_size"], # uncomment for smaller output
        col_names=["input_size", "output_size", "num_params", "trainable"],
        col_width=20,
        row_settings=["var_names"]
)

在使用类激活映射（CAM）的情况下，通常会选择网络中的最后一个卷积层或与最后一个卷积层紧密相关的层作为目标层。这是因为这些层通常包含关于目标类的空间信息，这对于理解网络如何“看到”和识别特定特征是非常有用的。

在我们提供的 DenseNet 模型中，应该将目标层指定为最后一个 DenseBlock 或其内部的最后一个 DenseLayer。具体来说，这可能是denseblock4 或 denseblock4 中的最后一个denselayer（例如 denselayer16）。这些层在空间分辨率上保留了足够的信息，同时包含了对模型决策至关重要的特征表示。

3）设置 CAM extractor

from torchcam.methods import SmoothGradCAMpp
cam_extractor = SmoothGradCAMpp(model, target_layer='features.denseblock4') # 默认情况下，检索CAM的层被设置为最后一个非简化卷积层。如果希望研究特定的层，请在构造函数中使用 target_layer 参数。

4）读取图像并可视化

from torchvision.io.image import read_image
from torchvision.transforms.functional import normalize, resize, to_pil_image

img = read_image("../../data/GSV/default_image.png")
img.shape

OUT：

torch.Size([4, 512, 512])

我们读取的 png 也有 alpha 通道，需要删除：

# 移除png的alpha通道
img = img[:3, :, :]
img.shape

OUT：

torch.Size([3, 512, 512])

转换图像为 tensor

# Preprocess it for your chosen model
input_tensor = normalize(resize(img, (300, 300)) / 255., [0.485, 0.456, 0.406], [0.229, 0.224, 0.225])
input_tensor.shape

OUT：

torch.Size([3, 300, 300])

在 Pytorch 中最容易出错的地方之一就是 tenor 的形状。我们需要时刻注意。

5）执行 CAM

with SmoothGradCAMpp(model, target_layer='features.denseblock4') as cam_extractor:

    input_tensor = input_tensor.to(device)

    # 把图像的tenso喂给模型
    out = model(input_tensor.unsqueeze(0))
    # 通过传递类索引和模型输出来检索CAM（类激活映射）。
    activation_map = cam_extractor(out.squeeze(0).argmax().item(), out)
activation_map

activation_map[0].shape为torch.Size([1, 9, 9])，第一个维度代表 cam 的值，在 0-1 之间，越高的代表模型的“关注度”越高，后两个通道是图片的分辨率：9x9 的图。

# 可视化你的热图 将其覆盖在输入图像上：
import matplotlib.pyplot as plt
from torchcam.utils import overlay_mask

# 调整CAM的大小使其能覆盖图像
result = overlay_mask(to_pil_image(img), to_pil_image(activation_map[0].cpu().squeeze(0), mode='F'), alpha=0.5)
# 绘制
plt.imshow(result); plt.axis('off'); plt.tight_layout(); plt.show()

3）我们将 CAM 嵌入评估流程中

参考以下代码：

详见4.1.3-建筑年代模型评价，可以当做练习。

with SmoothGradCAMpp(model, target_layer='features.denseblock4') as cam_extractor:
        # 将数据移动到GPU上（如果可用）
        images, labels = next_samples.to(device), next_labels.to(device)

        # 运行模型以获取预测（向前传递）
        outputs = model(images)

        # 使用argmax获取最大值的索引
        test_pred_labels = outputs.argmax(dim=1)

        # 选择分析CAM
        for j in range(images.size()[0]):
            if images_so_far == num_images:
                break  # 如果我们达到了所需的图像数量，就停止

            images_so_far += 1
            ax = plt.subplot(num_images//2, 2, images_so_far)
            ax.axis('off')

            # 获取原始图像和对应的输出
            img = images[j] # torch.Size([3, 300, 300])
            output = outputs[j] # torch.Size([9])
            # print(img.shape, output.shape)

            # 获取CAM并将其应用到原始图像上
            activation_map = cam_extractor(output.squeeze(0).argmax().item(), output)
            # print(" to_pil_image(img):", to_pil_image(img)) # <PIL.Image.Image image mode=RGB size=300x300 at 0x1F3C85FFFD0>
            # print("activation_map:", activation_map) # tensor列表
            # print("activation_map[0]:", activation_map[0].shape) # torch.Size([8, 9, 9])
            # [0, :, :] 选择第一个元素，并丢弃第一个维度
            #print("activation_map[0][0, :, :]：", activation_map[0][0, :, :].shape) # torch.Size([9, 9])
            result = overlay_mask(to_pil_image(img), to_pil_image(activation_map[0][0, :, :].cpu(), mode='F'), alpha=0.5)

            plt.imshow(result)

4）分别绘制 9 个年代的 CAM 图

参考论文的图 10（下图）：

1. 左侧四张小图是将 CAM 叠加在原始图像上。图像的红色区域主要覆盖一楼和二楼之间的窗户或门。
2. 右侧四张小图：根据 CAM 裁剪的图像显示了窗户的演变。早期的窗户通常框架较宽，装饰较多，而且较窄。最近的窗户样式以方形和水平形状为特点，框架更薄，装饰更少，深度更小。

总的来说，图 10 显示了 的 CAM 结果突出显示了模型关注的对象集中在一楼和二楼之间的窗户或门，而不是建筑物的随机部分。街景图像前面的汽车、自行车和行人通常会被忽略，因为它们与建筑年龄无关。本质上，模型从图像中学习有效的特征，并自动忽略图像中的不相关信息。

制作这张图的方法很简单，我们挑选一些照片之后，通过 PS 绘制出下图(图 10 利用类激活（CAM）图像观察不同类别模型的关注点）：

六、空间分布

6.1 建筑年代预测结果的空间分布

1）绘制思路

参考文中 图 7 绘制市中心的建筑年代预测结果图。

蓝色表示旧建筑被预测为新建筑，而粉色表示模型将新建筑预测为旧建筑。灰色表示预测正确。

我们可以参考上图进行制作，流程大概是：对所有的建筑进行预测——对真实年代和预测的年代的类别进行差值计算——将上一步的结果和建筑足迹的空间数据进行连接——提取出市中心的范围，设置符号系统然后出图。

2）处理预测结果

a.对训练集进行预测

我们利用“三、进行预测”的方法对训练集进行预测

b.合并预测结果

## 读取数据
import pandas as pd

df1 = pd.read_csv('predictions_with_building_age_model_6_on_test_data.csv')
df2 = pd.read_csv('predictions_with_building_age_model_6_on_train_data.csv')
df = pd.concat([df1, df2])
df.head()

c.将合并结果保存

# 保存合并后的原始结果
df.to_csv('predictions_with_building_age_model_6_on_all_data.csv', index=False)

predictions_with_building_age_model_6_on_all_data.csv 可以点击👉 此处下载

d.重新排列类别顺序

按照 4.1.3 的方法重新排列类别（我们将“pre-1652”的索引从 8 更改为 0，其他类别相应地向后移动），方便计算:

#我们需要更新标签以反映这个新顺序
# 这意味着我们需要将所有的8替换为0，然后将其他所有数字加1（因为我们把'–1652'放在了最前面）
df['true_label'] = df['true_label'].apply(lambda x: x + 1 if x < 8 else 0)
df['prediction'] = df['prediction'].apply(lambda x: x + 1 if x < 8 else 0)
df.head()

e.计算差值

计算prediction和true_label的差值：

## 新增一列差值
df['diff'] = df['prediction'] - df['true_label']
df.head()

注意：❗ 负数为预测为旧建筑，正数为预测为新建筑，0 为预测正确。差值为偏离真实值的程度。

d.将类别索引更换为类别名称

updated_class_dict = {
    0: '–1652',  # 这个现在是第一个
    1: '-1706',
    2: '-1765',
    3: '-1846',
    4: '-1911',
    5: '-1944',
    6: '-1978',
    7: '–1995',
    8: '–2023'  # 这个现在是最后一个
}
df['true_label'] = df['true_label'].apply(lambda x: updated_class_dict[x])
df['pre_label'] = df['prediction'].apply(lambda x: updated_class_dict[x])

df.drop(columns=['prediction'], inplace=True)
df.head()

e.在合并前处理id列

通过后续检查发现，预测结果 df 中的id与建筑足迹中identificatie的特征有所不同:预测结果 df 中的 id 列是整数，而建筑足迹Amsterdam_buildings_Project中的 id 列数据类型是 16 个字符，并在不足 16 位时用前导零填充：

df.id

# 将 id 转换为字符串，确保其长度为 16 个字符，必要时用前导零填充。
df['id'] = df['id'].apply(lambda x: f"{int(x):016}")
df['id']

保存：

df.to_csv('predictions_with_building_age_diff_citywide.csv', index=False)

3）读取空间数据

import geopandas as gpd
gdb = "../../5-ArcgisPro工程/建筑风格和年代深度学习.gdb"
lr_name = 'Amsterdam_buildings_Project'
gdf = gpd.read_file(filename=gdb, layer=lr_name)
# 选取需要的列
gdf = gdf[['identificatie', 'geometry']]

4）合并预测结果和建筑足迹数据

df_merge = pd.merge(left=gdf, right=df, left_on='identificatie', right_on='id', how='outer')
df_merge.head()

说明一下，GeoDataFrame和一个DataFrame合并需要用pandas中的merge方法，输出是否为GeoDataFrame取决于left字段是否为GeoDataFrame对象。

# 看一下diff非空值所占的比例，即有效数据的比例
df_merge['diff'].notnull().sum() / df_merge.shape[0]

OUT：

0.4870035261276262

真正被预测的建筑占比 48%，不是很高，大多数都被街景筛选掉了，进一步说明了我们的项目需要重新获取缺失的街景图片。

# 看一下diff的分布
df_merge['diff'].hist()

绝大多数建筑都被正确预测了。负数为预测为旧建筑，正数为预测为新建筑

5）ArcGIS Pro 绘图空间分布图

我喜欢在 ArcGIS Pro 中绘图，能实时可视化图面效果，并且最终的渲染效果也比在 python 中号，所以我们把gdf_merge导出到 GIS 的文件地理数据库：

df_merge.to_file(filename=gdb, layer='predictions_with_building_age_diff_city_wide')

调整符号系统，将地图放入布局中，加上图例：

粉色代表新建筑被预测为旧建筑，蓝色代表旧建筑被预测为新建筑。

---

6.2 绘制建筑年代预测结果在 150 米网格范围的准确度

我们要复现论文中的图 8：

图片上表现的是预测的精准程度在 150m 的网格上的空间分布，图中可以看出：市中心的错误率高于郊区。这可能是由于建筑年龄的高度多样性、市中心旧建筑的频繁改造以及上述对改造建筑的严格规定所致。阿姆斯特丹郊区没有明显的空间格局，这表明分类结果的空间相关性很小。为了证明空间相关性小，作者还计算了莫兰指数，城市郊区结果的 Moran’s I 为 0.27。 Moran’s I 的范围为 -1 到 1，其中 -1 表示完美分散dispersion ，0 表示完美随机性 randomness，1 表示完美聚类clustering。我们预测的阿姆斯特丹郊区的莫兰指数为 0.27，建筑年代表现出较弱的空间自相关性。

1）获取含建筑位置信息和预测准确度的建筑数据

import geopandas as gpd
import pandas as pd

gdb = r"../../5-ArcgisPro工程/建筑风格和年代深度学习.gdb" #
lr_name = 'predictions_with_building_age_diff_city_wide' # 在6.1中保存的全市范围建筑准确度的结果
gdf = gpd.read_file(filename=gdb, layer=lr_name)
gdf.head()
# 如果你在同一个notebook中不用重新读取

2）利用 ArcPy 的创建渔网工具创建 150 米的空间网格

a.创建渔网函数解析

官方帮助文档在：创建渔网 (数据管理)，同时在《利用 ArcGIS_Python 制作考虑路况的交通等时圈》文也讲过，可以点击文章查看，本文中的代码是在此基础上修改的。

我们回忆一下关键信息：

1️⃣ 创建渔网函数的参数：

• out_feature_class：包含由矩形像元组成的渔网的输出要素类。
• origin_coord：矩形框的左下端点为原点
• y_axis_coord：此点与原点的连线用于判断旋转的角度 我们不用旋转所以定义为原点正上方的点
• cell_width：网格的宽度
• cell_height：网格的高度
• number_rows：填’0’,代表留空，由 cell_width 和 cell_height 决定
• number_columns：填’0’，由 cell_width 和 cell_height 决定
• corner_coord：填’0’
• labels：‘LABELS’
• template：以空格分隔的坐标字符串 - 将使用指定渔网的范围。坐标以 x-min，y-min，x-max，y-max 的顺序表示。
• geometry_type：生成面

2️⃣ 创建渔网返回的结果：

• out_feature_class：包含由矩形像元组成的渔网的输出要素类。
• out_label：创建一个新的点要素类，其中每个渔网像元中心都具有标注点。如果选中了创建标注点参数（Python 中的 labels = ‘LABELS’），则会创建一个新的点要素类，其中每个渔网像元中心都具有标注点。此要素类的名称以 _label 为后缀并与输出要素类相同，且创建于同一位置。

b.定义 origin_coord 和 y_axis_coord

从上文 gdf 中获取边界:

from shapely.geometry import box

# total_bounds返回一个元组，包含最小的x，最小的y，最大的x，最大的y
x_min, y_min, x_max, y_max = gdf.total_bounds
# 创建一个边界框
bbox_geometry = box(x_min, y_min, x_max, y_max)
# 创建一个GeoDataFrame，将边界框作为几何图形
gdf_bbox = gpd.GeoDataFrame({'geometry': bbox_geometry}, index=[0]).set_crs(32631)
gdf_bbox

得到了一个覆盖整个建筑足迹的多边形。

定义 origin_coord 和 y_axis_coord 空格分割：

origin_coord = f"{x_min}  {y_min}"
y_axis_coord = f"{x_min}  {y_max}"
origin_coord, y_axis_coord

('617784.8387999535  5798425.739600182',
 '617784.8387999535  5810551.856800079')

c.导入 arcpy 创建渔网

我们可以用 python 创建渔网，这样就不需要借助 ArcGIS Pro 了，但是相比于 ArcGIS Pro 提供的渔网工具，后者能同时生成多边形和多边形的中心点，所以我选择了最简单的工具。

在此之前：

• 你需要安装好 ArcGIS Pro 并打开软件中的Jupyter Notebook窗口。

  "1、ArcGIS Pro 的安装

  对于新手，可以选择方式一试用。

  • 方式一：官方试用 21 天

  • 方式二：参考麻辣 GIS分享的ArcGIS Pro 3.0 完整安装教程和安装包。

• 需要克隆一个新的 arcgis pro 的环境然后安装geopandas和shapely。

以上两点在文章：一、Arcpy 介绍和安装【ArcGIS Python 系列】有具体说明。

# 导入arcpy
import arcpy

# 设置工作空间
arcpy.env.workspace = gdb
arcpy.env.overwriteOutput = True # 开启工作空间输出可覆盖选项

scales = 150 # 网格的单元尺度 单位为米
extent = f"{x_min} {y_min} {x_max} {y_max}" # '617784.8387999535 5798425.739600182 640252.4324002266 5810551.856800079'
# 设置空间参考对象
arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(32631)
out_fcs = 'fishnet' # 有bug 不能数字开头

# 设置工作空间
arcpy.management.CreateFishnet(out_feature_class = out_fcs, # 包含由矩形像元组成的渔网的输出要素类。
                               origin_coord = origin_coord, # 矩形框的左下端点为原点
                               y_axis_coord = y_axis_coord, # 此点与原点的连线用于判断旋转的角度 我们不用旋转所以定义为原点正上方的点
                               cell_width = scales,
                               cell_height = scales,
                               number_rows = "0", # 留空，由cell_width和cell_height决定
                               number_columns = "0", # 留空，由cell_width和cell_height决定
                               corner_coord = None, # 对角坐标不填写
                               labels = "LABELS",
                               template = extent, # 以空格分隔的坐标字符串 - 将使用指定渔网的范围。坐标以 x-min，y-min，x-max，y-max 的顺序表示。
                               geometry_type = "POLYGON" # 生成面
                               )

不报错并且软件左侧有生产了两个新要素就成功了。

# 定义新产生的点要素的名称
out_label = out_fcs + "_label"

你可以尝试一下用 python 的 geopandas 和 shapely 如何绘制渔网。

3）空间链接

geopandas和ArcPy都有空间连接的功能，但是geopandas的空间连接功能更强大，而且方便进行数据统计，所以我们使用geopandas的空间连接功能。

在选择处理工具的时候，我也不喜欢用ArcPy的游标进行数据处理，因为ArcPy的游标处理数据的速度太慢了，而且代码也不够简洁。

所以我更喜欢用geopandas基于dataframe的处理方式，索引、切片、查询等操作都很方便。但是ArcPy中提供一些实用的工具，比如同时创建渔网和渔网标注点，这个功能在 geopandas 中没有，所以我们用ArcPy创建渔网，然后用geopandas进行空间连接。另外，制图的时候，我们也会用到 ArcGIS Pro，可以实时看到制图的效果，并且渲染效果更好。

ArcPy我真的是又爱又恨啊！😑

我们来用 geopandas 读取刚刚创建的渔网，因为他在 gdb 数据库中，我们可以用read_file()去读取：

# 查看geopandas的版本
gpd.__version__ # 保证你的是0.12版本以上

# 读取渔网
gdf_fishnet = gpd.read_file(filename=gdb, layer=out_fcs)

# 读取渔网标注点
gdf_fishnet_label = gpd.read_file(filename=gdb, layer=out_label)

len(gdf_fishnet), len(gdf_fishnet_label)

OUT：

(12150, 12150)

12150 个数据，插个眼记录一下。

我们先新建一个 id 列，方便记录位置：

# 给渔网和标注点添加id列
gdf_fishnet['id'] = gdf_fishnet.index
gdf_fishnet_label['id'] = gdf_fishnet_label.index

空间连接的两种方式，一种是使用geopandas.GeoDataFrame.sjoin_nearest 同时搜寻和连接，一种是使用geopandas.GeoDataFrame.sjoin 直接进行空间连接。我们选择后者，这个函数很常用，建议访问 spatial-joins 官方帮助文档查看更多信息。

我们已经有方格网了，直接使用geopandas.GeoDataFrame.sjoin空间连接：

gdf_fishnet = gpd.sjoin(gdf_fishnet, gdf, how='left', predicate='intersects') # op参数在将来的版本中将被弃用，并建议使用predicate参数代替。

有两个参数需要注意：

• how=‘left’: 这指定了连接类型为"left"，这意味着结果中将包含gdf_fishnet中的所有几何图形。对于那些与gdf中的任何几何图形都没有交集的gdf_fishnet中的几何图形，连接的结果将为NaN。

• predicate=‘intersects’: 这是连接的谓词，指定了两个几何图形之间的空间关系。在这种情况下，它指的是"交集"。这意味着只要gdf_fishnet中的一个几何图形与gdf中的一个几何图形有任何形式的交集，它们就会被连接起来。

4）清洗网格并计算准确度

扔掉空网格数据：

gdf_fishnet = gdf_fishnet.dropna(subset=['id_right'])
gdf_fishnet.head()

计算预测的准确度

# 我们对diff列取绝对值，然后
gdf_fishnet['diff'].abs() # 取绝对值

正则化到 0-1

gdf_fishnet['diff'].abs() / gdf_fishnet['diff'].abs().max()

对 geodataframe 进行操作：

# 取值0-1 diff是差值为0预测准确，为1预测不准确，但是我们需要的是预测准确度，所以取1-diff
gdf_fishnet['accuracy'] = 1 - gdf_fishnet['diff'].abs() / gdf_fishnet['diff'].abs().max()
gdf_fishnet['accuracy']

5）聚合建筑准确度到网格上

两种方式：geopandas的dissolve和groupby。两者都可以用于在特定的列上执行聚合操作。dissolve的主要特点是它可以执行空间聚合。这意味着具有相同属性的邻近几何图形可以被合并成一个几何图形。例如，如果您有多个相邻的多边形，并且它们具有相同的属性值，dissolve可以将它们合并成一个大的多边形。

在使用sjoin()函数时返回的结果中，同一个小渔网会被和它相交的建筑多边形所相连，所以我们通过同一个渔网中的建筑物都拥有渔网要素的 id 来判断，这个 id 就是初始定义的 id 列，不过在使用sjoin()函数被重命名为"id_left"。

最后，我们选择一个聚合函数，根据求准确度accuracy的均值mean。

# 使用dissolve进行空间聚合并计算accuracy的平均值
accuracy_150m = gdf_fishnet.dissolve(by='id_left', aggfunc={"accuracy": "mean"})
accuracy_150m

使用dissolve进行空间聚合的同时，无法保留其他的列。因为dissolve的设计是为了合并那些具有相同键值的几何图形，并聚合其他列的值。我们还需要重新将 gdf_fishnet 中的部分列合并到dissolve结果中，为了使代码简介，我们使用join合并数据：

cols_to_keep = ['identificatie', 'true_label', 'pre_label', 'diff']

# 使用dissolve进行空间聚合并计算accuracy的平均值
accuracy_150m = gdf_fishnet.dissolve(by='id_left', aggfunc={"accuracy": "mean"})

# 直接使用join来连接其他需要的列
accuracy_150m = accuracy_150m.join(gdf_fishnet.set_index('id_left')[cols_to_keep], on='id_left')
accuracy_150m # join，默认使用索引作为连接键，更适合基于索引的连接。

6）将数据合并到渔网的点要素上

我们需要将 gdf_fishnet 中的预测标签和真实标签合并到 accuracy_150m 中。

# 我们需要的是点的准确度，所以将渔网的准确度赋值给渔网标注点 并且需要扔掉空值 我们使用merge
gdf_fishnet_point = gdf_fishnet_label.merge(accuracy_150m, left_on='id',right_on='id_left' , how='right')
gdf_fishnet_point

结果是包含两列几何信息的 gdf 对象，geometry_x是点，也就是渔网的点要素，geometry_y是渔网多边形。因为 merge 方法的左侧 gdf_fishnet_label，所以现在使用 gdf.geometry 获取的是gdf_fishnet_label的几何属性。

我们分别保存它们，面要素用于计算莫兰指数，点要素用于制作本课程的空间分布。

7）保存面要素

accuracy_polygon = gdf_fishnet_point.drop(columns=['id', 'geometry_x']).rename(columns={'geometry_y': 'geometry'})

accuracy_polygon_gdf = gpd.GeoDataFrame(accuracy_polygon, geometry='geometry').set_crs(32631)

# 保存一份geojson 用于空间自相关研究
accuracy_polygon_gdf.to_file(filename='../../data/output/accuracy_150m_polygon.geojson')

8）保存点要素

accuracy_points = gdf_fishnet_point.drop(columns=['id',  'geometry_y']).rename(
    columns={'geometry_x': 'geometry'})

gpd.GeoDataFrame(accuracy_points, geometry='geometry').set_crs(32631).to_file(filename=gdb, layer='accuracy_150m_points')

9）绘制

用 GIS 绘制出建筑年代预测结果在 150 米网格范围的准确度图：

上一篇：Part4.对建筑年代进行深度学习训练和预测(上)——《通过深度学习了解建筑年代和风格》

下一篇：Part5.对建筑风格进行深度学习训练和预测以及分析——《通过深度学习了解建筑年代和风格》

---

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

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

其他平台账号：

• 我的知乎
• 我的博客
• 我的 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)