wangtong/rscder-pub
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

1. 采用注册表机制注册算法

Browse Source 
2. 修改UI
3. 调整bug
This commit is contained in:
copper 2022-09-24 20:50:35 +08:00
parent 26dcc3a3e5
commit 06407ff7da
39 changed files with 1374 additions and 1182 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
ECD.py
View File

@ -1,13 +1,10 @@
import os
import sys
import argparse
sys.path.insert(0, os.path.join(os.path.dirname(__file__), 'libs'))
os.environ['PROJ_LIB'] = os.path.join(os.path.dirname(__file__), 'share/proj')
os.environ['GDAL_DATA'] = os.path.join(os.path.dirname(__file__), 'share')
os.environ['ECD_BASEDIR'] = os.path.dirname(__file__)
BASE_DIR = os.path.dirname(__file__)
import distutils
import distutils.version
from rscder import MulStart
import logging

1
plugins/In_one/__init__.py
View File

@ -1 +0,0 @@
from In_one.main import *

555
plugins/In_one/main.py
View File

@ -1,555 +0,0 @@
import os
from threading import Thread
import numpy as np
# from plugins.basic_change.main import MyDialog
from rscder.gui.actions import ActionManager
from rscder.plugins.basic import BasicPlugin
from PyQt5.QtWidgets import QAction, QDialog, QHBoxLayout, QVBoxLayout, QPushButton,QWidget,QLabel,QLineEdit,QPushButton,QComboBox,QDialogButtonBox
from PyQt5.QtGui import QPixmap
from rscder.gui.layercombox import PairLayerCombox
from rscder.utils.icons import IconInstance
from rscder.utils.geomath import geo2imageRC
from rscder.utils.project import Project, PairLayer,ResultPointLayer,MultiBandRasterLayer
from In_one.otsu import OTSU
from osgeo import gdal
from In_one import pic
import math
from skimage.filters import rank
from skimage.morphology import disk, rectangle
from In_one.scripts.UnsupervisedCD import LSTS,CVA,acd,aht,ocd,lhba,sh
def Meanfilter(x_size,y_size,layer:MultiBandRasterLayer):
x_size = int(x_size)
y_size = int(y_size)
pth = layer.path
if pth is None:
return
ds = gdal.Open(pth)
band_count = ds.RasterCount
out_path = os.path.join(Project().other_path, '{}_mean_filter.tif'.format(layer.name))
out_ds = gdal.GetDriverByName('GTiff').Create(out_path, ds.RasterXSize, ds.RasterYSize, band_count, ds.GetRasterBand(1).DataType)
out_ds.SetProjection(ds.GetProjection())
out_ds.SetGeoTransform(ds.GetGeoTransform())
for i in range(band_count):
band = ds.GetRasterBand(i+1)
data = band.ReadAsArray()
data = rank.mean(data, rectangle(y_size, x_size))
out_band = out_ds.GetRasterBand(i+1)
out_band.WriteArray(data)
out_ds.FlushCache()
del out_ds
del ds
return out_path
def basic_cd(pth1:str,pth2:str,layer_parent:PairLayer,send_message):
ds1:gdal.Dataset=gdal.Open(pth1)
ds2:gdal.Dataset=gdal.Open(pth2)
cell_size = layer_parent.cell_size
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
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)
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
out_ds.SetGeoTransform(geo)
out_ds.SetProjection(proj)
max_diff = 0
min_diff = math.inf
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
for j in range(yblocks + 1):#该改这里了
send_message.emit(f'计算{j}/{yblocks}')
block_xy1 = (start1x, start1y+j * cell_size[1])
block_xy2 = (start2x,start2y+j*cell_size[1])
block_xy=(0,j * cell_size[1])
if block_xy1[1] > end1y or block_xy2[1] > end2y:
break
block_size=(xsize, cell_size[1])
block_size1 = (xsize, cell_size[1])
block_size2 = (xsize,cell_size[1])
if block_xy[1] + block_size[1] > ysize:
block_size = (xsize, ysize - block_xy[1])
if block_xy1[1] + block_size1[1] > end1y:
block_size1 = (xsize,end1y - block_xy1[1])
if block_xy2[1] + block_size2[1] > end2y:
block_size2 = (xsize, end2y - block_xy2[1])
block_data1 = ds1.ReadAsArray(*block_xy1, *block_size1)
block_data2 = ds2.ReadAsArray(*block_xy2, *block_size2)
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)
send_message.emit(f'完成{j}/{yblocks}')
del ds2
del ds1
out_ds.FlushCache()
del out_ds
send_message.emit('归一化概率中...')
temp_in_ds = gdal.Open(out_tif)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(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(geo)
out_normal_ds.SetProjection(proj)
# 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
del out_normal_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('完成计算变化表格')
send_message.emit('差分法计算完成')
return out_normal_tif
def otsu(pth,name,send_message):
ds = gdal.Open(pth)
band = ds.GetRasterBand(1)
# band_count = ds.RasterCount
hist = np.zeros(256, dtype=np.int)
xsize = ds.RasterXSize
ysize = ds.RasterYSize
max_pixels = 1e7
max_rows = max_pixels // xsize
if max_rows < 1:
max_rows = 1
if max_rows > ysize:
max_rows = ysize
block_count = ysize // max_rows + 1
for i in range(block_count):
start_row = i * max_rows
end_row = min((i + 1) * max_rows, ysize)
block = band.ReadAsArray(0, start_row, xsize, end_row - start_row)
hist += np.histogram(block.flatten(), bins=256, range=(0, 255))[0]
hist = hist.astype(np.float32)
gap = OTSU(hist)
send_message.emit('阈值为：{}'.format(gap))
out_th = os.path.join(Project().bcdm_path, '{}_otsu_bcdm.tif'.format(name))
out_ds = gdal.GetDriverByName('GTiff').Create(out_th, xsize, ysize, 1, gdal.GDT_Byte)
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_ds.SetProjection(ds.GetProjection())
out_band = out_ds.GetRasterBand(1)
for i in range(block_count):
start_row = i * max_rows
end_row = min((i + 1) * max_rows, ysize)
block = band.ReadAsArray(0, start_row, xsize, end_row - start_row)
out_band.WriteArray(block > gap, 0, start_row)
out_band.FlushCache()
out_ds = None
ds = None
send_message.emit('OTSU阈值完成')
return out_th,gap
#otsu_layer = SingleBandRasterLayer(path = out_th, style_info={})
#layer.layer_parent.add_result_layer(otsu_layer)
def thresh(pth,gap,name,send_message):
ds = gdal.Open(pth)
band = ds.GetRasterBand(1)
# band_count = ds.RasterCount
xsize = ds.RasterXSize
ysize = ds.RasterYSize
max_pixels = 1e7
max_rows = max_pixels // xsize
if max_rows < 1:
max_rows = 1
if max_rows > ysize:
max_rows = ysize
block_count = ysize // max_rows + 1
# for i in range(block_count):
# start_row = i * max_rows
# end_row = min((i + 1) * max_rows, ysize)
# block = band.ReadAsArray(0, start_row, xsize, end_row - start_row)
# hist += np.histogram(block.flatten(), bins=256, range=(0, 255))[0]
# hist = hist.astype(np.float32)
send_message.emit('阈值为：{}'.format(gap))
out_th = os.path.join(Project().bcdm_path, '{}_thresh{}_bcdm.tif'.format(name,gap))
out_ds = gdal.GetDriverByName('GTiff').Create(out_th, xsize, ysize, 1, gdal.GDT_Byte)
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_ds.SetProjection(ds.GetProjection())
out_band = out_ds.GetRasterBand(1)
for i in range(block_count):
start_row = i * max_rows
end_row = min((i + 1) * max_rows, ysize)
block = band.ReadAsArray(0, start_row, xsize, end_row - start_row)
out_band.WriteArray(block > gap, 0, start_row)
out_band.FlushCache()
out_ds = None
ds = None
send_message.emit('自定义阈值分割完成')
return out_th
#otsu_layer = SingleBandRasterLayer(path = out_th, style_info={})
#layer.layer_parent.add_result_layer(otsu_layer)
def table_layer(pth,layer,name,send_message,dict):
send_message.emit('正在计算表格结果...')
cell_size = layer.layer_parent.cell_size
ds = gdal.Open(pth)
xsize = ds.RasterXSize
ysize = ds.RasterYSize
geo = ds.GetGeoTransform()
out_csv = os.path.join(Project().other_path, f'{name}_table_result.csv')
yblocks = ysize // cell_size[1] + 1
xblocks = xsize // cell_size[0] + 1
with open(out_csv, 'w') as f:
f.write('x,y,diff,status\n')
for j in range(yblocks):
block_xy = (0, j * cell_size[1])
block_size = (xsize, cell_size[1])
if block_xy[1] + block_size[1] > ysize:
block_size = (xsize, ysize - block_xy[1])
block_data = ds.ReadAsArray(*block_xy, *block_size)
for i in range(xblocks):
start_x = i * cell_size[0]
end_x = start_x + cell_size[0]
if end_x > xsize:
end_x = xsize
block_data_xy = block_data[:, start_x:end_x]
center_x = start_x + cell_size[0] // 2
center_y = j * cell_size[1] + cell_size[1] // 2
center_x = center_x * geo[1] + geo [0]
center_y = center_y * geo[5] + geo [3]
f.write(f'{center_x},{center_y},{block_data_xy.mean() * 100},{int(block_data_xy.mean() > 0.5)}\n')
result_layer = ResultPointLayer(out_csv, enable=True, proj=layer.proj, geo=layer.geo,result_path=dict)
# print(result_layer.result_path)
layer.layer_parent.add_result_layer(result_layer)
send_message.emit('计算完成')
class LockerButton(QPushButton):
def __init__(self,parent=None):
super(LockerButton,self).__init__(parent)
m_imageLabel = QLabel(self)
m_imageLabel.setFixedWidth(20)
m_imageLabel.setScaledContents(True)
m_imageLabel.setStyleSheet("QLabel{background-color:transparent;}")
m_textLabel = QLabel(self)
m_textLabel.setStyleSheet("QLabel{background-color:transparent;}")
self.m_imageLabel=m_imageLabel
self.m_textLabel=m_textLabel
self.hide_=1
mainLayout = QHBoxLayout()
mainLayout.addWidget(self.m_imageLabel)
mainLayout.addWidget(self.m_textLabel)
mainLayout.setContentsMargins(0,0,0,0)
mainLayout.setSpacing(0)
self.setLayout(mainLayout)
def SetImageLabel(self, pixmap:QPixmap):
self.m_imageLabel.setPixmap(pixmap)
def SetTextLabel(self, text):
self.m_textLabel.setText(text)
class selectCombox(QComboBox):
def __init__(self, parent,list:list,default='--') :
super(selectCombox,self).__init__(parent)
self.choose=None
self.list=list
self.default=default
self.addItem(default, None)
self.addItems(list)
self.currentIndexChanged.connect(self.on_change)
def on_change(self,index):
if index == 0:
self.choose=self.default
else:
self.choose=self.list[index-1]
# print(self.choose)
class AllInOne(QDialog):
def __init__(self, pre,cd,threshold,parent=None):
super(AllInOne, self).__init__(parent)
self.setWindowTitle('变化检测')
self.setWindowIcon(IconInstance().LOGO)
self.pre=pre#['均值滤波','test滤波']
self.cd=cd#['差分法','test法']
self.threshold=threshold#['OTSU']
self.initUI()
def initUI(self):
#图层
self.layer_combox = PairLayerCombox(self)
layerbox = QHBoxLayout()
layerbox.addWidget(self.layer_combox)
#预处理
filterWeight=QWidget(self)
filterlayout=QVBoxLayout()
filerButton =LockerButton(filterWeight)
filerButton.setObjectName("filerButton")
filerButton.SetTextLabel("预处理")
filerButton.SetImageLabel(QPixmap('plugins/In_one/pic/箭头_列表展开.png'))
filerButton.setStyleSheet("#filerButton{background-color:transparent;border:none;}"
"#filerButton:hover{background-color:rgba(195,195,195,0.4);border:none;}")
self.pre_select=selectCombox(self,self.pre)
x_size_input = QLineEdit(self)
x_size_input.setText('3')
y_size_input = QLineEdit(self)
y_size_input.setText('3')
size_label = QLabel(self)
size_label.setText('窗口大小:')
time_label = QLabel(self)
time_label.setText('X')
self.x_size_input = x_size_input
self.y_size_input = y_size_input
hlayout1 = QHBoxLayout()
hlayout1.addWidget(size_label)
hlayout1.addWidget(x_size_input)
hlayout1.addWidget(time_label)
hlayout1.addWidget(y_size_input)
vlayout = QVBoxLayout()
vlayout.addWidget(self.pre_select)
vlayout.addLayout(hlayout1)
filterWeight.setLayout(vlayout)
filterlayout.addWidget(filerButton)
filterlayout.addWidget(filterWeight)
#变化检测
changelayout=QVBoxLayout()
changeWeight=QWidget(self)
changeButton =LockerButton(changeWeight)
changeButton.setObjectName("changeButton")
changeButton.SetTextLabel("变化检测")
changeButton.SetImageLabel(QPixmap('plugins/In_one/pic/箭头_列表展开.png'))
changeButton.setStyleSheet("#changeButton{background-color:transparent;border:none;}"
"#changeButton:hover{background-color:rgba(195,195,195,0.4);border:none;}")
changeWeightlayout=QVBoxLayout()
self.cd_select=selectCombox(self,self.cd)
changeWeightlayout.addWidget(self.cd_select)
changeWeight.setLayout(changeWeightlayout)
changelayout.addWidget(changeButton)
changelayout.addWidget(changeWeight)
#阈值处理
thresholdlayout=QVBoxLayout()
thresholdWeight=QWidget(self)
thresholdButton =LockerButton(thresholdWeight)
thresholdButton.setObjectName("thresholdButton")
thresholdButton.SetTextLabel("阈值处理")
thresholdButton.SetImageLabel(QPixmap('plugins/In_one/pic/箭头_列表展开.png'))
thresholdButton.setStyleSheet("#thresholdButton{background-color:transparent;border:none;}"
"#thresholdButton:hover{background-color:rgba(195,195,195,0.4);border:none;}")
self.threshold_select=selectCombox(self,self.threshold,default='手动阈值')
self.threshold_input=QLineEdit(self)
self.threshold_input.setText('0.5')
self.threshold_select.currentIndexChanged.connect(lambda index:self.hide_(self.threshold_input,index==0))
thresholdWeightlayout=QVBoxLayout()
thresholdWeightlayout.addWidget(self.threshold_select)
thresholdWeightlayout.addWidget(self.threshold_input)
thresholdWeight.setLayout(thresholdWeightlayout)
thresholdlayout.addWidget(thresholdButton)
thresholdlayout.addWidget(thresholdWeight)
#确认
self.ok_button = QPushButton('确定', self)
self.ok_button.setIcon(IconInstance().OK)
self.ok_button.clicked.connect(self.accept)
self.ok_button.setDefault(True)
self.cancel_button = QPushButton('取消', self)
self.cancel_button.setIcon(IconInstance().CANCEL)
self.cancel_button.clicked.connect(self.reject)
self.cancel_button.setDefault(False)
buttonbox=QDialogButtonBox(self)
buttonbox.addButton(self.ok_button,QDialogButtonBox.NoRole)
buttonbox.addButton(self.cancel_button,QDialogButtonBox.NoRole)
buttonbox.setCenterButtons(True)
#buttonbox.setContentsMargins(,)
totalvlayout=QVBoxLayout()
totalvlayout.addLayout(layerbox)
totalvlayout.addLayout(filterlayout)
totalvlayout.addLayout(changelayout)
totalvlayout.addLayout(thresholdlayout)
totalvlayout.addWidget(buttonbox)
totalvlayout.addStretch()
self.setLayout(totalvlayout)
filerButton.clicked.connect(lambda: self.hide(filerButton,filterWeight))
changeButton.clicked.connect(lambda: self.hide(changeButton,changeWeight))
thresholdButton.clicked.connect(lambda: self.hide(thresholdButton,thresholdWeight))
def hide(self,button:LockerButton,weight:QWidget):
if ((button.hide_)%2)==1:
weight.setVisible(False)
button.SetImageLabel(QPixmap('plugins/In_one/pic/箭头_列表向右.png'))
else:
weight.setVisible(True)
button.SetImageLabel(QPixmap('plugins/In_one/pic/箭头_列表展开.png'))
button.hide_=(button.hide_)%2+1
def hide_(self,widget:QWidget,h:bool):
if h:
widget.setVisible(True)
else:
widget.setVisible(False)
def hideWidget(self,widget:QWidget):
if widget.isVisible:
widget.setVisible(False)
else:
widget.setVisible(True)
@property
def param(self):
class param(object):
pass
p=param()
p.pre=self.pre
p.cd=self.cd
p.threshold=self.threshold
p.layer_combox=self.layer_combox
p.pre_select=self.pre_select
p.x_size_input=self.x_size_input
p.y_size_input=self.y_size_input
p.threshold_select=self.threshold_select
p.threshold_input=self.threshold_input
p.cd_select=self.cd_select
return p
class InOnePlugin(BasicPlugin):
pre={"均值滤波":Meanfilter}#可添加其他方法
cd={'差分法':basic_cd,'LSTS':LSTS,'CVA':CVA,'ACD':acd,'AHT':aht,'OCD':ocd,'LHBA':lhba,'SH':sh}#可添加其他方法
threshold={'OTSU阈值':otsu}#可添加其他方法
@staticmethod
def info():
return {
'name': 'AllinOne',
'description': 'AllinOne',
'author': 'RSCDER',
'version': '1.0.0',
}
def set_action(self):
basic_diff_method_in_one = QAction(IconInstance().UNSUPERVISED, '&无监督变化检测')
ActionManager().change_detection_menu.addAction(basic_diff_method_in_one)
# ActionManager().menubar.addAction(basic_diff_method_in_one)
self.basic_diff_method_in_one = basic_diff_method_in_one
basic_diff_method_in_one.triggered.connect(self.run)
def run(self):
myDialog=AllInOne(list(self.pre.keys()),list(self.cd.keys()),list(self.threshold.keys()),self.mainwindow)
myDialog.show()
if myDialog.exec_()==QDialog.Accepted:
w=myDialog.param
t=Thread(target=self.run_alg,args=(w,))
t.start()
def run_alg(self,w:AllInOne):
dict={}
layer1=w.layer_combox.layer1
pth1 = w.layer_combox.layer1.path
pth2 = w.layer_combox.layer2.path
name=layer1.layer_parent.name
# 预处理
# 若添加的预处理函数接口相同，则无需判断是哪种方法
# if w.pre_select.choose==self.pre.keys()[0]:
# pass
# el
preKey=w.pre_select.choose
pth1=self.pre[preKey](w.x_size_input.text(),w.y_size_input.text(),w.layer_combox.layer1)
self.send_message.emit('{}图像{}'.format(preKey,w.layer_combox.layer1.name))
pth2=self.pre[preKey](w.x_size_input.text(),w.y_size_input.text(),w.layer_combox.layer2)
self.send_message.emit('{}图像{}'.format(preKey,w.layer_combox.layer2.name))
name=name+'_'+preKey
dict['预处理']=[preKey,'|'.format(pth1,pth2)]
cdpth=None
#变化检测
# if w.cd_select.choose==self.cd[0]:
# if w.cd_select.choose=='ACD':
# cdKey='ACD'
# self.send_message.emit('ACD计算中...')
# cdpth=os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(w.layer_combox.layer1.layer_parent.name, int(np.random.rand() * 100000)))
# print(pth1,pth2,cdpth)
# ACD(pth1,pth2,cdpth)
# else:
# pass
cdKey=w.cd_select.choose
cdpth=self.cd[cdKey](pth1,pth2,w.layer_combox.layer1.layer_parent,self.send_message)
name += '_'+cdKey
dict['变化检测算法']=[cdKey,cdpth]
#阈值处理
#例如手动阈值和otsu参数不同，则要做区分
thpth=None
if w.threshold_select.choose=='手动阈值':
thpth=thresh(cdpth,float(w.threshold_input.text()),w.layer_combox.layer1.layer_parent.name,self.send_message)
dict['后处理']=['手动阈值',[float(w.threshold_input.text())],thpth]
else:
thpth,gap=otsu(cdpth,w.layer_combox.layer1.layer_parent.name,self.send_message)
name+='_otsu'
dict['后处理']=['OTSU阈值',gap,cdpth]
table_layer(thpth,layer1,name,self.send_message,dict)

42
plugins/In_one/otsu.py
View File

@ -1,42 +0,0 @@
import numpy as np
def OTSU(hist):
u1=0.0#背景像素的平均灰度值
u2=0.0#前景像素的平均灰度值
th=0.0
#总的像素数目
PixSum= np.sum(hist)
#各灰度值所占总像素数的比例
PixRate=hist / PixSum
#统计各个灰度值的像素个数
Max_var = 0
#确定最大类间方差对应的阈值
GrayScale = len(hist)
for i in range(1,len(hist)):#从1开始是为了避免w1为0.
u1_tem=0.0
u2_tem=0.0
#背景像素的比列
w1=np.sum(PixRate[:i])
#前景像素的比例
w2=1.0-w1
if w1==0 or w2==0:
pass
else:#背景像素的平均灰度值
for m in range(i):
u1_tem=u1_tem+PixRate[m]*m
u1 = u1_tem * 1.0 / w1
#前景像素的平均灰度值
for n in range(i,GrayScale):
u2_tem = u2_tem + PixRate[n]*n
u2 = u2_tem / w2
#print(u1)
#类间方差公式：G=w1*w2*(u1-u2)**2
tem_var=w1*w2*np.power((u1-u2),2)
#print(tem_var)
#判断当前类间方差是否为最大值。
if Max_var<tem_var:
Max_var=tem_var#深拷贝，Max_var与tem_var占用不同的内存空间。
th=i
return th

2
plugins/In_one/pic/箭头 下.svg
View File

@ -1,2 +0,0 @@
<?xml version="1.0" standalone="no"?><!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd"><svg t="1655105353713" class="icon" viewBox="0 0 1024 1024" version="1.1" xmlns="http://www.w3.org/2000/svg" p-id="2292" xmlns:xlink="http://www.w3.org/1999/xlink" width="128" height="128"><defs><style type="text/css">@font-face { font-family: feedback-iconfont; src: url("//at.alicdn.com/t/font_1031158_u69w8yhxdu.woff2?t=1630033759944") format("woff2"), url("//at.alicdn.com/t/font_1031158_u69w8yhxdu.woff?t=1630033759944") format("woff"), url("//at.alicdn.com/t/font_1031158_u69w8yhxdu.ttf?t=1630033759944") format("truetype"); }
</style></defs><path d="M517.688889 796.444444c-45.511111 0-85.333333-17.066667-119.466667-51.2L73.955556 381.155556c-22.755556-22.755556-17.066667-56.888889 5.688888-79.644445 22.755556-22.755556 56.888889-17.066667 79.644445 5.688889l329.955555 364.088889c5.688889 5.688889 17.066667 11.377778 28.444445 11.377778s22.755556-5.688889 34.133333-17.066667l312.888889-364.088889c22.755556-22.755556 56.888889-28.444444 79.644445-5.688889 22.755556 22.755556 28.444444 56.888889 5.688888 79.644445L637.155556 739.555556c-28.444444 39.822222-68.266667 56.888889-119.466667 56.888888 5.688889 0 0 0 0 0z" p-id="2293"></path></svg>
Side by Side

Before

Width:  |  Height:  |  Size: 1.3 KiB

2
plugins/In_one/pic/箭头 右.svg
View File

@ -1,2 +0,0 @@
<?xml version="1.0" standalone="no"?><!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd"><svg t="1655105348118" class="icon" viewBox="0 0 1024 1024" version="1.1" xmlns="http://www.w3.org/2000/svg" p-id="2137" xmlns:xlink="http://www.w3.org/1999/xlink" width="128" height="128"><defs><style type="text/css">@font-face { font-family: feedback-iconfont; src: url("//at.alicdn.com/t/font_1031158_u69w8yhxdu.woff2?t=1630033759944") format("woff2"), url("//at.alicdn.com/t/font_1031158_u69w8yhxdu.woff?t=1630033759944") format("woff"), url("//at.alicdn.com/t/font_1031158_u69w8yhxdu.ttf?t=1630033759944") format("truetype"); }
</style></defs><path d="M312.888889 995.555556c-17.066667 0-28.444444-5.688889-39.822222-17.066667-22.755556-22.755556-17.066667-56.888889 5.688889-79.644445l364.088888-329.955555c11.377778-11.377778 17.066667-22.755556 17.066667-34.133333 0-11.377778-5.688889-22.755556-17.066667-34.133334L273.066667 187.733333c-22.755556-22.755556-28.444444-56.888889-5.688889-79.644444 22.755556-22.755556 56.888889-28.444444 79.644444-5.688889l364.088889 312.888889c34.133333 28.444444 56.888889 73.955556 56.888889 119.466667s-17.066667 85.333333-51.2 119.466666l-364.088889 329.955556c-11.377778 5.688889-28.444444 11.377778-39.822222 11.377778z" p-id="2138"></path></svg>
Side by Side

Before

Width:  |  Height:  |  Size: 1.3 KiB

BIN
plugins/In_one/pic/箭头_列表向右.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 1.2 KiB

BIN
plugins/In_one/pic/箭头_列表展开.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 1.8 KiB

553
plugins/In_one/scripts/UnsupervisedCD.py
View File

@ -1,553 +0,0 @@
from osgeo import gdal
import math,os
import time
from sklearn.cluster import k_means
from rscder.utils.geomath import geo2imageRC, imageRC2geo
from rscder.utils.project import Project, PairLayer
import numpy as np
from .ACD import ACD
from .AHT import AHT
from .OCD import OCD
from .LHBA import LHBA
from .SH import SH
def warp(file,ds:gdal.Dataset,srcWin=[0,0,0,0]):
driver = gdal.GetDriverByName('GTiff')
xsize=ds.RasterXSize
ysize=ds.RasterYSize
geo=ds.GetGeoTransform()
orj=ds.GetProjection()
band=ds.RasterCount
if os.path.exists(file):
os.remove(file)
out_ds:gdal.Dataset=driver.Create(file, xsize, ysize, band, gdal.GDT_Byte)
out_ds.SetGeoTransform(geo)
out_ds.SetProjection(orj)
for b in range(1,band+1):
out_ds.GetRasterBand(b).WriteArray(ds.ReadAsArray(*srcWin,band_list=[b]),*(0,0))
del out_ds
def basic_cd(pth1:str,pth2:str,layer_parent:PairLayer,send_message):
ds1:gdal.Dataset=gdal.Open(pth1)
ds2:gdal.Dataset=gdal.Open(pth2)
cell_size = layer_parent.cell_size
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
band = ds1.RasterCount
yblocks = ysize // cell_size[1]
driver = gdal.GetDriverByName('GTiff')
out_tif = os.path.join(Project().other_path, 'temp.tif')
out_ds = driver.Create(out_tif, xsize, ysize, 1, gdal.GDT_Float32)
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
out_ds.SetGeoTransform(geo)
out_ds.SetProjection(proj)
max_diff = 0
min_diff = math.inf
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
for j in range(yblocks + 1):#该改这里了
send_message.emit(f'计算{j}/{yblocks}')
block_xy1 = (start1x, start1y+j * cell_size[1])
block_xy2 = (start2x,start2y+j*cell_size[1])
block_xy=(0,j * cell_size[1])
if block_xy1[1] > end1y or block_xy2[1] > end2y:
break
block_size=(xsize, cell_size[1])
block_size1 = (xsize, cell_size[1])
block_size2 = (xsize,cell_size[1])
if block_xy[1] + block_size[1] > ysize:
block_size = (xsize, ysize - block_xy[1])
if block_xy1[1] + block_size1[1] > end1y:
block_size1 = (xsize,end1y - block_xy1[1])
if block_xy2[1] + block_size2[1] > end2y:
block_size2 = (xsize, end2y - block_xy2[1])
block_data1 = ds1.ReadAsArray(*block_xy1, *block_size1)
block_data2 = ds2.ReadAsArray(*block_xy2, *block_size2)
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)
send_message.emit(f'完成{j}/{yblocks}')
del ds2
del ds1
out_ds.FlushCache()
del out_ds
send_message.emit('归一化概率中...')
temp_in_ds = gdal.Open(out_tif)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(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(geo)
out_normal_ds.SetProjection(proj)
# 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
del out_normal_ds
try:
os.remove(out_tif)
except:
pass
send_message.emit('差分法计算完成')
return out_normal_tif
def LSTS(pth1:str,pth2:str,layer_parent:PairLayer,send_message,n=5,w_size=(3,3)):
ds1:gdal.Dataset=gdal.Open(pth1)
ds2:gdal.Dataset=gdal.Open(pth2)
cell_size = layer_parent.cell_size
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
band = ds1.RasterCount
yblocks = ysize // cell_size[1]
driver = gdal.GetDriverByName('GTiff')
out_tif = os.path.join(Project().other_path, 'temp.tif')
out_ds = driver.Create(out_tif, xsize, ysize, 1, gdal.GDT_Float32)
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
out_ds.SetGeoTransform(geo)
out_ds.SetProjection(proj)
pixnum=w_size[0]*w_size[1]
print('pixnum',pixnum)
max_diff = 0
min_diff = math.inf
win_h=w_size[0]//2 #half hight of window
win_w=w_size[1]//2 #half width of window
a=[[(i+1)**j for j in range(n+1)] for i in range(pixnum)]
A=np.array(a).astype(np.float64)#
k_=np.array(range(1,n+1))
df1=np.zeros(pixnum).astype(np.float64)
df2=np.zeros(pixnum).astype(np.float64)
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
for j in range(yblocks + 1):
send_message.emit(f'计算{j}/{yblocks}')
block_xy1 = (start1x, start1y+j * cell_size[1])
block_xy2 = (start2x,start2y+j*cell_size[1])
block_xy=(0,j * cell_size[1])
if block_xy1[1] > end1y or block_xy2[1] > end2y:
break
block_size=(xsize, cell_size[1])
block_size1 = (xsize, cell_size[1])
block_size2 = (xsize,cell_size[1])
if block_xy[1] + block_size[1] > ysize:
block_size = (xsize, ysize - block_xy[1])
if block_xy1[1] + block_size1[1] > end1y:
block_size1 = (xsize,end1y - block_xy1[1])
if block_xy2[1] + block_size2[1] > end2y:
block_size2 = (xsize, end2y - block_xy2[1])
block_data1 = ds1.ReadAsArray(*block_xy1, *block_size1)
block_data2 = ds2.ReadAsArray(*block_xy2, *block_size2)
if band == 1:
block_data1 = block_data1[None, ...]
block_data2 = block_data2[None, ...]
# pdb.set_trace()
else:
block_data1=np.mean(block_data1,0)
block_data2=np.mean(block_data2,0)
block_diff=np.zeros(block_data1.shape).astype(np.float64)
for i in range(win_h,block_size1[1]-win_h):
for j_ in range(win_w,block_size1[0]-win_w):
pix=0
#get b
# b1=block_data[i+win_h:i+win_h] c in range(j_-win_w,j_+win_w+1)
b1=block_data1[i-win_h:i+win_h+1,j_-win_w:j_+win_w+1]
b2=block_data2[i-win_h:i+win_h+1,j_-win_w:j_+win_w+1]
b1=[b if (r+1)//2 else b[::-1] for r,b in enumerate(b1)]
b2=[b if (r+1)//2 else b[::-1] for r,b in enumerate(b2)]
b1=np.expand_dims(np.concatenate(b1,0),1)
b2=np.expand_dims(np.concatenate(b2,0),1)
#end get b
x1=np.squeeze(np.linalg.pinv(A).dot(b1))
x2=np.squeeze(np.linalg.pinv(A).dot(b2))
#df
k_=range(1,n+1)
for pix in range(1,pixnum+1):
df1[pix-1]=x1[1:n+1].dot(np.array([k*(pix**(k-1)) for k in k_]))
df2[pix-1]=x2[1:n+1].dot(np.array([k*(pix**(k-1)) for k in k_]))
#distance 欧式距离
block_diff[i][j_]=np.dot(df1-df2,df1-df2)**0.5
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)
send_message.emit(f'完成{j}/{yblocks}')
del ds2
del ds1
out_ds.FlushCache()
del out_ds
send_message.emit('归一化概率中...')
temp_in_ds = gdal.Open(out_tif)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(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(geo)
out_normal_ds.SetProjection(proj)
# 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
del out_normal_ds
try:
os.remove(out_tif)
except:
pass
send_message.emit('LSTS法计算完成')
return out_normal_tif
def CVA(pth1:str,pth2:str,layer_parent:PairLayer,send_message):
ds1:gdal.Dataset=gdal.Open(pth1)
ds2:gdal.Dataset=gdal.Open(pth2)
cell_size = layer_parent.cell_size
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
band = ds1.RasterCount
yblocks = ysize // cell_size[1]
driver = gdal.GetDriverByName('GTiff')
out_tif = os.path.join(Project().other_path, 'temp.tif')
out_ds = driver.Create(out_tif, xsize, ysize, 1, gdal.GDT_Float32)
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
out_ds.SetGeoTransform(geo)
out_ds.SetProjection(proj)
max_diff = 0
min_diff = math.inf
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
for j in range(yblocks + 1):
send_message.emit(f'计算{j}/{yblocks}')
block_xy1 = (start1x, start1y+j * cell_size[1])
block_xy2 = (start2x,start2y+j*cell_size[1])
block_xy=(0,j * cell_size[1])
if block_xy1[1] > end1y or block_xy2[1] > end2y:
break
block_size=(xsize, cell_size[1])
block_size1 = (xsize, cell_size[1])
block_size2 = (xsize,cell_size[1])
if block_xy[1] + block_size[1] > ysize:
block_size = (xsize, ysize - block_xy[1])
if block_xy1[1] + block_size1[1] > end1y:
block_size1 = (xsize,end1y - block_xy1[1])
if block_xy2[1] + block_size2[1] > end2y:
block_size2 = (xsize, end2y - block_xy2[1])
block_data1 = ds1.ReadAsArray(*block_xy1, *block_size1)
block_data2 = ds2.ReadAsArray(*block_xy2, *block_size2)
if band == 1:
block_data1 = block_data1[None, ...]
block_data2 = block_data2[None, ...]
# pdb.set_trace()
block_diff=np.sum((block_data1-block_data2)**2,0)**0.5
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)
send_message.emit(f'完成{j}/{yblocks}')
del ds2
del ds1
out_ds.FlushCache()
del out_ds
send_message.emit('归一化概率中...')
temp_in_ds = gdal.Open(out_tif)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(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(geo)
out_normal_ds.SetProjection(proj)
# 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
del out_normal_ds
try:
os.remove(out_tif)
except:
pass
send_message.emit('欧式距离计算完成')
return out_normal_tif
def acd(pth1:str,pth2:str,layer_parent:PairLayer,send_message):
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
#提取公共部分
send_message.emit('提取重叠区域数据.....')
ds2:gdal.Dataset=gdal.Open(pth2)
temp_tif2 = os.path.join(Project().other_path,'temp2.tif')
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize])
del ds2
send_message.emit('图像二提取完成')
ds1:gdal.Dataset=gdal.Open(pth1)
temp_tif1 = os.path.join(Project().other_path, 'temp1.tif')
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize])
del ds1
send_message.emit('图像一提取完成')
#运算
send_message.emit('开始ACD计算.....')
time.sleep(0.1)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000)))
ACD(temp_tif1,temp_tif2,out_normal_tif)
#添加投影
send_message.emit('录入投影信息.....')
time.sleep(0.1)
ds=gdal.Open(out_normal_tif,1)
ds.SetGeoTransform(geo)
ds.SetProjection(proj)
del ds
return out_normal_tif
def aht(pth1:str,pth2:str,layer_parent:PairLayer,send_message):
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
#提取公共部分
send_message.emit('提取重叠区域数据.....')
ds2:gdal.Dataset=gdal.Open(pth2)
temp_tif2 = os.path.join(Project().other_path,'temp2.tif')
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize])
del ds2
send_message.emit('图像二提取完成')
ds1:gdal.Dataset=gdal.Open(pth1)
temp_tif1 = os.path.join(Project().other_path, 'temp1.tif')
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize])
del ds1
send_message.emit('图像一提取完成')
#运算
send_message.emit('开始AHT计算.....')
time.sleep(0.1)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000)))
AHT(temp_tif1,temp_tif2,out_normal_tif)
#添加投影
send_message.emit('录入投影信息.....')
time.sleep(0.1)
ds=gdal.Open(out_normal_tif,1)
ds.SetGeoTransform(geo)
ds.SetProjection(proj)
del ds
return out_normal_tif
def ocd(pth1:str,pth2:str,layer_parent:PairLayer,send_message):
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
#提取公共部分
send_message.emit('提取重叠区域数据.....')
ds2:gdal.Dataset=gdal.Open(pth2)
temp_tif2 = os.path.join(Project().other_path,'temp2.tif')
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize])
del ds2
send_message.emit('图像二提取完成')
ds1:gdal.Dataset=gdal.Open(pth1)
temp_tif1 = os.path.join(Project().other_path, 'temp1.tif')
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize])
del ds1
send_message.emit('图像一提取完成')
#运算
send_message.emit('开始OCD计算.....')
time.sleep(0.1)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000)))
OCD(temp_tif1,temp_tif2,out_normal_tif,Project().other_path)
#添加投影
send_message.emit('录入投影信息.....')
time.sleep(0.1)
ds=gdal.Open(out_normal_tif,1)
ds.SetGeoTransform(geo)
ds.SetProjection(proj)
del ds
return out_normal_tif
def lhba(pth1:str,pth2:str,layer_parent:PairLayer,send_message):
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
#提取公共部分
send_message.emit('提取重叠区域数据.....')
ds2:gdal.Dataset=gdal.Open(pth2)
temp_tif2 = os.path.join(Project().other_path,'temp2.tif')
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize])
del ds2
send_message.emit('图像二提取完成')
ds1:gdal.Dataset=gdal.Open(pth1)
temp_tif1 = os.path.join(Project().other_path, 'temp1.tif')
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize])
del ds1
send_message.emit('图像一提取完成')
#运算
send_message.emit('开始LHBA计算.....')
time.sleep(0.1)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000)))
LHBA(temp_tif1,temp_tif2,out_normal_tif)
#添加投影
send_message.emit('录入投影信息.....')
time.sleep(0.1)
ds=gdal.Open(out_normal_tif,1)
ds.SetGeoTransform(geo)
ds.SetProjection(proj)
del ds
return out_normal_tif
def sh(pth1:str,pth2:str,layer_parent:PairLayer,send_message):
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
#提取公共部分
send_message.emit('提取重叠区域数据.....')
ds2:gdal.Dataset=gdal.Open(pth2)
temp_tif2 = os.path.join(Project().other_path,'temp2.tif')
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize])
del ds2
send_message.emit('图像二提取完成')
ds1:gdal.Dataset=gdal.Open(pth1)
temp_tif1 = os.path.join(Project().other_path, 'temp1.tif')
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize])
del ds1
send_message.emit('图像一提取完成')
#运算
send_message.emit('开始SH计算.....')
time.sleep(0.1)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000)))
SH(temp_tif1,temp_tif2,out_normal_tif)
#添加投影
send_message.emit('录入投影信息.....')
time.sleep(0.1)
ds=gdal.Open(out_normal_tif,1)
ds.SetGeoTransform(geo)
ds.SetProjection(proj)
del ds
return out_normal_tif

