From e032838c399f536fffabf3515ebb3b4b1ac659ab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E9=99=88=E5=A2=9E=E8=BE=89?= <3045316072@qq.com> Date: Mon, 22 Sep 2025 15:07:12 +0800 Subject: [PATCH] =?UTF-8?q?=E6=8B=86=E5=88=86=E6=B8=AF=E5=8F=A3=E5=B7=A5?= =?UTF-8?q?=E5=85=B7?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- SplitPortRasterTools/SplitPortRasterTools.py | 0 .../SplitShipPortRasterTools.py | 287 ++++++++++++++++++ 2 files changed, 287 insertions(+) delete mode 100644 SplitPortRasterTools/SplitPortRasterTools.py create mode 100644 SplitPortRasterTools/SplitShipPortRasterTools.py diff --git a/SplitPortRasterTools/SplitPortRasterTools.py b/SplitPortRasterTools/SplitPortRasterTools.py deleted file mode 100644 index e69de29..0000000 diff --git a/SplitPortRasterTools/SplitShipPortRasterTools.py b/SplitPortRasterTools/SplitShipPortRasterTools.py new file mode 100644 index 0000000..a52a00f --- /dev/null +++ b/SplitPortRasterTools/SplitShipPortRasterTools.py @@ -0,0 +1,287 @@ +from osgeo import ogr +import os +import argparse +from osgeo import ogr, gdal +import os +import argparse +import numpy as np +from scipy.spatial import KDTree +from tools.DotaOperator import DotaObj,readDotaFile,writerDotaFile,createDota +from glob import glob +from pathlib import Path +import shutil + +""" +1. 港口，5000x5000，军港与民港切片， + a. 预标注框给他们 + b. 重叠率25%重复 + c. 军民2000m 以内，影像抽出来 + d. 原影像如果只有 1 个港口，-》 港口类型 + 多 个港口，-》 远近，如果港口距离1000m ,单独拉出来 + +""" + + +MLCName="MLC" +JLCName="JLC" +MJLCName="MJLC" +NOLCName="NOLC" + +def find_tif_files_pathlib(directory): + path = Path(directory) + # 使用rglob递归匹配所有.tif和.tiff文件 + tif_files = list(path.rglob('*.tiff'))+list(path.rglob('*.tif')) + # 将Path对象转换为字符串路径 + return [str(file) for file in tif_files] + + +def find_srcPath(srcFolder): + root_path = Path(srcFolder) + target_path = [folderpath for folderpath in root_path.rglob("*") if folderpath.is_dir() and folderpath.name=="0-原图"] + tiff_files = [] + for folderpath in target_path: + tiff_files=tiff_files+find_tif_files_pathlib(folderpath) + + tiff_dict={} + for filepath in tiff_files: + rootname=Path(filepath).stem + tiff_dict[rootname]=filepath + + return tiff_dict + + + + + +def read_tifInfo(path): + dataset = gdal.Open(path) # 打开TIF文件 + if dataset is None: + print("无法打开文件,读取文件信息") + return None, None, None + + cols = dataset.RasterXSize # 图像宽度 + rows = dataset.RasterYSize # 图像高度 + bands = dataset.RasterCount + im_proj = dataset.GetProjection() # 获取投影信息 + im_Geotrans = dataset.GetGeoTransform() # 获取仿射变换信息 + # im_data = dataset.ReadAsArray(0, 0, cols, rows) # 读取栅格数据为NumPy数组 + print("行数：", rows) + print("列数：", cols) + print("波段：", bands) + x1=im_Geotrans[0]+im_Geotrans[1]*0 + x2=im_Geotrans[0]+im_Geotrans[1]*cols + + y1=im_Geotrans[3]+im_Geotrans[5]*0 + y2=im_Geotrans[3]+im_Geotrans[5]*rows + + xmin=min(x1,x2) + xmax=max(x1,x2) + ymin=min(y1,y2) + ymax=max(y1,y2) + + geoExtend=[xmin,ymin,xmax,ymax] + del dataset # 关闭数据集 + return im_proj, im_Geotrans,geoExtend + +def getshapefileInfo(shp_path): + """ + 将Shapefile转换为DOTA格式 + :param shp_path: Shapefile文件路径 + """ + + geom_points=[] + + print("shapefile: ",shp_path) + + # 注册所有驱动 + ogr.RegisterAll() + + # 打开Shapefile文件 + driver = ogr.GetDriverByName('ESRI Shapefile') + datasource = driver.Open(shp_path, 0) + if datasource is None: + print("无法打开Shapefile文件") + return + + print("layer count: ",datasource.GetLayerCount()) + for layerid in range(datasource.GetLayerCount()): + print("layer id: ",layerid) + # 获取图层 + layer = datasource.GetLayer(layerid) + layer_defn=layer.GetLayerDefn() + field_count=layer_defn.GetFieldCount() + + print("field_count:", field_count) + for i in range (field_count): + field_defn=layer_defn.GetFieldDefn(i) + field_name=field_defn.GetName() + field_type=field_defn.GetType() + field_type_name=field_defn.GetFieldTypeName(field_type) + print("field_name:", field_name, field_type_name, field_type_name) + + for feature in layer: + geom = feature.GetGeometryRef() + if geom.GetGeometryName() == 'POINT': + x=geom.GetX() + y=geom.GetY() + geom_points.append([x,y]) + return np.array(geom_points) + + +def getTiffsInfo(tiffnames,folderpath): + """ + 获取所有影像的几何信息 + Args: + tiff_paths: tiff列表 + + Returns: + + """ + tiffdict={} + for tiff_name in tiffnames: + if tiff_name.endswith(".tiff"): + tiff_path=os.path.join(folderpath,tiff_name) + im_proj, im_Geotrans, geoExtend=read_tifInfo(tiff_path) + tiffdict[tiff_name]={"geoExtend":geoExtend,"geoTrans":im_Geotrans,"imProj":im_proj} + return tiffdict + + +def getMJSignal(tiffpath,shipPortTree): + im_proj, im_Geotrans, geoExtend = read_tifInfo(tiffpath) # geoExtend : [xmin,ymin,xmax,ymax] + [xmin, ymin, xmax, ymax]=geoExtend + center_x = (xmin + xmax) / 2.0 + center_y = (ymin + ymax) / 2.0 + center_point = [center_x, center_y] + # 2. 计算能够覆盖整个矩形区域的最小半径（中心点到任一角点的最大距离） + radius_to_corner = np.sqrt((xmax - center_x) ** 2 + (ymax - center_y) ** 2) + + MLCFlag=False + JLCFlag=False + ## MLC + if MLCName in shipPortTree and not shipPortTree[MLCName] is None: + # 3. 使用 query_ball_point 查找以中心点为圆心，radius_to_corner 为半径的圆内的所有点的索引 + potential_indices = shipPortTree[MLCName].query_ball_point(center_point, r=radius_to_corner) + + # 4. 获取这些潜在点的实际坐标 + # 假设你的 KDTree 是从 data_points 构建的：MLCTree = KDTree(data_points) + potential_points = shipPortTree[MLCName].data[potential_indices] # 这是所有潜在点的坐标数组 + + # 5. 进行精确的矩形范围过滤 + # 条件判断：x 坐标在 xmin 和 xmax 之间，且 y 坐标在 ymin 和 ymax 之间 + x_in_range = (potential_points[:, 0] >= xmin) & (potential_points[:, 0] <= xmax) + y_in_range = (potential_points[:, 1] >= ymin) & (potential_points[:, 1] <= ymax) + within_rect_indices_mask = x_in_range & y_in_range + + # 6. 获取最终在矩形范围内的点的坐标（在原始 data_points 中的索引是 potential_indices[within_rect_indices_mask]） + final_points = potential_points[within_rect_indices_mask] + + # final_points 就是你要的矩形范围内的点 + # 如果你需要的是这些点在原始数据中的索引，而不是坐标本身： + final_indices = np.array(potential_indices)[within_rect_indices_mask] + + MLCFlag=True + if JLCName in shipPortTree and not shipPortTree[JLCName] is None: + # 3. 使用 query_ball_point 查找以中心点为圆心，radius_to_corner 为半径的圆内的所有点的索引 + potential_indices = shipPortTree[JLCName].query_ball_point(center_point, r=radius_to_corner) + + # 4. 获取这些潜在点的实际坐标 + # 假设你的 KDTree 是从 data_points 构建的：MLCTree = KDTree(data_points) + potential_points = shipPortTree[JLCName].data[potential_indices] # 这是所有潜在点的坐标数组 + + # 5. 进行精确的矩形范围过滤 + # 条件判断：x 坐标在 xmin 和 xmax 之间，且 y 坐标在 ymin 和 ymax 之间 + x_in_range = (potential_points[:, 0] >= xmin) & (potential_points[:, 0] <= xmax) + y_in_range = (potential_points[:, 1] >= ymin) & (potential_points[:, 1] <= ymax) + within_rect_indices_mask = x_in_range & y_in_range + + # 6. 获取最终在矩形范围内的点的坐标（在原始 data_points 中的索引是 potential_indices[within_rect_indices_mask]） + final_points = potential_points[within_rect_indices_mask] + + # final_points 就是你要的矩形范围内的点 + # 如果你需要的是这些点在原始数据中的索引，而不是坐标本身： + final_indices = np.array(potential_indices)[within_rect_indices_mask] + JLCFlag=True + # 处理软件 + return MLCFlag,JLCFlag + + + + +def getTiffInPort(shipPortTree,srcFolderPath_0img,outTiffInfoFilePath): + tiffpaths=find_tif_files_pathlib(srcFolderPath_0img) + tiffLCPort={ + MLCName:[], + JLCName:[], + MJLCName:[], + NOLCName:[] + } + for tiffpath in tiffpaths: + MLCFlag,JLCFlag=getMJSignal(tiffpath,shipPortTree) + + if MLCFlag and JLCFlag: + tiffLCPort[MJLCName].append(tiffpath) + elif MLCFlag: + tiffLCPort[MJLCName].append(tiffpath) + elif JLCFlag: + tiffLCPort[JLCName].append(tiffpath) + else: + tiffLCPort[NOLCName].append(tiffpath) + + # 输出文件 + with open(outTiffInfoFilePath,'w',encoding="utf-8") as f: + for k in tiffLCPort: + for tiffpath in tiffLCPort[k]: + f.write("{} {}\n".format(k,tiffpath)) + + +def SpliteProcess(srcfolderpath,outfolderpath,MLCPath,JLCPath,JMLCPath): + shipPort={ + MLCName:getshapefileInfo(MLCPath), + JLCName:getshapefileInfo(JLCPath), + # "JMLC":getshapefileInfo(JMLCPath), # 舰船不区分 居民一体 + } + + shipPortTree={ + MLCName:KDTree(shipPort[MLCName]), + JLCName:KDTree(shipPort[JLCName]), + # "JMLC":KDTree(shipPort["JMLC"]), + } + srcFolderPath_0img=os.path.join(srcfolderpath,"0-原图") # 0-原图 文件路径 + outTiffInfoFilePath=os.path.join(srcfolderpath,"JMPort.txt") + return True + pass + + +def getParams(): + parser = argparse.ArgumentParser() + parser.add_argument('-s','--srcfolder',type=str,default=r'R:\TYSAR-德清院\TYSAR-条带模式(SM)\港口\20250903-不分类', help='输入shapefile文件') + parser.add_argument('-o', '--outfolder',type=str,default=r'D:\TYSAR-德清院\TYSAR-条带模式(SM)\港口\20250903-不分类', help='输出geojson文件') + parser.add_argument('-m', '--mLC',type=str,help=r'MLC', default=r'D:\TYSAR-德清院\目标点位信息更新\0828目标点位\港口（民船）.shp') + parser.add_argument('-j', '--jLC',type=str,help=r'JLC' ,default=r'D:\TYSAR-德清院\目标点位信息更新\0828目标点位\军港.shp') + parser.add_argument('-jm', '--jmLC',type=str,help=r'MJLC', default=r'D:\TYSAR-德清院\目标点位信息更新\0828目标点位\军民一体港口.shp') + args = parser.parse_args() + return args + +if __name__ == '__main__': + try: + parser = getParams() + srcfolder=parser.srcfolder + outfolder=parser.outfolder + mLCPath=parser.mLC + jLCPath=parser.jLC + jmLCPath=parser.jmLC + print('srcfolder=',srcfolder) + print('outfolder=',outfolder) + print('mLCPath=',mLCPath) + print('jLCPath=',jLCPath) + print('jmLCPath=',jmLCPath) + SpliteProcess(srcfolder,outfolder,mLCPath,jLCPath,jmLCPath) + exit(2) + except Exception as e: + print(e) + exit(3) + + + + +