Add files via upload

sen+mk/mk.py

@@ -0,0 +1,136 @@
+'''
+Created on 2021年4月11日
+
+@author: SunStrong
+'''
+#coding:utf-8
+import numpy as np
+import pymannkendall as mk
+import os 
+import rasterio as ras
+
+
+
+def sen_mk_test(path1,outputPath):
+
+    #image_path:影像的存储路径
+    #outputPath:结果输出路径
+
+    filepaths=[]
+    for file in os.listdir(path1):
+        filepath1=os.path.join(path1,file)
+        filepaths.append(filepath1)
+
+    #获取影像数量
+    num_images=len(filepaths)
+    #读取影像数据
+    img1=ras.open(filepaths[0])
+    #获取影像的投影，高度和宽度
+    transform1=img1.transform
+    height1=img1.height
+    width1=img1.width 
+    array1=img1.read()
+    img1.close()
+
+    #读取所有影像
+    for path1 in filepaths[1:]:
+        if path1[-4:]=='tiff':
+            print(path1)
+            img2=ras.open(path1)
+            array2=img2.read()
+            array1=np.vstack((array1,array2))
+            img2.close()
+
+    nums,width,height=array1.shape 
+    #写影像
+    def writeImage(image_save_path,height1,width1,para_array,bandDes,transform1):
+        with ras.open(
+               image_save_path,
+               'w',
+               driver='GTiff',
+               height=height1,
+               width=width1,
+               count=1,
+               dtype=para_array.dtype,
+               crs='+proj=latlong',
+               transform=transform1,
+        ) as dst:
+                   dst.write_band(1,para_array)
+                   dst.set_band_description(1,bandDes)
+        del dst
+
+    #输出矩阵，无值区用-9999填充    
+    slope_array=np.full([width,height],-9999.0000) 
+    z_array=np.full([width,height],-9999.0000)
+    Trend_array=np.full([width,height],-9999.0000) 
+    Tau_array=np.full([width,height],-9999.0000)
+    s_array=np.full([width,height],-9999.0000)
+    p_array=np.full([width,height],-9999.0000)
+    #只有有值的区域才进行mk检验
+    c1=np.isnan(array1)
+    sum_array1=np.sum(c1,axis=0)
+    nan_positions=np.where(sum_array1==num_images)
+
+    positions=np.where(sum_array1!=num_images) 
+
+
+    #输出总像元数量
+    print("all the pixel counts are {0}".format(len(positions[0])))
+    #mk test
+    for i in range(len(positions[0])):
+        print(i)
+        x=positions[0][i]
+        y=positions[1][i]    
+        mk_list1=array1[:,x,y]
+        trend, h, p, z, Tau, s, var_s, slope, intercept  = mk.original_test(mk_list1)
+        '''        
+        trend: tells the trend (increasing, decreasing or no trend)
+                h: True (if trend is present) or False (if trend is absence)
+                p: p-value of the significance test
+                z: normalized test statistics
+                Tau: Kendall Tau
+                s: Mann-Kendal's score
+                var_s: Variance S
+                slope: Theil-Sen estimator/slope
+                intercept: intercept of Kendall-Theil Robust Line
+        '''
+
+
+        if trend=="decreasing":
+            trend_value=-1
+        elif trend=="increasing":
+            trend_value=1
+        else:
+            trend_value=0
+        slope_array[x,y]=slope#senslope
+        s_array[x,y]=s
+        z_array[x,y]=z
+        Trend_array[x,y]=trend_value
+        p_array[x,y]=p
+        Tau_array[x,y]=Tau 
+
+
+    all_array=[slope_array,Trend_array,p_array,s_array,Tau_array,z_array]   
+
+    slope_save_path=os.path.join(result_path,"slope.tiff")
+    Trend_save_path=os.path.join(result_path,"Trend.tiff")
+    p_save_path=os.path.join(result_path,"p.tiff")
+    s_save_path=os.path.join(result_path,"s.tiff")
+    tau_save_path=os.path.join(result_path,"tau.tiff")
+    z_save_path=os.path.join(result_path,"z.tiff")
+    image_save_paths=[slope_save_path,Trend_save_path,p_save_path,s_save_path,tau_save_path,z_save_path]
+    band_Des=['slope','trend','p_value','score','tau','z_value']
+    for i in range(len(all_array)):
+        writeImage(image_save_paths[i], height1, width1, all_array[i], band_Des[i],transform1)
+
+#调用
+if __name__ == '__main__':
+
+    path1="H:\\001zone\\00张家口承德地区\\cal_zc_evf\\evf_year"
+    result_path="H:\\001zone\\00张家口承德地区\\cal_zc_evf\\result"
+    sen_mk_test(path1, result_path)
+
+
+
+
+

sen+mk/month_bands.py

