191 lines
6.4 KiB
Python
191 lines
6.4 KiB
Python
import math
|
|
import os
|
|
import pdb
|
|
from rscder.gui.actions import ActionManager
|
|
from rscder.plugins.basic import BasicPlugin
|
|
from PyQt5.QtWidgets import QAction, QDialog, QHBoxLayout, QVBoxLayout, QPushButton
|
|
from PyQt5.QtCore import pyqtSignal
|
|
from PyQt5.QtGui import QIcon
|
|
from rscder.utils.icons import IconInstance
|
|
from rscder.utils.project import BasicLayer, Project, RasterLayer, SingleBandRasterLayer
|
|
from rscder.gui.layercombox import PairLayerCombox
|
|
from osgeo import gdal
|
|
from threading import Thread
|
|
import numpy as np
|
|
|
|
class MyDialog(QDialog):
|
|
|
|
def __init__(self, parent=None):
|
|
super().__init__(parent)
|
|
|
|
self.setWindowTitle('差分法')
|
|
self.setWindowIcon(IconInstance().LOGO)
|
|
|
|
# self.setFixedWidth(500)
|
|
|
|
self.layer_select = PairLayerCombox(self)
|
|
|
|
# self.number_input = QLineEdit(self)
|
|
|
|
self.ok_button = QPushButton('确定', self)
|
|
self.ok_button.setIcon(IconInstance().OK)
|
|
self.ok_button.clicked.connect(self.on_ok)
|
|
|
|
self.cancel_button = QPushButton('取消', self)
|
|
self.cancel_button.setIcon(IconInstance().CANCEL)
|
|
self.cancel_button.clicked.connect(self.on_cancel)
|
|
|
|
self.button_layout = QHBoxLayout()
|
|
self.button_layout.addWidget(self.ok_button)
|
|
self.button_layout.addWidget(self.cancel_button)
|
|
|
|
self.main_layout = QVBoxLayout()
|
|
self.main_layout.addWidget(self.layer_select)
|
|
|
|
self.main_layout.addLayout(self.button_layout)
|
|
self.setLayout(self.main_layout)
|
|
|
|
def on_ok(self):
|
|
if self.layer_select.layer1 is None or self.layer_select.layer2 is None:
|
|
self.reject()
|
|
return
|
|
self.accept()
|
|
|
|
def on_cancel(self):
|
|
self.reject()
|
|
|
|
|
|
|
|
class BasicMethod(BasicPlugin):
|
|
|
|
result_ok = pyqtSignal(dict)
|
|
|
|
@staticmethod
|
|
def info():
|
|
return {
|
|
'name': 'BasicMethod',
|
|
'description': 'BasicMethod',
|
|
'author': 'RSCDER',
|
|
'version': '1.0.0',
|
|
}
|
|
|
|
def set_action(self):
|
|
|
|
basic_diff_method = QAction('差分法')
|
|
ActionManager().unsupervised_menu.addAction(basic_diff_method)
|
|
basic_diff_method.setEnabled(False)
|
|
self.basic_diff_method = basic_diff_method
|
|
basic_diff_method.triggered.connect(self.basic_diff_alg)
|
|
|
|
def setup(self):
|
|
self.basic_diff_method.setEnabled(True)
|
|
|
|
def run_basic_diff_alg(self, layer1:RasterLayer, layer2:RasterLayer):
|
|
|
|
pth1 = layer1.path
|
|
pth2 = layer2.path
|
|
|
|
|
|
|
|
cell_size = layer1.layer_parent.cell_size
|
|
|
|
self.send_message.emit('开始计算差分法')
|
|
|
|
ds1 = gdal.Open(pth1)
|
|
ds2 = gdal.Open(pth2)
|
|
if not layer1.compare(layer2):
|
|
self.send_message.emit('两个图层的尺寸不同')
|
|
return
|
|
xsize = ds1.RasterXSize
|
|
ysize = ds1.RasterYSize
|
|
band = ds1.RasterCount
|
|
yblocks = ysize // cell_size[1]
|
|
|
|
driver = gdal.GetDriverByName('GTiff')
|
|
out_tif = os.path.join(Project().cmi_path, 'temp.tif')
|
|
out_ds = driver.Create(out_tif, xsize, ysize, 1, gdal.GDT_Float32)
|
|
out_ds.SetGeoTransform(ds1.GetGeoTransform())
|
|
out_ds.SetProjection(ds1.GetProjection())
|
|
max_diff = 0
|
|
min_diff = math.inf
|
|
for j in range(yblocks + 1):
|
|
|
|
self.message_send.emit(f'计算{j}/{yblocks}')
|
|
block_xy = (0, j * cell_size[1])
|
|
if block_xy[1] > ysize:
|
|
break
|
|
block_size = (xsize, cell_size[1])
|
|
if block_xy[1] + block_size[1] > ysize:
|
|
block_size = (xsize, ysize - block_xy[1])
|
|
block_data1 = ds1.ReadAsArray(*block_xy, *block_size)
|
|
block_data2 = ds2.ReadAsArray(*block_xy, *block_size)
|
|
|
|
if band == 1:
|
|
block_data1 = block_data1[None, ...]
|
|
block_data2 = block_data2[None, ...]
|
|
# pdb.set_trace()
|
|
block_diff = block_data1.sum(0) - block_data2.sum(0)
|
|
block_diff = block_diff.astype(np.float32)
|
|
block_diff = np.abs(block_diff)
|
|
|
|
min_diff = min(min_diff, block_diff[block_diff > 0].min())
|
|
max_diff = max(max_diff, block_diff.max())
|
|
out_ds.GetRasterBand(1).WriteArray(block_diff, *block_xy)
|
|
|
|
self.send_message.emit(f'完成{j}/{yblocks}')
|
|
|
|
out_ds.FlushCache()
|
|
del out_ds
|
|
self.send_message.emit('归一化概率中...')
|
|
temp_in_ds = gdal.Open(out_tif)
|
|
|
|
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer1.layer_parent.name, int(np.random.rand() * 100000)))
|
|
out_normal_ds = driver.Create(out_normal_tif, xsize, ysize, 1, gdal.GDT_Byte)
|
|
out_normal_ds.SetGeoTransform(ds1.GetGeoTransform())
|
|
out_normal_ds.SetProjection(ds1.GetProjection())
|
|
# hist = np.zeros(256, dtype=np.int32)
|
|
for j in range(yblocks+1):
|
|
block_xy = (0, j * cell_size[1])
|
|
if block_xy[1] > ysize:
|
|
break
|
|
block_size = (xsize, cell_size[1])
|
|
if block_xy[1] + block_size[1] > ysize:
|
|
block_size = (xsize, ysize - block_xy[1])
|
|
block_data = temp_in_ds.ReadAsArray(*block_xy, *block_size)
|
|
block_data = (block_data - min_diff) / (max_diff - min_diff) * 255
|
|
block_data = block_data.astype(np.uint8)
|
|
out_normal_ds.GetRasterBand(1).WriteArray(block_data, *block_xy)
|
|
# hist_t, _ = np.histogram(block_data, bins=256, range=(0, 256))
|
|
# hist += hist_t
|
|
# print(hist)
|
|
del temp_in_ds
|
|
try:
|
|
os.remove(out_tif)
|
|
except:
|
|
pass
|
|
|
|
raster_result_layer = SingleBandRasterLayer(None, True, out_normal_tif, BasicLayer.BOATH_VIEW)
|
|
|
|
# layer1.layer_parent.add_result_layer(point_result_lalyer)
|
|
layer1.layer_parent.add_result_layer(raster_result_layer)
|
|
|
|
# self.send_message.emit('完成计算变化表格')
|
|
|
|
self.send_message.emit('差分法计算完成')
|
|
|
|
def basic_diff_alg(self):
|
|
# layer_select =
|
|
|
|
layer_select = MyDialog(self.mainwindow)
|
|
if layer_select.exec_():
|
|
layer1 = layer_select.layer_select.layer1
|
|
layer2 = layer_select.layer_select.layer2
|
|
else:
|
|
return
|
|
# layer:PairLayer = list(self.project.layers.values())[0]
|
|
|
|
t = Thread(target=self.run_basic_diff_alg, args=(layer1, layer2))
|
|
t.start()
|
|
|
|
|
|
|