BIN
plugins/In_one/scripts/_ACD.pyd
View File

Binary file not shown.

BIN
plugins/In_one/scripts/_AHT.pyd
View File

Binary file not shown.

BIN
plugins/In_one/scripts/_LHBA.pyd
View File

Binary file not shown.

BIN
plugins/In_one/scripts/_OCD.pyd
View File

Binary file not shown.

BIN
plugins/In_one/scripts/_SH.pyd
View File

Binary file not shown.

BIN
plugins/In_one/scripts/opencv_ffmpeg3415_64.dll
View File

Binary file not shown.

BIN
plugins/In_one/scripts/opencv_world3415.dll
View File

Binary file not shown.

3
plugins/In_one/test.py
View File

@ -1,3 +0,0 @@
from In_one.scripts.USCD import ACD
from osgeo import gdal

5
plugins/filter_collection/__init__.py Normal file
View File

@ -0,0 +1,5 @@
from misc import Register
FILTER = Register('filter')
from filter_collection.main import *

98
plugins/some_filter/main.py → plugins/filter_collection/main.py
View File

@ -1,3 +1,4 @@
from datetime import datetime
import os
from threading import Thread
from PyQt5.QtWidgets import QDialog, QAction
@ -10,19 +11,98 @@ from rscder.plugins.basic import BasicPlugin
from rscder.gui.layercombox import RasterLayerCombox
from osgeo import gdal, gdal_array
from skimage.filters import rank
from skimage.morphology import disk, rectangle
from skimage.morphology import rectangle
from filter_collection import FILTER
from misc import AlgFrontend
@FILTER.register
class MainFilter(AlgFrontend):
@staticmethod
def get_name():
return '均值滤波'
@staticmethod
def get_widget(parent=None):
widget = QtWidgets.QWidget(parent)
x_size_input = QtWidgets.QLineEdit(widget)
x_size_input.setText('3')
x_size_input.setValidator(QtGui.QIntValidator())
x_size_input.setObjectName('xinput')
y_size_input = QtWidgets.QLineEdit(widget)
y_size_input.setValidator(QtGui.QIntValidator())
y_size_input.setObjectName('yinput')
y_size_input.setText('3')
size_label = QtWidgets.QLabel(widget)
size_label.setText('窗口大小:')
time_label = QtWidgets.QLabel(widget)
time_label.setText('X')
hlayout1 = QtWidgets.QHBoxLayout()
hlayout1.addWidget(size_label)
hlayout1.addWidget(x_size_input)
hlayout1.addWidget(time_label)
hlayout1.addWidget(y_size_input)
widget.setLayout(hlayout1)
return widget
@staticmethod
def get_params(widget:QtWidgets.QWidget=None):
if widget is None:
return dict(x_size=3, y_size=3)
x_input = widget.findChild(QtWidgets.QLineEdit, 'xinput')
y_input = widget.findChild(QtWidgets.QLineEdit, 'yinput')
if x_input is None or y_input is None:
return dict(x_size=3, y_size=3)
x_size = int(x_input.text())
y_size = int(y_input.text())
return dict(x_size=x_size, y_size=y_size)
@staticmethod
def run_alg(pth, x_size, y_size, *args, **kargs):
x_size = int(x_size)
y_size = int(y_size)
# pth = layer.path
if pth is None:
return
ds = gdal.Open(pth)
band_count = ds.RasterCount
out_path = os.path.join(Project().other_path, 'mean_filter_{}.tif'.format(int(datetime.now().timestamp() * 1000)))
out_ds = gdal.GetDriverByName('GTiff').Create(out_path, ds.RasterXSize, ds.RasterYSize, band_count, ds.GetRasterBand(1).DataType)
out_ds.SetProjection(ds.GetProjection())
out_ds.SetGeoTransform(ds.GetGeoTransform())
for i in range(band_count):
band = ds.GetRasterBand(i+1)
data = band.ReadAsArray()
data = rank.mean(data, rectangle(y_size, x_size))
out_band = out_ds.GetRasterBand(i+1)
out_band.WriteArray(data)
out_ds.FlushCache()
del out_ds
del ds
return out_path
class FilterSetting(QDialog):
def __init__(self, parent=None):
super(FilterSetting, self).__init__(parent)
self.setWindowTitle('滤波设置')
# self.setWindowFlags(Qt.WindowStaysOnTopHint)
# self.setFixedSize(300, 200)
# self.setStyleSheet("QDialog{background-color:rgb(255,255,255);}")
self.setWindowIcon(IconInstance().FILTER)
# self.setWindowIconText('Filter Setting')
# self.setWindowModality(Qt.ApplicationModal)
self.initUI()
# self.show()
def initUI(self):
self.layer_combox = RasterLayerCombox(self)
@ -78,10 +158,10 @@ class MainPlugin(BasicPlugin):
@staticmethod
def info():
return {
'name': 'mean_filter',
'name': 'FilterCollection',
'author': 'rscder',
'version': '0.0.1',
'description': 'Mean Filter'
'description': 'Filter Collections'
}
def set_action(self):

