Part4-2.对建筑年代的预测结果进行分析：绘制混淆矩阵、计算分类报告，绘制空间分布

# 重新加载
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

# 数据集的类别名称
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}

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)

# 只需要调整尺寸和转换为张量
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")

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])

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)

BATCH_SIZE = 128

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

166

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

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)

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)

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) # 可以打印具体结果

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

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

{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

'1653–1705'

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

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)   

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()

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]

from sklearn.metrics import confusion_matrix

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

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()

# 借助混淆矩阵计算各种分类指标(召回率、精确度和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)

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

# 设置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"]
) 

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

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

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

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

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

# 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

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

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

# 可视化你的热图 将其覆盖在输入图像上：
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()

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)

## 读取数据
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()

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

#我们需要更新标签以反映这个新顺序
# 这意味着我们需要将所有的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()

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

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()

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)

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']]

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

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

0.4870035261276262

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

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

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中不用重新读取

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 = 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')

# 导入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"

# 查看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)

(12150, 12150)

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

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

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

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

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

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

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

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，默认使用索引作为连接键，更适合基于索引的连接。

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

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')

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')