diff --git a/leafAreaIndex/product.xml b/leafAreaIndex/product.xml new file mode 100644 index 00000000..017cb135 --- /dev/null +++ b/leafAreaIndex/product.xml @@ -0,0 +1,60 @@ + + + 后向散射系数 + BackScattering + 5 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + None + 参考产品介绍PDF + + + + + 德清 + + + + \ No newline at end of file diff --git a/leafAreaIndex/sample_process.py b/leafAreaIndex/sample_process.py new file mode 100644 index 00000000..c32132e5 --- /dev/null +++ b/leafAreaIndex/sample_process.py @@ -0,0 +1,202 @@ +# +# 样本处理的相关的库 +# + +from tool.algorithm.image.ImageHandle import ImageHandler +import math +import numpy as np +import random +import scipy +# 最小二乘求解非线性方程组 +from scipy.optimize import leastsq,fsolve,root +from osgeo import gdal,gdalconst +import pandas as pds +from scipy import interpolate +from multiprocessing import pool +# 常量声明区域 +imageHandler=ImageHandler() + + +# python 的函数类 +def read_sample_csv(csv_path): + """ 读取样本的csv + Args: + csv_path (string): 样本csv的地址,绝对路径 + return: + [ + ['日期','样方编号','经度','纬度','叶面积指数',"后向散射系数"], + ['日期','样方编号','经度','纬度','叶面积指数',"后向散射系数"],...... + ] + """ + lai_csv=pds.read_csv(csv_path)# 代码测试区域 + lai_csv=lai_csv.loc[:,['id','lon','lat','leaf',"cal"]] + result=[] + for i in range(len(lai_csv)): + result.append([ + 0, + lai_csv.loc[i,'id'], + lai_csv.loc[i,'lon'], # lon,x + lai_csv.loc[i,'lat'], # lat,y + lai_csv.loc[i,'leaf'], + 10**(float(lai_csv.loc[i,'cal'])/10), + ]) + return result +def read_tiff(tiff_path): + """ 从文件中读取影像 + + Args: + tiff_path (string): 文件影像路径 + """ + im_proj, im_geotrans, im_arr=imageHandler.read_img(tiff_path) + return { + 'proj':im_proj, + 'geotrans':im_geotrans, + 'data':im_arr + } + +def ReprojectImages2(in_tiff_path,ref_tiff_path,out_tiff_path,resampleAlg=gdalconst.GRA_Bilinear): + """ 将输入影像重采样到参考影像的范围内 + + Args: + in_tiff_path (string): 输入影像 + ref_tiff_path (string): 参考影像 + out_tiff_path (string): 输出地址 + resampleAlg (gadlconst): 插值方法 + """ + # 若采用gdal.Warp()方法进行重采样 + # 获取输出影像信息 + inputrasfile = gdal.Open(in_tiff_path, gdal.GA_ReadOnly) + inputProj = inputrasfile.GetProjection() + # 获取参考影像信息 + referencefile = gdal.Open(ref_tiff_path, gdal.GA_ReadOnly) + referencefileProj = referencefile.GetProjection() + referencefileTrans = referencefile.GetGeoTransform() + bandreferencefile = referencefile.GetRasterBand(1) + x = referencefile.RasterXSize + y = referencefile.RasterYSize + nbands = referencefile.RasterCount + # 创建重采样输出文件(设置投影及六参数) + driver = gdal.GetDriverByName('GTiff') + output = driver.Create(out_tiff_path, x, y, nbands, bandreferencefile.DataType) + output.SetGeoTransform(referencefileTrans) + output.SetProjection(referencefileProj) + options = gdal.WarpOptions(srcSRS=inputProj, dstSRS=referencefileProj, resampleAlg=gdalconst.GRA_Bilinear) + gdal.Warp(output, in_tiff_path, options=options) + +def combine_sample_attr(sample_list,attr_tiff): + """ 构建样本 + + Args: + sample_list (list): 原样本 + attr_tiff (string): 添加的属性数据 + + Returns: + list:[sample,new_attr] + """ + result=[] + # 因为soil_tiff 的影像的 影像分辨率较低 + inv_gt=gdal.InvGeoTransform(attr_tiff['geotrans']) + for sample_item in sample_list: + sample_lon=sample_item[2] + sample_lat=sample_item[3] + sample_in_tiff_x=inv_gt[0]+inv_gt[1]*sample_lon+inv_gt[2]*sample_lat # x + sample_in_tiff_y=inv_gt[3]+inv_gt[4]*sample_lon+inv_gt[5]*sample_lat # y + x_min=int(np.floor(sample_in_tiff_x)) + x_max=int(np.ceil(sample_in_tiff_x)) + y_min=int(np.floor(sample_in_tiff_y)) + y_max=int(np.ceil(sample_in_tiff_y)) + if x_min<0 or y_min<0 or x_max>=attr_tiff['data'].shape[1] or y_max>=attr_tiff['data'].shape[0]: + continue + # + """ + f = interpolate.interp2d([0,0,1,1], [0,1,1,0], + [attr_tiff['data'][y_min,x_min], + attr_tiff['data'][y_max,x_min], + attr_tiff['data'][y_max,x_max], + attr_tiff['data'][y_min,x_min] + ], kind='linear') + interp_value=f(sample_in_tiff_x-x_min,sample_in_tiff_y-y_min) + sample_item.append(interp_value[0]) + """ + # 9x9 + x_min=x_min-4 if x_min-9>=0 else 0 + y_min=y_min-4 if y_min-9>=0 else 0 + x_max=x_max+4 if x_max+490: + continue + if sample_soil<=0 or sample_soil>=1: + continue + if sample_lai<=0 or sample_lai>=20: + continue + result.append(item) + # 绘制分布图 + # lai=[] + # sigma=[] + # csv_sigmas=[] + # text_label=[] + # for item in result: + # if len(item)==10: + # sample_time,sample_code,sample_lon,sample_lat,sample_lai,csv_sigma,sample_soil,sample_inc,sample_sigma,sample_NDVI=item + # else: + # sample_time,sample_code,sample_lon,sample_lat,sample_lai,csv_sigma,sample_soil,sample_inc,sample_sigma=item + # text_label.append(sample_code) + # lai.append(sample_lai) + # sigma.append(sample_sigma) + # csv_sigmas.append(csv_sigma) + # from matplotlib import pyplot as plt + # plt.scatter(np.array(lai),np.array(sigma),label="lai-tiff_sigma") + # for i in range(len(sigma)): + # plt.annotate(text_label[i], xy = (lai[i], sigma[i])) # 这里xy是需要标记的坐标,xytext是对应的标签坐标 + # + # plt.scatter(np.array(lai),np.array(csv_sigmas),label="lai-csv_sigmas") + # for i in range(len(csv_sigmas)): + # plt.annotate(text_label[i], xy = (lai[i],csv_sigmas[i])) # 这里xy是需要标记的坐标,xytext是对应的标签坐标 + # plt.legend() + # plt.show() + + return result + + + +def split_sample_list(sample_list,train_ratio): + """ 切分样本比值 + + Args: + sample_list (list): 样本列表 + train_ratio (double): 训练样本的比重 + + Returns: + list: [sample_train,sample_test] + """ + sample_train=[] + sample_test=[] + n=len(sample_list) + for i in range(n): + if random.random()<=train_ratio: + sample_train.append(sample_list[i]) + else: + sample_test.append(sample_list[i]) + return [sample_train,sample_test] \ No newline at end of file diff --git a/soilMoistureTop/SoilMoistureMain.py b/soilMoistureTop/SoilMoistureMain.py index 7bb6512c..871e112b 100644 --- a/soilMoistureTop/SoilMoistureMain.py +++ b/soilMoistureTop/SoilMoistureMain.py @@ -249,7 +249,7 @@ class MoistureMain: bp = BlockProcess() # block_size = bp.get_block_size(self.__rows, self.__cols,block_size_config) block_size = bp.get_block_size(self.__rows, self.__cols) - bp.cut(self.__workspace_preprocessed_path, self.__workspace_block_tif_path, ['tif', 'tiff'], 'tif', block_size) + bp.cut_new(self.__workspace_preprocessed_path, self.__workspace_block_tif_path, ['tif', 'tiff'], 'tif', block_size) logger.info('blocking tifs success!') img_dir, img_name = bp.get_file_names(self.__workspace_block_tif_path, ['tif']) @@ -339,7 +339,7 @@ class MoistureMain: h = ImageHandler.get_img_height(self.__preprocessed_paras['HH']) out_path = self.__workspace_processing_path[0:-1] tif_type = bp.get_tif_dtype(self.__preprocessed_paras['HH']) - bp.combine(data_dir, w, h, out_path, file_type=['tif'], datetype=tif_type) + bp.combine_new(data_dir, w, h, out_path, file_type=['tif'], datetype=tif_type) # 添加地理信息 soil_moisture_path = self.__workspace_processing_path + 'soil_moisture.tif' diff --git a/soilMoistureTop/product.xml b/soilMoistureTop/product.xml new file mode 100644 index 00000000..017cb135 --- /dev/null +++ b/soilMoistureTop/product.xml @@ -0,0 +1,60 @@ + + + 后向散射系数 + BackScattering + 5 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + None + 参考产品介绍PDF + + + + + 德清 + + + + \ No newline at end of file diff --git a/soilSalinity/lee_refined_filter_T3.exe b/soilSalinity/lee_refined_filter_T3.exe new file mode 100644 index 00000000..3b4839c1 Binary files /dev/null and b/soilSalinity/lee_refined_filter_T3.exe differ diff --git a/soilSalinity/product.xml b/soilSalinity/product.xml new file mode 100644 index 00000000..017cb135 --- /dev/null +++ b/soilSalinity/product.xml @@ -0,0 +1,60 @@ + + + 后向散射系数 + BackScattering + 5 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + None + 参考产品介绍PDF + + + + + 德清 + + + + \ No newline at end of file diff --git a/surfaceRoughness_oh2004/product.xml b/surfaceRoughness_oh2004/product.xml new file mode 100644 index 00000000..017cb135 --- /dev/null +++ b/surfaceRoughness_oh2004/product.xml @@ -0,0 +1,60 @@ + + + 后向散射系数 + BackScattering + 5 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + None + 参考产品介绍PDF + + + + + 德清 + + + + \ No newline at end of file diff --git a/vegetationPhenology/product.xml b/vegetationPhenology/product.xml new file mode 100644 index 00000000..017cb135 --- /dev/null +++ b/vegetationPhenology/product.xml @@ -0,0 +1,60 @@ + + + 后向散射系数 + BackScattering + 5 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 + None + 参考产品介绍PDF + + + + + 德清 + + + + \ No newline at end of file diff --git a/vegetationPhenology/testxmlreading.py b/vegetationPhenology/testxmlreading.py new file mode 100644 index 00000000..cc3a0377 --- /dev/null +++ b/vegetationPhenology/testxmlreading.py @@ -0,0 +1,79 @@ +#encoding=utf-8 +import xml.etree.ElementTree as ET +import pandas as pd +import csv + +def xml2csv(xmlpath): + tree_obj = ET.parse(xmlpath) + + # 得到所有匹配Region 标签的Element对象的list集合 + list_Region = tree_obj.findall("Region") + for Region in list_Region: + # 几何面对应的类(phenology_name)在标签 + Region_dict = Region.attrib + phenology_name = Region_dict.get("name") + print(phenology_name) + list_GeometryDef = Region.findall("GeometryDef") + list_Polygon = list_GeometryDef[0].findall("Polygon") # 获得该类下的几何面list + for polygon in list_Polygon: + # 从polygon list中获取得到标签的数据 注意是空格分隔的和csv中不同 + Coordinates_list = coordinates = polygon.find('.//Coordinates').text.strip().split() + # POLYGON((119.035 31.51,119.035 31.50,119.033 31.50)) csv中 + print("value") + +# 向csv中写 +def csvfile(csvpath,data): + + with open(csvpath, 'a', newline='') as file: + # 2. step + writer = csv.writer(file) + # data example + #data = ["This", "is", "a", "Test"] + writer.writerow(data) + + +# Define the structure of the data +#data示例 student_header = ['name', 'age', 'major', 'minor'] +def csvcreateTitile(csvpath,data): + # 1. Open a new CSV file + with open(csvpath, 'w') as file: + # 2. Create a CSV writer + writer = csv.writer(file) + # 3. Write data to the file + writer.writerow(data) + +# 将列表中的坐标对转换为字符串 +def createcsv_roi_polygon(coordinates): + coord_str = ', '.join([f'{coordinates[i]} {coordinates[i + 1]}' for i in range(0, len(coordinates), 2)]) + # 构建最终的POLYGON字符串 + polygon_str = f'POLYGON(({coord_str}))' + print(polygon_str) + return polygon_str + +if __name__ == '__main__': + xmlpath = r"E:\MicroWorkspace\GF3A_nanjing\input-ortho\test_shp\test.xml" + + tree_obj = ET.parse(xmlpath) + csv_header = ['sar_img_name', 'phenology_id', 'phenology_name', 'roi_polygon'] + csvpath = r"E:\MicroWorkspace\GF3A_nanjing\input-ortho\test_shp\test.csv" + # csvcreateTitile(csvpath,csv_header) + csvfile(csvpath,csv_header) + # 得到所有匹配Region 标签的Element对象的list集合 + list_Region = tree_obj.findall("Region") + for Region in list_Region: + # 几何面对应的类(phenology_name)在标签 + Region_dict = Region.attrib + phenology_name = Region_dict.get("name") + print(phenology_name) + # list_GeometryDef = Region.findall("GeometryDef") + list_Polygon = Region.findall(".//Polygon") # 获得该类下的几何面list + count = 1 + for polygon in list_Polygon: + # 从polygon list中获取得到标签的数据 注意是空格分隔的和csv中不同 + Coordinates_list = coordinates = polygon.find('.//Coordinates').text.strip().split() + # 将坐标和ploygon对应的写入到.csv中 + polygon_str=createcsv_roi_polygon(Coordinates_list) + # POLYGON((119.035 31.51,119.035 31.50,119.033 31.50)) csv中 + data = ['0', count, phenology_name, polygon_str] + csvfile(csvpath,data) + count += 1