288 lines
10 KiB
Python
288 lines
10 KiB
Python
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)
|
||
|
||
|
||
|
||
|
||
|