3
plugins/misc/__init__.py Normal file
View File

@ -0,0 +1,3 @@
from misc.main import *
from misc.utils import *
from misc.table_layer import *

80
plugins/misc/main.py Normal file
View File

@ -0,0 +1,80 @@
from threading import Thread
from plugins.misc.utils import Register
from rscder.plugins.basic import BasicPlugin
from PyQt5.QtWidgets import QWidget, QComboBox, QVBoxLayout
class MISCPlugin(BasicPlugin):
@staticmethod
def info():
return {
'name': 'MISC Plugin',
'description': '基本开发工具',
'author': 'Wangtong',
'version': '1.0.0',
}
def set_action(self):
pass
class AlgFrontend(object):
@staticmethod
def get_name():
return None
@staticmethod
def get_widget(parent=None):
return QWidget(parent)
@staticmethod
def get_params(widget=None):
return dict()
@staticmethod
def run_alg(*args, **kargs):
pass
class AlgSelectWidget(QWidget):
def __init__(self, parent= None, registery:Register = None) -> None:
super().__init__(parent)
self.reg = registery
if registery is None:
return
self.selectbox = QComboBox(self)
self.selectbox.addItem('--', None)
for key in self.reg.keys():
alg:AlgFrontend = self.reg[key]
if alg.get_name() is None:
name = key
else:
name = alg.get_name()
self.selectbox.addItem(name, key)
self.vbox = QVBoxLayout(self)
self.vbox.addWidget(self.selectbox)
self.params_widget = QWidget(self)
self.vbox.addWidget(self.params_widget)
self.selectbox.currentIndexChanged.connect(self.on_changed)
def on_changed(self, index):
if index == 0:
return
self.alg:AlgFrontend = self.reg[self.selectbox.itemData(index)]
self.vbox.removeWidget(self.params_widget)
self.params_widget = self.alg.get_widget(self)
self.vbox.addWidget(self.params_widget)
def get_alg_and_params(self):
if self.selectbox.currentIndex() == 0:
return None, None
return self.alg, self.alg.get_params(self.params_widget)

