Python遥感变化监测

图像来源于LEVIR建筑物变化检测数据集（LEVIR-CD）

# 参考《python地理空间分析指南》
import numpy as np
import matplotlib.pyplot as plt
from osgeo import gdal, gdal_array
import cv2

# 加载两幅图像
im_1 = 'levir2/A.tif'
im_2 = 'levir2/B.tif'
# 载入数组
arr_1 = gdal_array.LoadFile(im_1).astype(np.int8)
arr_2 = gdal_array.LoadFile(im_2).astype(np.int8)
# 显示两幅图像
plt.figure(figsize=(20, 40))
plt.subplot(121);plt.imshow(cv2.merge((arr_1[0], arr_1[1], arr_1[2])).astype(np.uint8))
plt.subplot(122);plt.imshow(cv2.merge((arr_2[0], arr_2[1], arr_2[2])).astype(np.uint8))
plt.show()

# 进行变化检测
diff = arr_1 - arr_2
# 建立类别架构并将变化特征隔离
classes = np.histogram(diff, bins=5)[1]
# 用黑色遮罩遮住不明显的特征
lut = [[0,0,0], [0,0,0], [0,0,0], [255,0,0], [0,255,0], [100,50,255]]
# 类别初始值
start = 1
# 创建输出图片
rgb = np.zeros((3, diff.shape[1], diff.shape[2]), np.int8)
# 处理并分配颜色
for i in range(len(classes)):
    mask = np.logical_and(start <= diff, diff <= classes[i])
    masked = ((mask[0] == mask[1]) == mask[2])
    for j in range(len(lut[i])):
        rgb[j] = np.choose(masked, (rgb[j], lut[i][j]))
    strat = classes[i] + 1
# 保存
out = gdal_array.SaveArray(rgb, 'levir2/change.tif', format='GTiff', prototype=im_2)
out = None

# 查看差值图片
ima = 'levir2/change.tif'
arr = gdal_array.LoadFile(ima).astype(np.uint8)
plt.figure(figsize=(8, 8))
plt.imshow(cv2.merge((arr[0], arr[1], arr[2])))
plt.show()