from rscder.utils.geomath import geo2imageRC, imageRC2geo from rscder.utils.project import Project, PairLayer from rscder.utils.icons import IconInstance from misc import Register, AlgFrontend from osgeo import gdal import os import math import numpy as np from PyQt5.QtGui import QIntValidator from PyQt5.QtWidgets import QWidget, QLabel, QPushButton, QLineEdit, QHBoxLayout, QVBoxLayout from . import VEG_CD @VEG_CD.register class VFCCD(AlgFrontend): @staticmethod def get_name(): return '植被覆盖度变化' @staticmethod def get_icon(): return IconInstance().ARITHMETIC2 @staticmethod def get_widget(parent=None): widget = QWidget(parent) layout = QVBoxLayout() redlabel = QLabel('红波段：') redinput = QLineEdit() redinput.setValidator(QIntValidator(redinput)) redinput.setObjectName('redinput') redlayout = QHBoxLayout() redlayout.addWidget(redlabel) redlayout.addWidget(redinput) nirlabel = QLabel('近红外：') nirinput = QLineEdit() nirinput.setValidator(QIntValidator(0, 100, nirinput)) nirinput.setObjectName('nirinput') nirlayout = QHBoxLayout() nirlayout.addWidget(nirlabel) nirlayout.addWidget(nirinput) layout.addLayout(redlayout) layout.addLayout(nirlayout) widget.setLayout(layout) return widget @staticmethod def get_params(widget:QWidget=None): if widget is None: return dict(nir=3, red=0) redinput = widget.findChild(QLineEdit, 'redinput') nirinput = widget.findChild(QLineEdit, 'nirinput') if redinput is None or nirinput is None: return dict(nir=3, red=0) nir = int(nirinput.text()) red = int(redinput.text()) return dict(red=red, nir=nir) @staticmethod def run_alg(pth1:str,pth2:str,layer_parent:PairLayer, red=0, nir=3, send_message = None,*args, **kargs): ds1:gdal.Dataset=gdal.Open(pth1) ds2:gdal.Dataset=gdal.Open(pth2) if ds1.RasterCount <= 2: if send_message is not None: send_message.emit('至少包含两个波段') return if ds1.RasterCount < (max(red, nir) + 1): if send_message is not None: send_message.emit(f'{max(nir, red)}超过波段数量') return 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) ndvi1 = (block_data1[nir] - block_data1[red])/((block_data1[nir] + block_data1[red]) + 1e-6) ndvi2 = (block_data2[nir] - block_data2[red])/((block_data2[nir] + block_data2[red]) + 1e-6) ndvi1[ndvi1 < 0] = 0 ndvi2[ndvi2 < 0] = 0 # pdb.set_trace() block_diff = ndvi1 - ndvi2 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