46
plugins/misc/table_layer.py Normal file
View File

@ -0,0 +1,46 @@
from datetime import datetime
import math
from random import random
from osgeo import gdal
from rscder.utils.project import Project
import os
from rscder.utils.project import BasicLayer, ResultPointLayer
def table_layer(pth:str,layer:BasicLayer, name, send_message = None):
if send_message is not None:
send_message.emit('正在计算表格结果...')
cell_size = layer.layer_parent.cell_size
ds = gdal.Open(pth)
xsize = ds.RasterXSize
ysize = ds.RasterYSize
geo = ds.GetGeoTransform()
out_csv = os.path.join(Project().other_path, f'{name}_table_result_{ int(datetime.now().timestamp() * 1000) }.csv')
yblocks = ysize // cell_size[1] + 1
xblocks = xsize // cell_size[0] + 1
with open(out_csv, 'w') as f:
f.write('x,y,diff,status\n')
for j in range(yblocks):
block_xy = (0, j * cell_size[1])
block_size = (xsize, cell_size[1])
if block_xy[1] + block_size[1] > ysize:
block_size = (xsize, ysize - block_xy[1])
block_data = ds.ReadAsArray(*block_xy, *block_size)
for i in range(xblocks):
start_x = i * cell_size[0]
end_x = start_x + cell_size[0]
if end_x > xsize:
end_x = xsize
block_data_xy = block_data[:, start_x:end_x]
center_x = start_x + cell_size[0] // 2
center_y = j * cell_size[1] + cell_size[1] // 2
center_x = center_x * geo[1] + geo [0]
center_y = center_y * geo[5] + geo [3]
f.write(f'{center_x},{center_y},{block_data_xy.mean() * 100},{int(block_data_xy.mean() > 0.5)}\n')
result_layer = ResultPointLayer(out_csv, enable=True, proj=layer.proj, geo=layer.geo,result_path={})
# print(result_layer.result_path)
layer.layer_parent.add_result_layer(result_layer)
if send_message is not None:
send_message.emit('计算完成')

