#coding:utf-8
# 参考《python地理空间分析指南》
#插值算法
#执行变化监测的步骤大致为:用gdal_array将图片导入numpy数据中——变化前后影像作差——图像分类着色——保存图片。具体的代码为
#实验结果中，绿色表示添加的元素，红色表示减少的元素。示例中采用的差值算法的处理结果虽然不够完美，但也能提供一种变化检测的处理思路。最终的实验结果为:
import numpy as np
import matplotlib.pyplot as plt
from osgeo import gdal, gdal_array
import cv2
import sys
# 加载两幅图像
# img1 = 'D:\\BaiduNetdiskWorkspace\\pythonAlgorTest\\24A.png'
# img2 = 'D:\\BaiduNetdiskWorkspace\\pythonAlgorTest\\24B.png'
# outimg='D:\BaiduNetdiskWorkspace\pythonAlgorTest\change.tif'

def main(img1,img2,outimg): 
 im_1 =img1
 im_2 = img2
 outurl = outimg
 print("输入路径1："+im_1)
 print("输入路径1："+im_2)
# 载入数组
 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, outurl, format='GTiff', prototype=im_2)
    out1 = None
 print("输出路径1："+outurl)
 return outurl

# # 查看差值图片
# ima = 'D:\BaiduNetdiskWorkspace\pythonAlgorTest\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()
if __name__ == '__main__':
    url = sys.argv[1]
    url2 =sys.argv[2]
    url3 =sys.argv[3]
    main(url,url2,url3)