remove unused method

plugins/veg_method/scripts/__init__.py

@@ -10,35 +10,14 @@ 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]):
-    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
 
 @VEG_CD.register
 class BasicCD(AlgFrontend):
 
     @staticmethod
     def get_name():
-        return '差分法'
+        return '植被覆盖度变化'
     
     @staticmethod
     def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None,*args, **kargs):

plugins/veg_method/scripts/vfc.py

@@ -0,0 +1,113 @@
+from rscder.utils.geomath import geo2imageRC, imageRC2geo
+from rscder.utils.project import Project, PairLayer
+from misc import Register, AlgFrontend
+
+from . import VEG_CD
+
+@VEG_CD.register
+class VFCCD(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