42
plugins/misc/utils.py Normal file
View File

@ -0,0 +1,42 @@
import os
import random
import logging
# def generate_temp_file
class Register:
def __init__(self, registry_name):
self._dict = {}
self._name = registry_name
def __setitem__(self, key, value):
if not callable(value):
raise Exception(f"Value of a Registry must be a callable!\nValue: {value}")
if key is None:
key = value.__name__
if key in self._dict:
logging.warning("Key %s already in registry %s." % (key, self._name))
self._dict[key] = value
def register(self, target):
"""Decorator to register a function or class."""
def add(key, value):
self[key] = value
return value
if callable(target):
# @reg.register
return add(None, target)
# @reg.register('alias')
return lambda x: add(target, x)
def __getitem__(self, key):
return self._dict[key]
def __contains__(self, key):
return key in self._dict
def keys(self):
"""key"""
return self._dict.keys()

16
plugins/plugins.yaml
View File

@ -13,11 +13,11 @@
path: ./plugin\export_to
version: 1.0.0
- author: RSCDER
description: MeanFilter
description: FilterCollection
enabled: true
module: some_filter
name: MeanFilter
path: ./plugin\some_filter
module: filter_collection
name: FilterCollection
path: ./plugin\filter_collection
version: 1.0.0
- author: RSCDER
description: Evaluation
@ -27,11 +27,11 @@
path: ./plugin\evaluation
version: 1.0.0
- author: RSCDER
description: AllinOne
description: UnsupervisedPlugin
enabled: true
module: In_one
name: basic_diff_AllinOne
path: ./plugin\In_one
module: unsupervised_method
name: UnsupervisedPlugin
path: ./plugin\unsupervsied_method
version: 1.0.0
- author: RSCDER
description: set Change Rate

1
plugins/some_filter/__init__.py
View File

@ -1 +0,0 @@
from some_filter.main import *

5
plugins/thres/__init__.py Normal file
View File

@ -0,0 +1,5 @@
from misc import Register
THRES = Register('Thres')
from thres.main import *

178
plugins/thres/main.py Normal file
View File

