179 lines
6.6 KiB
Python
179 lines
6.6 KiB
Python
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 |