@@ -0,0 +1,118 @@
+
+# -*- coding: utf-8 -*-
+import glob
+
+from osgeo import gdal
+import numpy as np
+import os
+
+'''
+说明：
+    提取img数据中的红、绿、近红外波段数据
+        Rrs_555	绿波段（62）
+        Rrs_645	红波段（63）
+        Rrs_859	近红外波段（64）
+        注：下划线后的数字为对应波段的中心波长
+'''
+
+def writeTiff(im_data, im_width, im_height, im_bands, im_geotrans, im_proj, path):  # 定义函数writeTiff
+    if 'int8' in im_data.dtype.name:
+        datatype = gdal.GDT_Byte
+    elif 'int16' in im_data.dtype.name:
+        datatype = gdal.GDT_UInt16
+    else:
+        datatype = gdal.GDT_Float32
+
+    if len(im_data.shape) == 3:
+        im_bands, im_height, im_width = im_data.shape
+    elif len(im_data.shape) == 2:
+        im_data = np.array([im_data])
+    else:
+        im_bands, (im_height, im_width) = 1, im_data.shape
+    # 创建文件
+    driver = gdal.GetDriverByName("GTiff")
+    dataset = driver.Create(path, im_width, im_height, im_bands, datatype)
+    if (dataset != None):
+        dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
+        dataset.SetProjection(im_proj)  # 写入投影
+    for i in range(im_bands):
+        dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
+    del dataset
+
+
+if __name__ == '__main__':
+    gdal.AllRegister()
+
+    year = 2002
+    while year < 2021:
+        path_mother = os.path.join('H:\\001zone\\00张家口承德地区\\cal_zc_evf\\evf_bands',str(year)) # 输入路径，文件夹
+        path_out = os.path.join('H:\\001zone\\00张家口承德地区\\cal_zc_evf\\month_bands',str(year))  # 输出路径，文件夹
+
+        # 获取该目录下所有文件，存入列表中
+        os.chdir(path_mother)
+        for file in glob.glob('*.tiff'):
+            img_file_name = path_mother + os.sep + file
+            dataset = gdal.Open(img_file_name)
+
+            filename = str(file.split('_')[0])
+            print(filename)
+
+            month = int(filename[4:6])
+            adfGeoTransform = dataset.GetGeoTransform()
+
+            band_1 = dataset.GetRasterBand(1)   # 读取数据集的属性，选择待提取波段.62,63,64代表的是波段其在img中的序列号（绿波段在img影像中位于第62个波段）
+            #band_2 = dataset.GetRasterBand(63)
+            #band_3 = dataset.GetRasterBand(64)
+
+            im_width = dataset.RasterXSize  # 栅格矩阵的列数
+            im_height = dataset.RasterYSize  # 栅格矩阵的行数
+
+            Rrs_1 = band_1.ReadAsArray(0, 0, im_width, im_height)  # 目标波段
+            #Rrs_2 = band_2.ReadAsArray(0, 0, im_width, im_height)
+            #Rrs_3 = band_3.ReadAsArray(0, 0, im_width, im_height)
+
+            im_geotrans = dataset.GetGeoTransform()  # 获取仿射矩阵信息
+            im_proj = dataset.GetProjection()  # 获取投影信息
+
+            im_bands = 1  # band_469.RasterCount  # 输出文件的波段数
+            if month == 1:
+                B1 = Rrs_1 * 31  # / 10000
+            elif month == 2:
+                B1 = Rrs_1 * 28
+            elif month == 3:
+                B1 = Rrs_1 * 31
+            elif month == 4:
+                B1 = Rrs_1 * 30
+            elif month == 5:
+                B1 = Rrs_1 * 31
+            elif month == 6:
+                B1 = Rrs_1 * 30
+            elif month == 7:
+                B1 = Rrs_1 * 31
+            elif month == 8:
+                B1 = Rrs_1 * 31
+            elif month == 9:
+                B1 = Rrs_1 * 30
+            elif month == 10:
+                B1 = Rrs_1 * 31
+            elif month == 11:
+                B1 = Rrs_1 * 30
+            elif month == 12:
+                B1 = Rrs_1 * 31
+
+            #B2 = Rrs_2 * 1.0  # / 10000
+            #B3 = Rrs_3 * 1.0  # / 10000
+
+            outpath_1 = path_out + os.sep + filename  + '.tiff'
+            writeTiff(B1, im_width, im_height, im_bands, im_geotrans, im_proj, outpath_1)
+
+            #outpath_2 = path_out + os.sep + filename + '_Rrs_645.Red' + '.tif'
+            #writeTiff(B2, im_width, im_height, im_bands, im_geotrans, im_proj, outpath_2)
+
+            #outpath_3 = path_out + os.sep + filename + '_Rrs_859.NIR' + '.tif'
+            #writeTiff(B3, im_width, im_height, im_bands, im_geotrans, im_proj, outpath_3)
+            print('this is ok')
+        year += 1
+        print(year)
+    print('------All is ok------')
+

sen+mk/year_bands.py