@ -0,0 +1,178 @@
from osgeo import gdal, gdalconst
from rscder.utils.project import Project
import os
from misc import AlgFrontend
from thres import THRES
from PyQt5.QtWidgets import QWidget, QLabel, QHBoxLayout, QLineEdit
from PyQt5.QtGui import QIntValidator, QDoubleValidator
@THRES.register
class ManulGapAlg(AlgFrontend):
@staticmethod
def get_name():
return '手动阈值'
@staticmethod
def get_params(widget:QWidget=None):
if widget is None:
return dict(gap=125)
lineedit:QLineEdit = widget.layout().findChild(QLineEdit, 'lineedit')
if lineedit is None:
return dict(gap=125)
gap = int(lineedit.text())
return dict(gap=gap)
@staticmethod
def get_widget(parent=None):
widget = QWidget(parent)
lineedit = QLineEdit(widget)
lineedit.setObjectName('lineedit')
lineedit.setValidator(QIntValidator(1, 254))
# lineedit.
label = QLabel('阈值：')
layout = QHBoxLayout()
layout.addWidget(label)
layout.addWidget(lineedit)
widget.setLayout(layout)
return widget
@staticmethod
def run_alg(pth,name = '',gap = 125, send_message = None):
ds = gdal.Open(pth)
band = ds.GetRasterBand(1)
xsize = ds.RasterXSize
ysize = ds.RasterYSize
max_pixels = 1e7
max_rows = max_pixels // xsize
if max_rows < 1:
max_rows = 1
if max_rows > ysize:
max_rows = ysize
block_count = ysize // max_rows + 1
if send_message is not None:
send_message.emit('阈值为：{}'.format(gap))
out_th = os.path.join(Project().bcdm_path, '{}_thresh{}_bcdm.tif'.format(name,gap))
out_ds = gdal.GetDriverByName('GTiff').Create(out_th, xsize, ysize, 1, gdal.GDT_Byte)
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_ds.SetProjection(ds.GetProjection())
out_band = out_ds.GetRasterBand(1)
for i in range(block_count):
start_row = i * max_rows
end_row = min((i + 1) * max_rows, ysize)
block = band.ReadAsArray(0, start_row, xsize, end_row - start_row)
out_band.WriteArray(block > gap, 0, start_row)
out_band.FlushCache()
out_ds = None
ds = None
if send_message is not None:
send_message.emit('自定义阈值分割完成')
return out_th
import numpy as np
def OTSU(hist):
u1=0.0#背景像素的平均灰度值
u2=0.0#前景像素的平均灰度值
th=0.0
#总的像素数目
PixSum= np.sum(hist)
#各灰度值所占总像素数的比例
PixRate=hist / PixSum
#统计各个灰度值的像素个数
Max_var = 0
#确定最大类间方差对应的阈值
GrayScale = len(hist)
for i in range(1,len(hist)):#从1开始是为了避免w1为0.
u1_tem=0.0
u2_tem=0.0
#背景像素的比列
w1=np.sum(PixRate[:i])
#前景像素的比例
w2=1.0-w1
if w1==0 or w2==0:
pass
else:#背景像素的平均灰度值
for m in range(i):
u1_tem=u1_tem+PixRate[m]*m
u1 = u1_tem * 1.0 / w1
#前景像素的平均灰度值
for n in range(i,GrayScale):
u2_tem = u2_tem + PixRate[n]*n
u2 = u2_tem / w2
#print(u1)
#类间方差公式：G=w1*w2*(u1-u2)**2
tem_var=w1*w2*np.power((u1-u2),2)
#print(tem_var)
#判断当前类间方差是否为最大值。
if Max_var<tem_var:
Max_var=tem_var#深拷贝，Max_var与tem_var占用不同的内存空间。
th=i
return th
@THRES.register
class OTSUAlg(AlgFrontend):
@staticmethod
def get_name():
return 'OTSU阈值'
@staticmethod
def run_alg(pth,name='',send_message=None):
ds = gdal.Open(pth)
band = ds.GetRasterBand(1)
# band_count = ds.RasterCount
hist = np.zeros(256, dtype=np.int)
xsize = ds.RasterXSize
ysize = ds.RasterYSize
max_pixels = 1e7
max_rows = max_pixels // xsize
if max_rows < 1:
max_rows = 1
if max_rows > ysize:
max_rows = ysize
block_count = ysize // max_rows + 1
for i in range(block_count):
start_row = i * max_rows
end_row = min((i + 1) * max_rows, ysize)
block = band.ReadAsArray(0, start_row, xsize, end_row - start_row)
hist += np.histogram(block.flatten(), bins=256, range=(0, 255))[0]
hist = hist.astype(np.float32)
gap = OTSU(hist)
send_message.emit('阈值为：{}'.format(gap))
out_th = os.path.join(Project().bcdm_path, '{}_otsu_bcdm.tif'.format(name))
out_ds = gdal.GetDriverByName('GTiff').Create(out_th, xsize, ysize, 1, gdal.GDT_Byte)
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_ds.SetProjection(ds.GetProjection())
out_band = out_ds.GetRasterBand(1)
for i in range(block_count):
start_row = i * max_rows
end_row = min((i + 1) * max_rows, ysize)
block = band.ReadAsArray(0, start_row, xsize, end_row - start_row)
out_band.WriteArray(block > gap, 0, start_row)
out_band.FlushCache()
out_ds = None
ds = None
send_message.emit('OTSU阈值完成')
return out_th
# ManulGapAlg.run_alg

1
plugins/unsupervised_method/__init__.py Normal file
View File

@ -0,0 +1 @@
from unsupervised_method.main import *

202
plugins/unsupervised_method/main.py Normal file
View File

@ -0,0 +1,202 @@
from threading import Thread
from plugins.misc.main import AlgFrontend
from rscder.gui.actions import ActionManager
from rscder.plugins.basic import BasicPlugin
from PyQt5.QtWidgets import QAction, QToolBar, QMenu, QDialog, QHBoxLayout, QVBoxLayout, QPushButton,QWidget,QLabel,QLineEdit,QPushButton,QComboBox,QDialogButtonBox
from rscder.gui.layercombox import PairLayerCombox
from rscder.utils.icons import IconInstance
from filter_collection import FILTER
from .scripts import UNSUPER_CD
from thres import THRES
from misc import table_layer, AlgSelectWidget
class UnsupervisedCDMethod(QDialog):
def __init__(self,parent=None, alg:AlgFrontend=None):
super(UnsupervisedCDMethod, self).__init__(parent)
self.alg = alg
self.setWindowTitle('无监督变化检测:{}'.format(alg.get_name()))
self.setWindowIcon(IconInstance().LOGO)
self.initUI()
self.setMinimumWidth(500)
def initUI(self):
#图层
self.layer_combox = PairLayerCombox(self)
layerbox = QHBoxLayout()
layerbox.addWidget(self.layer_combox)
self.filter_select = AlgSelectWidget(self, FILTER)
self.param_widget = self.alg.get_widget(self)
self.unsupervised_menu = self.param_widget
self.thres_select = AlgSelectWidget(self, THRES)
self.ok_button = QPushButton('确定', self)
self.ok_button.setIcon(IconInstance().OK)
self.ok_button.clicked.connect(self.accept)
self.ok_button.setDefault(True)
self.cancel_button = QPushButton('取消', self)
self.cancel_button.setIcon(IconInstance().CANCEL)
self.cancel_button.clicked.connect(self.reject)
self.cancel_button.setDefault(False)
buttonbox=QDialogButtonBox(self)
buttonbox.addButton(self.ok_button,QDialogButtonBox.NoRole)
buttonbox.addButton(self.cancel_button,QDialogButtonBox.NoRole)
buttonbox.setCenterButtons(True)
totalvlayout=QVBoxLayout()
totalvlayout.addLayout(layerbox)
totalvlayout.addWidget(self.filter_select)
if self.param_widget is not None:
totalvlayout.addWidget(self.param_widget)
totalvlayout.addWidget(self.thres_select)
totalvlayout.addWidget(buttonbox)
totalvlayout.addStretch()
self.setLayout(totalvlayout)
class UnsupervisedCD(QDialog):
def __init__(self,parent=None):
super(UnsupervisedCD, self).__init__(parent)
self.setWindowTitle('无监督变化检测')
self.setWindowIcon(IconInstance().LOGO)
self.setMinimumWidth(600)
self.initUI()
def initUI(self):
#图层
self.layer_combox = PairLayerCombox(self)
layerbox = QHBoxLayout()
layerbox.addWidget(self.layer_combox)
self.filter_select = AlgSelectWidget(self, FILTER)
self.unsupervised_select = AlgSelectWidget(self, UNSUPER_CD)
self.thres_select = AlgSelectWidget(self, THRES)
self.ok_button = QPushButton('确定', self)
self.ok_button.setIcon(IconInstance().OK)
self.ok_button.clicked.connect(self.accept)
self.ok_button.setDefault(True)
self.cancel_button = QPushButton('取消', self)
self.cancel_button.setIcon(IconInstance().CANCEL)
self.cancel_button.clicked.connect(self.reject)
self.cancel_button.setDefault(False)
buttonbox=QDialogButtonBox(self)
buttonbox.addButton(self.ok_button,QDialogButtonBox.NoRole)
buttonbox.addButton(self.cancel_button,QDialogButtonBox.NoRole)
buttonbox.setCenterButtons(True)
totalvlayout=QVBoxLayout()
totalvlayout.addLayout(layerbox)
totalvlayout.addWidget(self.filter_select)
totalvlayout.addWidget(self.unsupervised_select)
totalvlayout.addWidget(self.thres_select)
totalvlayout.addWidget(buttonbox)
totalvlayout.addStretch()
self.setLayout(totalvlayout)
class UnsupervisedPlugin(BasicPlugin):
@staticmethod
def info():
return {
'name': 'UnsupervisedPlugin',
'description': 'UnsupervisedPlugin',
'author': 'RSCDER',
'version': '1.0.0',
}
def set_action(self):
unsupervised_menu = QMenu('&无监督变化检测', self.mainwindow)
unsupervised_menu.setIcon(IconInstance().UNSUPERVISED)
ActionManager().change_detection_menu.addMenu(unsupervised_menu)
basic_diff_method_in_one = QAction(IconInstance().UNSUPERVISED, '&无监督变化检测', self.mainwindow)
basic_diff_method_in_one.triggered.connect(self.run)
toolbar:QToolBar = ActionManager().add_toolbar('无监督变化检测')
toolbar.addAction(basic_diff_method_in_one)
# self.mainwindow.addToolbar(toolbar)
# print(UNSUPER_CD.keys())
for key in UNSUPER_CD.keys():
alg:AlgFrontend = UNSUPER_CD[key]
if alg.get_name() is None:
name = key
else:
name = alg.get_name()
action = QAction(name, unsupervised_menu)
action.triggered.connect(lambda : self.run_cd(alg))
unsupervised_menu.addAction(action)
def run(self):
myDialog= UnsupervisedCD(self.mainwindow)
myDialog.show()
if myDialog.exec_()==QDialog.Accepted:
w=myDialog
t=Thread(target=self.run_alg,args=(w,))
t.start()
def run_cd(self, alg):
dialog = UnsupervisedCDMethod(self.mainwindow, alg)
dialog.show()
if dialog.exec_() == QDialog.Accepted:
t = Thread(target=self.run_cd_alg, args=(dialog,))
t.start()
def run_cd_alg(self, w:UnsupervisedCDMethod):
layer1=w.layer_combox.layer1
pth1 = w.layer_combox.layer1.path
pth2 = w.layer_combox.layer2.path
name = layer1.layer_parent.name
falg, fparams = w.filter_select.get_alg_and_params()
cdalg = w.alg
cdparams = w.alg.get_params()
thalg, thparams = w.thres_select.get_alg_and_params()
if cdalg is None or thalg is None:
return
if falg is not None:
pth1 = falg.run_alg(pth1, name=name, send_message=self.send_message, **fparams)
pth2 = falg.run_alg(pth2, name=name, send_message=self.send_message, **fparams)
cdpth = cdalg.run_alg(pth1, pth2, layer1.layer_parent, send_message=self.send_message,**cdparams)
thpth = thalg.run_alg(cdpth, name=name, send_message=self.send_message, **thparams)
table_layer(thpth,layer1,name,self.send_message)
def run_alg(self,w:UnsupervisedCD):
layer1=w.layer_combox.layer1
pth1 = w.layer_combox.layer1.path
pth2 = w.layer_combox.layer2.path
name = layer1.layer_parent.name
falg, fparams = w.filter_select.get_alg_and_params()
cdalg, cdparams = w.unsupervised_select.get_alg_and_params()
thalg, thparams = w.thres_select.get_alg_and_params()
if cdalg is None or thalg is None:
return
if falg is not None:
pth1 = falg.run_alg(pth1, name=name, send_message=self.send_message, **fparams)
pth2 = falg.run_alg(pth2, name=name, send_message=self.send_message, **fparams)
cdpth = cdalg.run_alg(pth1, pth2, layer1.layer_parent, send_message=self.send_message,**cdparams)
thpth = thalg.run_alg(cdpth, name=name, send_message=self.send_message, **thparams)
table_layer(thpth,layer1,name,self.send_message)

