diff --git a/ReadMe.md b/ReadMe.md index e3c54f5..e438c3d 100644 --- a/ReadMe.md +++ b/ReadMe.md @@ -6,7 +6,15 @@ CVEO团队 Python3.7 + PyQt + QGIS ## 配置方式 + +``` conda create -f conda.yaml +``` +或 + +``` +conda env create -n cevo_qt -f conda.yaml +``` # 打包方式 diff --git a/conda.yaml b/conda.yaml index c68b502..8310339 100644 --- a/conda.yaml +++ b/conda.yaml @@ -11,6 +11,7 @@ dependencies: - pyqtads=3.8.2=py37hf2a7229_0 - python=3.7.10=h7840368_101_cpython - python_abi=3.7=2_cp37m + - pip - pyyaml=5.4.1=py37hcc03f2d_1 - qgis=3.18.3=py37h3dc7164_2 - qt=5.12.9=h5909a2a_4 diff --git a/icons/Algorithm_icon/功能-01.png b/icons/Algorithm_icon/功能-01.png new file mode 100644 index 0000000..1b2ad58 Binary files /dev/null and b/icons/Algorithm_icon/功能-01.png differ diff --git a/icons/Algorithm_icon/功能-02.png b/icons/Algorithm_icon/功能-02.png new file mode 100644 index 0000000..d629be0 Binary files /dev/null and b/icons/Algorithm_icon/功能-02.png differ diff --git a/icons/Algorithm_icon/功能-03.png b/icons/Algorithm_icon/功能-03.png new file mode 100644 index 0000000..dccadc3 Binary files /dev/null and b/icons/Algorithm_icon/功能-03.png differ diff --git a/icons/植被变化.png b/icons/植被变化.png new file mode 100644 index 0000000..a565fc8 Binary files /dev/null and b/icons/植被变化.png differ diff --git a/lic/license.lic b/lic/license.lic index fcdb940..0b622e9 100644 --- a/lic/license.lic +++ b/lic/license.lic @@ -1 +1 @@ -vd4FiYncytyziGH9GNCAA8hGGr1/79Xmphtc5+PHPJDpxvqj1hP7+985QMojYO4M5Qn/aqEAvFgeDN3CA8x1YAK8SdCgSXSBJpRBK8wqPQjBY1ak96QfdPCrTLunr+xuPxK3Gxe772adTTsee2+ot7WePYUsC4y4NcS5+rlP1if87xtYqVeSwx3c64cOmAGP \ No newline at end of file +IEhM2L0c7TWXeIEGUm8n1nMToPxkWiapMjEFtLFLHad4a1+2IKcrD2wbRobBLmTIdux5iOV+AV+t2JY75+w5NAU+umtYPWmpbAiphAQ1S8s8/h7fclRBU6Ym8oXDxCmZoWpCLCpJd8lE/R3nfKBoexDRYXecz/lw58wkdZbHktpNVKfjoe827/aKa7gmYNXO \ No newline at end of file diff --git a/plugins/unsupervised_method/scripts/__init__.py b/plugins/unsupervised_method/scripts/__init__.py index f0fe433..c432083 100644 --- a/plugins/unsupervised_method/scripts/__init__.py +++ b/plugins/unsupervised_method/scripts/__init__.py @@ -1,6 +1,14 @@ +from rscder.utils.icons import IconInstance +from .SH import SH +from .LHBA import LHBA +from .OCD import OCD +from .AHT import AHT +from .ACD import ACD +import numpy as np from datetime import datetime from osgeo import gdal -import math,os +import math +import os import time from PyQt5 import QtWidgets from sklearn.cluster import k_means @@ -10,42 +18,43 @@ from misc import Register, AlgFrontend UNSUPER_CD = Register('无监督变化检测方法') -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]): +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 + 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: 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)) + 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) - + def get_icon(): + return IconInstance().ARITHMETIC3 + + @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] @@ -56,48 +65,52 @@ class BasicCD(AlgFrontend): 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 + 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]) + 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]) - for j in range(yblocks + 1):#该改这里了 + 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]) + 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]) + 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]) + block_size1 = (xsize, end1y - block_xy1[1]) if block_xy2[1] + block_size2[1] > end2y: - block_size2 = (xsize, end2y - block_xy2[1]) + 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, ...] + 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) @@ -110,10 +123,12 @@ class BasicCD(AlgFrontend): 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) + 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) @@ -141,6 +156,7 @@ class BasicCD(AlgFrontend): send_message.emit('差分法计算完成') return out_normal_tif + @UNSUPER_CD.register class LSTS(AlgFrontend): @@ -148,6 +164,10 @@ class LSTS(AlgFrontend): def get_name(): return 'LSTS' + @staticmethod + def get_icon(): + return IconInstance().ARITHMETIC3 + @staticmethod def get_widget(parent=None): @@ -157,13 +177,13 @@ class LSTS(AlgFrontend): @staticmethod def get_params(widget=None): - return dict(n=5, w_size=(3,3)) - + 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) - + 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] @@ -172,84 +192,91 @@ class LSTS(AlgFrontend): 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_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 + 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] + 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) + 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) - 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]) + k_ = np.array(range(1, n+1)) + df1 = np.zeros(pixnum).astype(np.float64) + df2 = np.zeros(pixnum).astype(np.float64) - 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]) + 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]) + 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]) + 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]) + block_size1 = (xsize, end1y - block_xy1[1]) if block_xy2[1] + block_size2[1] > end2y: - block_size2 = (xsize, end2y - block_xy2[1]) + 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, ...] + 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) + 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 + 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) + 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 - 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) @@ -261,10 +288,12 @@ class LSTS(AlgFrontend): 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) + 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) @@ -288,7 +317,7 @@ class LSTS(AlgFrontend): os.remove(out_tif) except: pass - + if send_message is not None: send_message.emit('LSTS法计算完成') return out_normal_tif @@ -302,11 +331,15 @@ class CVAAlg(AlgFrontend): return 'CVA' @staticmethod - def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None, *args, **kargs): + def get_icon(): + return IconInstance().ARITHMETIC3 + + @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) - 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] @@ -317,45 +350,48 @@ class CVAAlg(AlgFrontend): 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 + 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]) + 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]) + 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]) + 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]) + 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]) + block_size1 = (xsize, end1y - block_xy1[1]) if block_xy2[1] + block_size2[1] > end2y: - block_size2 = (xsize, end2y - block_xy2[1]) + 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, ...] + 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 + 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) @@ -367,10 +403,12 @@ class CVAAlg(AlgFrontend): 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) + 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) @@ -398,19 +436,24 @@ class CVAAlg(AlgFrontend): 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): - + def get_icon(): + return IconInstance().ARITHMETIC3 + + @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() @@ -418,43 +461,45 @@ class ACDAlg(AlgFrontend): xsize = layer_parent.size[0] ysize = layer_parent.size[1] - geo=layer_parent.grid.geo - proj=layer_parent.grid.proj - #提取公共部分 + 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]) + 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) + 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]) + 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) - #添加投影 + 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 = gdal.Open(out_normal_tif, 1) ds.SetGeoTransform(geo) ds.SetProjection(proj) del ds - + return out_normal_tif @@ -464,56 +509,62 @@ 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): - + def get_icon(): + return IconInstance().ARITHMETIC3 + + @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 - #提取公共部分 + 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]) + 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) + 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]) + 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) - #添加投影 + 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 = gdal.Open(out_normal_tif, 1) ds.SetGeoTransform(geo) ds.SetProjection(proj) del ds - + return out_normal_tif @@ -523,109 +574,123 @@ 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): - + def get_icon(): + return IconInstance().ARITHMETIC3 + + @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 - #提取公共部分 + 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]) + 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) + 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]) + 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) - #添加投影 + 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 = 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): - + def get_icon(): + return IconInstance().ARITHMETIC3 + + @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 - #提取公共部分 + 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]) + 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) + 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]) + 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) - #添加投影 + 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 = gdal.Open(out_normal_tif, 1) ds.SetGeoTransform(geo) ds.SetProjection(proj) del ds @@ -638,52 +703,58 @@ 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): - + def get_icon(): + return IconInstance().ARITHMETIC3 + + @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 - #提取公共部分 + 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]) + 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) + 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]) + 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) - #添加投影 + 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 = gdal.Open(out_normal_tif, 1) ds.SetGeoTransform(geo) ds.SetProjection(proj) del ds diff --git a/plugins/veg_method/main.py b/plugins/veg_method/main.py index 18cc6b7..94b8034 100644 --- a/plugins/veg_method/main.py +++ b/plugins/veg_method/main.py @@ -3,7 +3,7 @@ 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 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 @@ -13,8 +13,9 @@ from thres import THRES from misc import table_layer, AlgSelectWidget from follow import FOLLOW + class VegtationCDMethod(QDialog): - def __init__(self,parent=None, alg:AlgFrontend=None): + def __init__(self, parent=None, alg: AlgFrontend = None): super(VegtationCDMethod, self).__init__(parent) self.alg = alg self.setWindowTitle('植被变化检测:{}'.format(alg.get_name())) @@ -23,11 +24,11 @@ class VegtationCDMethod(QDialog): 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 @@ -42,12 +43,12 @@ class VegtationCDMethod(QDialog): 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 = QDialogButtonBox(self) + buttonbox.addButton(self.ok_button, QDialogButtonBox.NoRole) + buttonbox.addButton(self.cancel_button, QDialogButtonBox.NoRole) buttonbox.setCenterButtons(True) - totalvlayout=QVBoxLayout() + totalvlayout = QVBoxLayout() totalvlayout.addLayout(layerbox) totalvlayout.addWidget(self.filter_select) if self.param_widget is not None: @@ -59,9 +60,10 @@ class VegtationCDMethod(QDialog): hbox.addWidget(buttonbox) totalvlayout.addLayout(hbox) # totalvlayout.addStretch() - + self.setLayout(totalvlayout) + @FOLLOW.register class VegetationCDFollow(AlgFrontend): @@ -71,14 +73,14 @@ class VegetationCDFollow(AlgFrontend): @staticmethod def get_icon(): - return IconInstance().UNSUPERVISED - + return IconInstance().VEGETATION + @staticmethod def get_widget(parent=None): widget = QWidget(parent) layer_combox = PairLayerCombox(widget) layer_combox.setObjectName('layer_combox') - + filter_select = AlgSelectWidget(widget, FILTER) filter_select.setObjectName('filter_select') unsupervised_select = AlgSelectWidget(widget, VEG_CD) @@ -86,83 +88,84 @@ class VegetationCDFollow(AlgFrontend): thres_select = AlgSelectWidget(widget, THRES) thres_select.setObjectName('thres_select') - totalvlayout=QVBoxLayout() + totalvlayout = QVBoxLayout() totalvlayout.addWidget(layer_combox) totalvlayout.addWidget(filter_select) totalvlayout.addWidget(unsupervised_select) totalvlayout.addWidget(thres_select) totalvlayout.addStretch() - + widget.setLayout(totalvlayout) return widget @staticmethod - def get_params(widget:QWidget=None): + def get_params(widget: QWidget = None): if widget is None: return dict() - + layer_combox = widget.findChild(PairLayerCombox, 'layer_combox') filter_select = widget.findChild(AlgSelectWidget, 'filter_select') - unsupervised_select = widget.findChild(AlgSelectWidget, 'unsupervised_select') + unsupervised_select = widget.findChild( + AlgSelectWidget, 'unsupervised_select') thres_select = widget.findChild(AlgSelectWidget, 'thres_select') - layer1=layer_combox.layer1 + layer1 = layer_combox.layer1 pth1 = layer_combox.layer1.path pth2 = layer_combox.layer2.path - falg, fparams = filter_select.get_alg_and_params() + falg, fparams = filter_select.get_alg_and_params() cdalg, cdparams = unsupervised_select.get_alg_and_params() thalg, thparams = thres_select.get_alg_and_params() - + if cdalg is None or thalg is None: return dict() return dict( layer1=layer1, - pth1 = pth1, - pth2 = pth2, - falg = falg, - fparams = fparams, - cdalg = cdalg, - cdparams = cdparams, - thalg = thalg, - thparams = thparams, + pth1=pth1, + pth2=pth2, + falg=falg, + fparams=fparams, + cdalg=cdalg, + cdparams=cdparams, + thalg=thalg, + thparams=thparams, ) @staticmethod def run_alg(layer1=None, - pth1 = None, - pth2 = None, - falg = None, - fparams = None, - cdalg = None, - cdparams = None, - thalg = None, - thparams = None, - send_message = None): - + pth1=None, + pth2=None, + falg=None, + fparams=None, + cdalg=None, + cdparams=None, + thalg=None, + thparams=None, + send_message=None): + if cdalg is None or thalg is None: return name = layer1.name if falg is not None: - pth1 = falg.run_alg(pth1, name=name, send_message= send_message, **fparams) - pth2 = falg.run_alg(pth2, name=name, send_message= send_message, **fparams) - - - cdpth = cdalg.run_alg(pth1, pth2, layer1.layer_parent, send_message= send_message,**cdparams) - thpth = thalg.run_alg(cdpth, name=name, send_message= send_message, **thparams) - - table_layer(thpth,layer1,name, send_message) + pth1 = falg.run_alg( + pth1, name=name, send_message=send_message, **fparams) + pth2 = falg.run_alg( + pth2, name=name, send_message=send_message, **fparams) + cdpth = cdalg.run_alg(pth1, pth2, layer1.layer_parent, + send_message=send_message, **cdparams) + thpth = thalg.run_alg( + cdpth, name=name, send_message=send_message, **thparams) + table_layer(thpth, layer1, name, send_message) class VegtationPlugin(BasicPlugin): - @staticmethod def info(): return { @@ -178,19 +181,18 @@ class VegtationPlugin(BasicPlugin): # ActionManager().veg_menu.addMenu(veg_menu) for key in VEG_CD.keys(): - alg:AlgFrontend = VEG_CD[key] + alg: AlgFrontend = VEG_CD[key] if alg.get_name() is None: name = key else: name = alg.get_name() - + action = QAction(name, veg_menu) func = partial(self.run_cd, alg) action.triggered.connect(func) veg_menu.addAction(action) - def run_cd(self, alg): dialog = VegtationCDMethod(self.mainwindow, alg) dialog.show() @@ -198,29 +200,31 @@ class VegtationPlugin(BasicPlugin): if dialog.exec_() == QDialog.Accepted: t = Thread(target=self.run_cd_alg, args=(dialog,)) t.start() - - def run_cd_alg(self, w:VegtationCDMethod): - - layer1=w.layer_combox.layer1 + + def run_cd_alg(self, w: VegtationCDMethod): + + 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() + 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) + 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) + 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) - \ No newline at end of file + table_layer(thpth, layer1, name, self.send_message) diff --git a/plugins/veg_method/scripts/__init__.py b/plugins/veg_method/scripts/__init__.py index 2572358..743c7d0 100644 --- a/plugins/veg_method/scripts/__init__.py +++ b/plugins/veg_method/scripts/__init__.py @@ -1,51 +1,60 @@ +from .SH import SH +from .LHBA import LHBA +from .OCD import OCD +from .AHT import AHT +from .ACD import ACD +import numpy as np from datetime import datetime from osgeo import gdal -import math,os +import math +import os import time from PyQt5 import QtWidgets from sklearn.cluster import k_means from rscder.utils.geomath import geo2imageRC, imageRC2geo +from rscder.utils.icons import IconInstance from rscder.utils.project import Project, PairLayer from misc import Register, AlgFrontend VEG_CD = Register('植被变化检测方法') -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]): +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 + 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: 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)) + for b in range(1, band+1): + out_ds.GetRasterBand(b).WriteArray( + ds.ReadAsArray(*srcWin, band_list=[b]), *(0, 0)) del out_ds + @VEG_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) - + def get_icon(): + return IconInstance().VEGETATION + + @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] @@ -56,48 +65,52 @@ class BasicCD(AlgFrontend): 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 + 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]) + 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]) - for j in range(yblocks + 1):#该改这里了 + 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]) + 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]) + 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]) + block_size1 = (xsize, end1y - block_xy1[1]) if block_xy2[1] + block_size2[1] > end2y: - block_size2 = (xsize, end2y - block_xy2[1]) + 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, ...] + 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) @@ -110,10 +123,12 @@ class BasicCD(AlgFrontend): 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) + 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) @@ -141,6 +156,7 @@ class BasicCD(AlgFrontend): send_message.emit('差分法计算完成') return out_normal_tif + @VEG_CD.register class LSTS(AlgFrontend): @@ -157,13 +173,13 @@ class LSTS(AlgFrontend): @staticmethod def get_params(widget=None): - return dict(n=5, w_size=(3,3)) - + 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) - + 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] @@ -172,84 +188,91 @@ class LSTS(AlgFrontend): 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_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 + 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] + 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) + 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) - 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]) + k_ = np.array(range(1, n+1)) + df1 = np.zeros(pixnum).astype(np.float64) + df2 = np.zeros(pixnum).astype(np.float64) - 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]) + 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]) + 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]) + 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]) + block_size1 = (xsize, end1y - block_xy1[1]) if block_xy2[1] + block_size2[1] > end2y: - block_size2 = (xsize, end2y - block_xy2[1]) + 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, ...] + 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) + 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 + 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) + 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 - 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) @@ -261,10 +284,12 @@ class LSTS(AlgFrontend): 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) + 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) @@ -288,7 +313,7 @@ class LSTS(AlgFrontend): os.remove(out_tif) except: pass - + if send_message is not None: send_message.emit('LSTS法计算完成') return out_normal_tif @@ -302,11 +327,11 @@ class CVAAlg(AlgFrontend): return 'CVA' @staticmethod - def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None, *args, **kargs): + 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) - 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] @@ -317,45 +342,48 @@ class CVAAlg(AlgFrontend): 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 + 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]) + 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]) + 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]) + 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]) + 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]) + block_size1 = (xsize, end1y - block_xy1[1]) if block_xy2[1] + block_size2[1] > end2y: - block_size2 = (xsize, end2y - block_xy2[1]) + 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, ...] + 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 + 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) @@ -367,10 +395,12 @@ class CVAAlg(AlgFrontend): 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) + 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) @@ -398,19 +428,20 @@ class CVAAlg(AlgFrontend): send_message.emit('欧式距离计算完成') return out_normal_tif + @VEG_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): - + 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() @@ -418,43 +449,45 @@ class ACDAlg(AlgFrontend): xsize = layer_parent.size[0] ysize = layer_parent.size[1] - geo=layer_parent.grid.geo - proj=layer_parent.grid.proj - #提取公共部分 + 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]) + 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) + 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]) + 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) - #添加投影 + 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 = gdal.Open(out_normal_tif, 1) ds.SetGeoTransform(geo) ds.SetProjection(proj) del ds - + return out_normal_tif @@ -464,56 +497,58 @@ 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): - + 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 - #提取公共部分 + 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]) + 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) + 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]) + 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) - #添加投影 + 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 = gdal.Open(out_normal_tif, 1) ds.SetGeoTransform(geo) ds.SetProjection(proj) del ds - + return out_normal_tif @@ -523,109 +558,115 @@ 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): - + 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 - #提取公共部分 + 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]) + 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) + 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]) + 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) - #添加投影 + 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 = gdal.Open(out_normal_tif, 1) ds.SetGeoTransform(geo) ds.SetProjection(proj) del ds - + return out_normal_tif + @VEG_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): - + 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 - #提取公共部分 + 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]) + 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) + 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]) + 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) - #添加投影 + 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 = gdal.Open(out_normal_tif, 1) ds.SetGeoTransform(geo) ds.SetProjection(proj) del ds @@ -638,52 +679,54 @@ 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): - + 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 - #提取公共部分 + 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]) + 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) + 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]) + 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) - #添加投影 + 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 = gdal.Open(out_normal_tif, 1) ds.SetGeoTransform(geo) ds.SetProjection(proj) del ds diff --git a/rscder/utils/icons.py b/rscder/utils/icons.py index 4643057..995201e 100644 --- a/rscder/utils/icons.py +++ b/rscder/utils/icons.py @@ -2,6 +2,7 @@ from PyQt5.QtGui import QIcon from PyQt5.QtCore import QObject from .misc import singleton + @singleton class IconInstance(QObject): @@ -10,9 +11,9 @@ class IconInstance(QObject): self.GRID_ON = QIcon('./icons/格网开.png') self.AI_DETECT = QIcon('./icons/AI变化检测.png') - self.EVALUATION = QIcon('./icons/精度评估.png') + self.EVALUATION = QIcon('./icons/精度评估.png') self.SPLASH = QIcon('./icons/splash.png') - self.HELP = QIcon('./icons/帮助.png') + self.HELP = QIcon('./icons/帮助.png') self.LOGO = QIcon('./icons/变化检测.png') self.CHANGE = QIcon('./icons/变化检测.png') self.OK = QIcon('./icons/ok.svg') @@ -46,18 +47,22 @@ class IconInstance(QObject): self.RASTER = QIcon('./icons/影像.png') self.VEGETATION = QIcon('./icons/植被变化.png') self.NOISE = QIcon('./icons/噪声处理.png') + self.ARITHMETIC1 = QIcon('./icons/Algorithm_icon/功能-01.png') + self.ARITHMETIC2 = QIcon('./icons/Algorithm_icon/功能-02.png') + self.ARITHMETIC3 = QIcon('./icons/Algorithm_icon/功能-03.png') + self.ARITHMETIC3.QSize() self.DATA_LOAD = QIcon('./icons/数据加载.png') - + self.EXCIT = QIcon('./icons/退出.png') self.ZOOM_IN = QIcon('./icons/放大.png') self.ZOOM_OUT = QIcon('./icons/缩小.png') self.TABLE = QIcon('./icons/table.png') - + def __getattr__(self, name: str) -> QIcon: try: return getattr(self, name) except: - return self.LOGO \ No newline at end of file + return self.LOGO diff --git a/test-data/AAA.tif.ovr b/test-data/AAA.tif.ovr new file mode 100644 index 0000000..2fc0fa2 Binary files /dev/null and b/test-data/AAA.tif.ovr differ diff --git a/test-data/BBB.tif.ovr b/test-data/BBB.tif.ovr new file mode 100644 index 0000000..eb7f4fc Binary files /dev/null and b/test-data/BBB.tif.ovr differ diff --git a/test-data/test01/test01.prj b/test-data/test01/test01.prj new file mode 100644 index 0000000..37cf7ac --- /dev/null +++ b/test-data/test01/test01.prj @@ -0,0 +1,22 @@ +cell_size: +- 100 +- 100 +layers: +- enable: true + layers: [] + name: AAA.-BBB. + pth1: E:/xhong/CVEOdemo/fourth/rscder/test-data/AAA.tif + pth2: E:/xhong/CVEOdemo/fourth/rscder/test-data/BBB.tif + style_info1: + NIR: 4 + b: 3 + g: 2 + r: 1 + style_info2: + NIR: 4 + b: 3 + g: 2 + r: 1 +max_memory: 100 +max_threads: 4 +root: E:\xhong\CVEOdemo\fourth\rscder\test-data\test01