@@ -0,0 +1,119 @@
+
+# -*- coding: utf-8 -*-
+import glob
+
+from osgeo import gdal
+import numpy as np
+import os
+
+'''
+说明：
+    提取img数据中的红、绿、近红外波段数据
+        Rrs_555	绿波段（62）
+        Rrs_645	红波段（63）
+        Rrs_859	近红外波段（64）
+        注：下划线后的数字为对应波段的中心波长
+'''
+
+def writeTiff(im_data, im_width, im_height, im_bands, im_geotrans, im_proj, path):  # 定义函数writeTiff
+    if 'int8' in im_data.dtype.name:
+        datatype = gdal.GDT_Byte
+    elif 'int16' in im_data.dtype.name:
+        datatype = gdal.GDT_UInt16
+    else:
+        datatype = gdal.GDT_Float32
+
+    if len(im_data.shape) == 3:
+        im_bands, im_height, im_width = im_data.shape
+    elif len(im_data.shape) == 2:
+        im_data = np.array([im_data])
+    else:
+        im_bands, (im_height, im_width) = 1, im_data.shape
+    # 创建文件
+    driver = gdal.GetDriverByName("GTiff")
+    dataset = driver.Create(path, im_width, im_height, im_bands, datatype)
+    if (dataset != None):
+        dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
+        dataset.SetProjection(im_proj)  # 写入投影
+    for i in range(im_bands):
+        dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
+    del dataset
+
+
+if __name__ == '__main__':
+    gdal.AllRegister()
+
+    year = 2002
+    while year < 2003:
+        path_mother = os.path.join('H:\\001zone\\00张家口承德地区\\cal_zc_evf\\month_bands',str(year)) # 输入路径，文件夹
+        #path_out = os.path.join('H:\\001zone\\00张家口承德地区\\cal_zc_evf\\month_bands',str(year))  # 输出路径，文件夹
+        path_out = ('H:\\001zone\\00张家口承德地区\\cal_zc_evf\\evf_year') # 输出路径，文件夹
+
+        # 获取该目录下所有文件，存入列表中
+        os.chdir(path_mother)
+        for file in glob.glob('*.tiff'):
+            img_file_name = path_mother + os.sep + file
+            dataset = gdal.Open(img_file_name)
+
+            filename = str(file.split('_')[0])
+            print(filename)
+
+            month = int(filename[4:6])
+            adfGeoTransform = dataset.GetGeoTransform()
+
+            band_1 = dataset.GetRasterBand(1)   # 读取数据集的属性，选择待提取波段.62,63,64代表的是波段其在img中的序列号（绿波段在img影像中位于第62个波段）
+            #band_2 = dataset.GetRasterBand(63)
+            #band_3 = dataset.GetRasterBand(64)
+
+            im_width = dataset.RasterXSize  # 栅格矩阵的列数
+            im_height = dataset.RasterYSize  # 栅格矩阵的行数
+
+            Rrs_1 = band_1.ReadAsArray(0, 0, im_width, im_height)  # 目标波段
+            #Rrs_2 = band_2.ReadAsArray(0, 0, im_width, im_height)
+            #Rrs_3 = band_3.ReadAsArray(0, 0, im_width, im_height)
+
+            im_geotrans = dataset.GetGeoTransform()  # 获取仿射矩阵信息
+            im_proj = dataset.GetProjection()  # 获取投影信息
+
+            im_bands = 1  # band_469.RasterCount  # 输出文件的波段数
+            if month == 1:
+                B1 = Rrs_1 * 1  # / 10000
+            elif month == 2:
+                B2 = Rrs_1 * 1
+            elif month == 3:
+                B3 = Rrs_1 * 1
+            elif month == 4:
+                B4 = Rrs_1 * 1
+            elif month == 5:
+                B5 = Rrs_1 * 1
+            elif month == 6:
+                B6 = Rrs_1 * 1
+            elif month == 7:
+                B7 = Rrs_1 * 1
+            elif month == 8:
+                B8 = Rrs_1 * 1
+            elif month == 9:
+                B9 = Rrs_1 * 1
+            elif month == 10:
+                B10 = Rrs_1 * 1
+            elif month == 11:
+                B11 = Rrs_1 * 1
+            elif month == 12:
+                B12 = Rrs_1 * 1
+
+            #B2 = Rrs_2 * 1.0  # / 10000
+            #B3 = Rrs_3 * 1.0  # / 10000
+        B13 = B1 + B2 + B3 + B4 + B5 + B6 + B7 + B8 + B9 + B10 + B11 + B12
+        outpath_1 = path_out + os.sep + str(year)  + '_sum_year_evf.tiff'
+        writeTiff(B13, im_width, im_height, im_bands, im_geotrans, im_proj, outpath_1)
+
+            #outpath_2 = path_out + os.sep + filename + '_Rrs_645.Red' + '.tif'
+            #writeTiff(B2, im_width, im_height, im_bands, im_geotrans, im_proj, outpath_2)
+
+            #outpath_3 = path_out + os.sep + filename + '_Rrs_859.NIR' + '.tif'
+            #writeTiff(B3, im_width, im_height, im_bands, im_geotrans, im_proj, outpath_3)
+        print('this is ok')
+        year += 1
+        print(year)
+    print('------All is ok------')
+