0
plugins/In_one/pic.py → plugins/unsupervised_method/pic.py
View File

0
plugins/unsupervised_method/readme.md Normal file
View File

0
plugins/In_one/scripts/ACD.py → plugins/unsupervised_method/scripts/ACD.py
View File

0
plugins/In_one/scripts/AHT.py → plugins/unsupervised_method/scripts/AHT.py
View File

0
plugins/In_one/scripts/LHBA.py → plugins/unsupervised_method/scripts/LHBA.py
View File

0
plugins/In_one/scripts/OCD.py → plugins/unsupervised_method/scripts/OCD.py
View File

0
plugins/In_one/scripts/SH.py → plugins/unsupervised_method/scripts/SH.py
View File

688
plugins/unsupervised_method/scripts/__init__.py Normal file
View File

@ -0,0 +1,688 @@
from datetime import datetime
from osgeo import gdal
import math,os
import time
from PyQt5 import QtWidgets
from sklearn.cluster import k_means
from rscder.utils.geomath import geo2imageRC, imageRC2geo
from rscder.utils.project import Project, PairLayer
from misc import Register, AlgFrontend
UNSUPER_CD = Register('unsuper cd')
import numpy as np
from .ACD import ACD
from .AHT import AHT
from .OCD import OCD
from .LHBA import LHBA
from .SH import SH
def warp(file,ds:gdal.Dataset,srcWin=[0,0,0,0]):
driver = gdal.GetDriverByName('GTiff')
xsize=ds.RasterXSize
ysize=ds.RasterYSize
geo=ds.GetGeoTransform()
orj=ds.GetProjection()
band=ds.RasterCount
if os.path.exists(file):
os.remove(file)
out_ds:gdal.Dataset=driver.Create(file, xsize, ysize, band, gdal.GDT_Byte)
out_ds.SetGeoTransform(geo)
out_ds.SetProjection(orj)
for b in range(1,band+1):
out_ds.GetRasterBand(b).WriteArray(ds.ReadAsArray(*srcWin,band_list=[b]),*(0,0))
del out_ds
@UNSUPER_CD.register
class BasicCD(AlgFrontend):
@staticmethod
def get_name():
return '差分法'
@staticmethod
def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None,*args, **kargs):
ds1:gdal.Dataset=gdal.Open(pth1)
ds2:gdal.Dataset=gdal.Open(pth2)
cell_size = layer_parent.cell_size
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
band = ds1.RasterCount
yblocks = ysize // cell_size[1]
driver = gdal.GetDriverByName('GTiff')
out_tif = os.path.join(Project().other_path, 'temp.tif')
out_ds = driver.Create(out_tif, xsize, ysize, 1, gdal.GDT_Float32)
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
out_ds.SetGeoTransform(geo)
out_ds.SetProjection(proj)
max_diff = 0
min_diff = math.inf
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
for j in range(yblocks + 1):#该改这里了
if send_message is not None:
send_message.emit(f'计算{j}/{yblocks}')
block_xy1 = (start1x, start1y+j * cell_size[1])
block_xy2 = (start2x,start2y+j*cell_size[1])
block_xy=(0,j * cell_size[1])
if block_xy1[1] > end1y or block_xy2[1] > end2y:
break
block_size=(xsize, cell_size[1])
block_size1 = (xsize, cell_size[1])
block_size2 = (xsize,cell_size[1])
if block_xy[1] + block_size[1] > ysize:
block_size = (xsize, ysize - block_xy[1])
if block_xy1[1] + block_size1[1] > end1y:
block_size1 = (xsize,end1y - block_xy1[1])
if block_xy2[1] + block_size2[1] > end2y:
block_size2 = (xsize, end2y - block_xy2[1])
block_data1 = ds1.ReadAsArray(*block_xy1, *block_size1)
block_data2 = ds2.ReadAsArray(*block_xy2, *block_size2)
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)
if send_message is not None:
send_message.emit(f'完成{j}/{yblocks}')
del ds2
del ds1
out_ds.FlushCache()
del out_ds
if send_message is not None:
send_message.emit('归一化概率中...')
temp_in_ds = gdal.Open(out_tif)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(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(geo)
out_normal_ds.SetProjection(proj)
# 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
del out_normal_ds
try:
os.remove(out_tif)
except:
pass
if send_message is not None:
send_message.emit('差分法计算完成')
return out_normal_tif
@UNSUPER_CD.register
class LSTS(AlgFrontend):
@staticmethod
def get_name():
return 'LSTS'
@staticmethod
def get_widget(parent=None):
widget = QtWidgets.QWidget(parent)
def get_params(self):
return dict(n=5, w_size=(3,3))
@staticmethod
def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message=None,n=5,w_size=(3,3), *args, **kws):
ds1:gdal.Dataset=gdal.Open(pth1)
ds2:gdal.Dataset=gdal.Open(pth2)
cell_size = layer_parent.cell_size
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
band = ds1.RasterCount
yblocks = ysize // cell_size[1]
driver = gdal.GetDriverByName('GTiff')
out_tif = os.path.join(Project().other_path, '%d.tif'%(int(datetime.now().timestamp() * 1000)))
out_ds = driver.Create(out_tif, xsize, ysize, 1, gdal.GDT_Float32)
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
out_ds.SetGeoTransform(geo)
out_ds.SetProjection(proj)
pixnum=w_size[0]*w_size[1]
# send_message.emit('pixnum:'pixnum)
max_diff = 0
min_diff = math.inf
win_h=w_size[0]//2 #half hight of window
win_w=w_size[1]//2 #half width of window
a=[[(i+1)**j for j in range(n+1)] for i in range(pixnum)]
A=np.array(a).astype(np.float64)#
k_=np.array(range(1,n+1))
df1=np.zeros(pixnum).astype(np.float64)
df2=np.zeros(pixnum).astype(np.float64)
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
for j in range(yblocks + 1):
if send_message is not None:
send_message.emit(f'计算{j}/{yblocks}')
block_xy1 = (start1x, start1y+j * cell_size[1])
block_xy2 = (start2x,start2y+j*cell_size[1])
block_xy=(0,j * cell_size[1])
if block_xy1[1] > end1y or block_xy2[1] > end2y:
break
block_size=(xsize, cell_size[1])
block_size1 = (xsize, cell_size[1])
block_size2 = (xsize,cell_size[1])
if block_xy[1] + block_size[1] > ysize:
block_size = (xsize, ysize - block_xy[1])
if block_xy1[1] + block_size1[1] > end1y:
block_size1 = (xsize,end1y - block_xy1[1])
if block_xy2[1] + block_size2[1] > end2y:
block_size2 = (xsize, end2y - block_xy2[1])
block_data1 = ds1.ReadAsArray(*block_xy1, *block_size1)
block_data2 = ds2.ReadAsArray(*block_xy2, *block_size2)
if band == 1:
block_data1 = block_data1[None, ...]
block_data2 = block_data2[None, ...]
# pdb.set_trace()
else:
block_data1=np.mean(block_data1,0)
block_data2=np.mean(block_data2,0)
block_diff=np.zeros(block_data1.shape).astype(np.float64)
for i in range(win_h,block_size1[1]-win_h):
for j_ in range(win_w,block_size1[0]-win_w):
pix=0
#get b
# b1=block_data[i+win_h:i+win_h] c in range(j_-win_w,j_+win_w+1)
b1=block_data1[i-win_h:i+win_h+1,j_-win_w:j_+win_w+1]
b2=block_data2[i-win_h:i+win_h+1,j_-win_w:j_+win_w+1]
b1=[b if (r+1)//2 else b[::-1] for r,b in enumerate(b1)]
b2=[b if (r+1)//2 else b[::-1] for r,b in enumerate(b2)]
b1=np.expand_dims(np.concatenate(b1,0),1)
b2=np.expand_dims(np.concatenate(b2,0),1)
x1=np.squeeze(np.linalg.pinv(A).dot(b1))
x2=np.squeeze(np.linalg.pinv(A).dot(b2))
#df
k_=range(1,n+1)
for pix in range(1,pixnum+1):
df1[pix-1]=x1[1:n+1].dot(np.array([k*(pix**(k-1)) for k in k_]))
df2[pix-1]=x2[1:n+1].dot(np.array([k*(pix**(k-1)) for k in k_]))
#distance 欧式距离
block_diff[i][j_]=np.dot(df1-df2,df1-df2)**0.5
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)
send_message.emit(f'完成{j}/{yblocks}')
del ds2
del ds1
out_ds.FlushCache()
del out_ds
if send_message is not None:
send_message.emit('归一化概率中...')
temp_in_ds = gdal.Open(out_tif)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(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(geo)
out_normal_ds.SetProjection(proj)
# 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
del out_normal_ds
try:
os.remove(out_tif)
except:
pass
if send_message is not None:
send_message.emit('LSTS法计算完成')
return out_normal_tif
@UNSUPER_CD.register
class CVAAlg(AlgFrontend):
@staticmethod
def get_name():
return 'CVA'
@staticmethod
def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None, *args, **kargs):
ds1:gdal.Dataset=gdal.Open(pth1)
ds2:gdal.Dataset=gdal.Open(pth2)
cell_size = layer_parent.cell_size
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
band = ds1.RasterCount
yblocks = ysize // cell_size[1]
driver = gdal.GetDriverByName('GTiff')
out_tif = os.path.join(Project().other_path, 'temp.tif')
out_ds = driver.Create(out_tif, xsize, ysize, 1, gdal.GDT_Float32)
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
out_ds.SetGeoTransform(geo)
out_ds.SetProjection(proj)
max_diff = 0
min_diff = math.inf
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
for j in range(yblocks + 1):
if send_message is not None:
send_message.emit(f'计算{j}/{yblocks}')
block_xy1 = (start1x, start1y+j * cell_size[1])
block_xy2 = (start2x,start2y+j*cell_size[1])
block_xy=(0,j * cell_size[1])
if block_xy1[1] > end1y or block_xy2[1] > end2y:
break
block_size=(xsize, cell_size[1])
block_size1 = (xsize, cell_size[1])
block_size2 = (xsize,cell_size[1])
if block_xy[1] + block_size[1] > ysize:
block_size = (xsize, ysize - block_xy[1])
if block_xy1[1] + block_size1[1] > end1y:
block_size1 = (xsize,end1y - block_xy1[1])
if block_xy2[1] + block_size2[1] > end2y:
block_size2 = (xsize, end2y - block_xy2[1])
block_data1 = ds1.ReadAsArray(*block_xy1, *block_size1)
block_data2 = ds2.ReadAsArray(*block_xy2, *block_size2)
if band == 1:
block_data1 = block_data1[None, ...]
block_data2 = block_data2[None, ...]
# pdb.set_trace()
block_diff=np.sum((block_data1-block_data2)**2,0)**0.5
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)
if send_message is not None:
send_message.emit(f'完成{j}/{yblocks}')
del ds2
del ds1
out_ds.FlushCache()
del out_ds
if send_message is not None:
send_message.emit('归一化概率中...')
temp_in_ds = gdal.Open(out_tif)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(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(geo)
out_normal_ds.SetProjection(proj)
# 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
del out_normal_ds
try:
os.remove(out_tif)
except:
pass
if send_message is not None:
send_message.emit('欧式距离计算完成')
return out_normal_tif
@UNSUPER_CD.register
class ACDAlg(AlgFrontend):
@staticmethod
def get_name():
return 'ACD'
@staticmethod
def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None, *args, **kargs):
if send_message is None:
class Empty:
def emit(self, *args, **kws):
print(args)
send_message = Empty()
# send_message.emit = print
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
#提取公共部分
send_message.emit('提取重叠区域数据.....')
ds2:gdal.Dataset=gdal.Open(pth2)
temp_tif2 = os.path.join(Project().other_path,'temp2.tif')
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize])
del ds2
send_message.emit('图像二提取完成')
ds1:gdal.Dataset=gdal.Open(pth1)
temp_tif1 = os.path.join(Project().other_path, 'temp1.tif')
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize])
del ds1
send_message.emit('图像一提取完成')
#运算
send_message.emit('开始ACD计算.....')
time.sleep(0.1)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000)))
ACD(temp_tif1,temp_tif2,out_normal_tif)
#添加投影
send_message.emit('录入投影信息.....')
time.sleep(0.1)
ds=gdal.Open(out_normal_tif,1)
ds.SetGeoTransform(geo)
ds.SetProjection(proj)
del ds
return out_normal_tif
@UNSUPER_CD.register
class AHTAlg(AlgFrontend):
@staticmethod
def get_name():
return 'AHT'
@staticmethod
def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None, *args, **kargs):
if send_message is None:
class Empty:
def emit(self, *args, **kws):
print(args)
send_message = Empty()
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
#提取公共部分
send_message.emit('提取重叠区域数据.....')
ds2:gdal.Dataset=gdal.Open(pth2)
temp_tif2 = os.path.join(Project().other_path,'temp2.tif')
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize])
del ds2
send_message.emit('图像二提取完成')
ds1:gdal.Dataset=gdal.Open(pth1)
temp_tif1 = os.path.join(Project().other_path, 'temp1.tif')
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize])
del ds1
send_message.emit('图像一提取完成')
#运算
send_message.emit('开始AHT计算.....')
time.sleep(0.1)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000)))
AHT(temp_tif1,temp_tif2,out_normal_tif)
#添加投影
send_message.emit('录入投影信息.....')
time.sleep(0.1)
ds=gdal.Open(out_normal_tif,1)
ds.SetGeoTransform(geo)
ds.SetProjection(proj)
del ds
return out_normal_tif
@UNSUPER_CD.register
class OCDAlg(AlgFrontend):
@staticmethod
def get_name():
return 'OCD'
@staticmethod
def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None, *args, **kargs):
if send_message is None:
class Empty:
def emit(self, *args, **kws):
print(args)
send_message = Empty()
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
#提取公共部分
send_message.emit('提取重叠区域数据.....')
ds2:gdal.Dataset=gdal.Open(pth2)
temp_tif2 = os.path.join(Project().other_path,'temp2.tif')
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize])
del ds2
send_message.emit('图像二提取完成')
ds1:gdal.Dataset=gdal.Open(pth1)
temp_tif1 = os.path.join(Project().other_path, 'temp1.tif')
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize])
del ds1
send_message.emit('图像一提取完成')
#运算
send_message.emit('开始OCD计算.....')
time.sleep(0.1)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000)))
OCD(temp_tif1,temp_tif2,out_normal_tif,Project().other_path)
#添加投影
send_message.emit('录入投影信息.....')
time.sleep(0.1)
ds=gdal.Open(out_normal_tif,1)
ds.SetGeoTransform(geo)
ds.SetProjection(proj)
del ds
return out_normal_tif
@UNSUPER_CD.register
class LHBAAlg(AlgFrontend):
@staticmethod
def get_name():
return 'LHBA'
@staticmethod
def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None, *args, **kargs):
if send_message is None:
class Empty:
def emit(self, *args, **kws):
print(args)
send_message = Empty()
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
#提取公共部分
send_message.emit('提取重叠区域数据.....')
ds2:gdal.Dataset=gdal.Open(pth2)
temp_tif2 = os.path.join(Project().other_path,'temp2.tif')
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize])
del ds2
send_message.emit('图像二提取完成')
ds1:gdal.Dataset=gdal.Open(pth1)
temp_tif1 = os.path.join(Project().other_path, 'temp1.tif')
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize])
del ds1
send_message.emit('图像一提取完成')
#运算
send_message.emit('开始LHBA计算.....')
time.sleep(0.1)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000)))
LHBA(temp_tif1,temp_tif2,out_normal_tif)
#添加投影
send_message.emit('录入投影信息.....')
time.sleep(0.1)
ds=gdal.Open(out_normal_tif,1)
ds.SetGeoTransform(geo)
ds.SetProjection(proj)
del ds
return out_normal_tif
@UNSUPER_CD.register
class SHAlg(AlgFrontend):
@staticmethod
def get_name():
return 'SH'
@staticmethod
def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None, *args, **kargs):
if send_message is None:
class Empty:
def emit(self, *args, **kws):
print(args)
send_message = Empty()
xsize = layer_parent.size[0]
ysize = layer_parent.size[1]
geo=layer_parent.grid.geo
proj=layer_parent.grid.proj
#提取公共部分
send_message.emit('提取重叠区域数据.....')
ds2:gdal.Dataset=gdal.Open(pth2)
temp_tif2 = os.path.join(Project().other_path,'temp2.tif')
start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize])
del ds2
send_message.emit('图像二提取完成')
ds1:gdal.Dataset=gdal.Open(pth1)
temp_tif1 = os.path.join(Project().other_path, 'temp1.tif')
start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1])
end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3])
warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize])
del ds1
send_message.emit('图像一提取完成')
#运算
send_message.emit('开始SH计算.....')
time.sleep(0.1)
out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000)))
SH(temp_tif1,temp_tif2,out_normal_tif)
#添加投影
send_message.emit('录入投影信息.....')
time.sleep(0.1)
ds=gdal.Open(out_normal_tif,1)
ds.SetGeoTransform(geo)
ds.SetProjection(proj)
del ds
return out_normal_tif

27
rscder/gui/actions.py
View File

@ -2,7 +2,8 @@ import logging
import os
from pathlib import Path
from PyQt5 import QtCore, QtGui, QtWidgets
from PyQt5.QtWidgets import QAction, QActionGroup, QLabel, QFileDialog, QMenuBar
from PyQt5.QtCore import Qt, QSize, QSettings, pyqtSignal
from PyQt5.QtWidgets import QAction, QActionGroup, QLabel, QFileDialog, QMenuBar, QToolBar
from rscder.gui import project
from rscder.gui.project import Create
from rscder.utils.icons import IconInstance
@ -32,6 +33,8 @@ class ActionManager(QtCore.QObject):
self.action_groups = {}
self.action_group_actions = {}
self.toolbars = {}
self.double_map = double_map
self.layer_tree = layer_tree
self.follow_box = follow_box
@ -74,9 +77,21 @@ class ActionManager(QtCore.QObject):
self.help_menu = menubar.addMenu( '&帮助')
def set_toolbar(self, toolbar):
self.toolbar = toolbar
self.toolbar:QToolBar = toolbar
self.toolbar.setIconSize(QtCore.QSize(24, 24))
def add_toolbar(self, name=None):
toolbar = self.w_parent.addToolBar(name)
toolbar.setMovable(True)
toolbar.setFloatable(False)
toolbar.setIconSize(QSize(32, 32))
toolbar.setToolButtonStyle(Qt.ToolButtonTextBesideIcon)
toolbar.setContextMenuPolicy(Qt.PreventContextMenu)
toolbar.setLayoutDirection(Qt.LeftToRight)
return toolbar
def set_status_bar(self, status_bar):
self.status_bar = status_bar
@ -131,6 +146,9 @@ class ActionManager(QtCore.QObject):
zomm_out.setChecked(False)
zomm_in.setCheckable(True)
zomm_in.setChecked(False)
toolbar = self.add_toolbar('Map')
toolbar.addActions([pan, zomm_out, zomm_in, locate])
self.double_map.connect_map_tool(pan, zomm_in, zomm_out)
# self.double_map.connect_grid_show(grid_line)
@ -148,6 +166,9 @@ class ActionManager(QtCore.QObject):
plugin_list = self.add_action(QAction(IconInstance().PLUGINS,'&插件列表', self.w_parent), 'Plugin')
plugin_list.triggered.connect(self.plugin_list)
# toolbar = self.add_toolbar('Plugin')
# toolbar.addAction(plugin_list)
self.plugin_menu.addAction(plugin_list)
self.message_box.info('菜单初始化完成')

3
rscder/plugins/loader.py
View File

@ -62,7 +62,10 @@ class PluginLoader(QObject):
self.plugins.append(obj(self.ctx))
break
except Exception as e:
# import traceback
# traceback.print_exc()
self.ctx['message_box'].error(f'{plugin["name"]} load error: {e}')
self.plugin_loaded.emit()