447 lines
18 KiB
Python
447 lines
18 KiB
Python
from osgeo import ogr
|
||
import os
|
||
import argparse
|
||
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
|
||
|
||
|
||
MLCName="MLC" # M
|
||
JLCName="JLC" # J
|
||
MJLCName="MJLC" # JM 混合
|
||
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,outfolderPath):
|
||
rootname=Path(tiffpath).stem
|
||
portTxtpath=os.path.join(outfolderPath,rootname+".txt")
|
||
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]
|
||
if final_points.shape[0]>0:
|
||
MLCFlag=True
|
||
with open(portTxtpath,"w",encoding="utf-8") as f:
|
||
for i in range(final_points.shape[0]):
|
||
f.write("{}\t\t{},{}\n".format("MLC",final_points[i,0],final_points[i,1]))
|
||
pass
|
||
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]
|
||
if final_points.shape[0]>0:
|
||
JLCFlag=True
|
||
with open(portTxtpath,"a",encoding="utf-8") as f:
|
||
for i in range(final_points.shape[0]):
|
||
f.write("{}\t\t{},{}\n".format("JLC",final_points[i,0],final_points[i,1]))
|
||
pass
|
||
# 处理软件
|
||
return MLCFlag,JLCFlag
|
||
|
||
|
||
def getTiffInPort(shipPortTree,srcFolderPath_0img,outTiffInfoFilePath,outfolderPath):
|
||
tiffpaths=find_tif_files_pathlib(srcFolderPath_0img)
|
||
tiffLCPort={
|
||
MLCName:[],
|
||
JLCName:[],
|
||
MJLCName:[],
|
||
NOLCName:[]
|
||
}
|
||
for tiffpath in tiffpaths:
|
||
MLCFlag,JLCFlag=getMJSignal(tiffpath,shipPortTree,outfolderPath)
|
||
|
||
if MLCFlag and JLCFlag:
|
||
tiffLCPort[MJLCName].append(tiffpath)
|
||
elif MLCFlag:
|
||
tiffLCPort[MLCName].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("{}\t\t{}\n".format(k,tiffpath))
|
||
|
||
|
||
|
||
|
||
def getMJSignal(geoExtend,shipPortTree):
|
||
[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]
|
||
if final_points.shape[0]>0:
|
||
MLCFlag=True
|
||
# with open(portTxtpath,"w",encoding="utf-8") as f:
|
||
# for i in range(final_points.shape[0]):
|
||
# f.write("{}\t\t{},{}\n".format("MLC",final_points[i,0],final_points[i,1]))
|
||
# pass
|
||
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]
|
||
if final_points.shape[0]>0:
|
||
JLCFlag=True
|
||
# with open(portTxtpath,"a",encoding="utf-8") as f:
|
||
# for i in range(final_points.shape[0]):
|
||
# f.write("{}\t\t{},{}\n".format("JLC",final_points[i,0],final_points[i,1]))
|
||
# pass
|
||
# 处理软件
|
||
if MLCFlag and JLCFlag:
|
||
return "mix_airport"
|
||
# tiffLCPort[MJLCName].append(tiffpath)
|
||
elif MLCFlag:
|
||
return "civil_harbor"
|
||
# tiffLCPort[MLCName].append(tiffpath)
|
||
elif JLCFlag:
|
||
return "military_harbor"
|
||
# tiffLCPort[JLCName].append(tiffpath)
|
||
else:
|
||
return "no_harbor"
|
||
# tiffLCPort[NOLCName].append(tiffpath)
|
||
# return MLCFlag,JLCFlag
|
||
|
||
|
||
|
||
def shapefile_to_dota(shp_path, output_path, shipPortTree,difficulty_value=1):
|
||
"""
|
||
将Shapefile转换为DOTA格式
|
||
:param shp_path: Shapefile文件路径
|
||
:param output_path: 输出目录
|
||
:param class_field: 类别字段名
|
||
:param difficulty_value: 难度默认字段
|
||
"""
|
||
|
||
# 注册所有驱动
|
||
ogr.RegisterAll()
|
||
|
||
# 打开Shapefile文件
|
||
driver = ogr.GetDriverByName('ESRI Shapefile')
|
||
datasource = driver.Open(shp_path, 0)
|
||
if datasource is None:
|
||
print("无法打开Shapefile文件")
|
||
return
|
||
|
||
# 获取图层
|
||
layer = datasource.GetLayer()
|
||
|
||
output_file = output_path
|
||
|
||
with open(output_file, 'w',encoding="utf-8") as f:
|
||
# 写入DOTA格式头信息(可选)
|
||
# f.write('imagesource:unknown\n')
|
||
# f.write('gsd:1.0\n')
|
||
|
||
# 遍历所有要素
|
||
for feature in layer:
|
||
# 获取几何对象
|
||
geom = feature.GetGeometryRef()
|
||
if geom is None:
|
||
continue
|
||
# 获取类别和难度
|
||
try:
|
||
class_name = 'unknown'
|
||
except Exception as e:
|
||
class_name="MLC"
|
||
print(e)
|
||
difficulty = difficulty_value
|
||
|
||
# 处理不同类型的几何图形
|
||
if geom.GetGeometryName() == 'POLYGON':
|
||
# 获取多边形外环
|
||
ring = geom.GetGeometryRef(0)
|
||
# 获取所有点
|
||
points = []
|
||
for i in range(ring.GetPointCount()):
|
||
points.append(ring.GetPoint(i))
|
||
|
||
# 确保有足够的点(至少4个)
|
||
if len(points) >= 4:
|
||
# 取前4个点作为DOTA格式的四个角点
|
||
# 注意: DOTA要求按顺序排列(顺时针或逆时针)
|
||
x1, y1 = points[0][0], points[0][1]
|
||
x2, y2 = points[1][0], points[1][1]
|
||
x3, y3 = points[2][0], points[2][1]
|
||
x4, y4 = points[3][0], points[3][1]
|
||
|
||
xmin = min(x1, x2,x3,x4)
|
||
xmax = max(x1, x2,x3,x4)
|
||
ymin = min(y1, y2,y3,y4)
|
||
ymax = max(y1, y2,y3,y4)
|
||
# [xmin, ymin, xmax, ymax] = geoExtend
|
||
geoExtend = [xmin, ymin, xmax, ymax]
|
||
class_name=getMJSignal(geoExtend, shipPortTree)
|
||
# 写入DOTA格式行
|
||
line = f"{x1} {y1} {x2} {y2} {x3} {y3} {x4} {y4} {class_name} {difficulty}\n"
|
||
f.write(line)
|
||
|
||
# 释放资源
|
||
datasource.Destroy()
|
||
print("转换完毕")
|
||
|
||
|
||
def PortShapeProces(shp_path, output_path, MLCPath,JLCPath,JMLCPath,difficulty_value=1):
|
||
shipPort={
|
||
MLCName:getshapefileInfo(MLCPath),
|
||
JLCName:getshapefileInfo(JLCPath),
|
||
MJLCName:getshapefileInfo(JMLCPath), # 舰船不区分 居民一体
|
||
}
|
||
|
||
shipPortTree={
|
||
MLCName:KDTree(shipPort[MLCName]),
|
||
JLCName:KDTree(shipPort[JLCName]),
|
||
MJLCName:KDTree(shipPort[MJLCName]),
|
||
}
|
||
shapefile_to_dota(shp_path, output_path,shipPortTree, difficulty)
|
||
|
||
|
||
|
||
def getParams():
|
||
parser = argparse.ArgumentParser()
|
||
parser.add_argument('-i','--infile',type=str,default=r'D:\港口\港口\Geo_bc2-sm-org-vv-20231016t135315-008424-0020e8-01.military_harbor.shp', help='输入shapefile文件')
|
||
parser.add_argument('-o', '--outfile',type=str,default=r'D:\港口\港口dota\Geo_bc2-sm-org-vv-20231016t135315-008424-0020e8-01.military_harbor.txt', 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'JMLC', default=r'D:\TYSAR-德清院\目标点位信息更新\0828目标点位\君民一体港口.shp')
|
||
parser.add_argument('-d', '--difficulty',type=int,default=1, help='输出geojson文件')
|
||
args = parser.parse_args()
|
||
return args
|
||
|
||
if __name__ == '__main__':
|
||
try:
|
||
parser = getParams()
|
||
inFilePath=parser.infile
|
||
outpath=parser.outfile
|
||
mLCPath=parser.mLC
|
||
jLCPath=parser.jLC
|
||
jmLCPath=parser.jmlc
|
||
difficulty=parser.difficulty
|
||
print('infile=',inFilePath)
|
||
print('outfile=',outpath)
|
||
print('mLCPath=',mLCPath)
|
||
print('jLCPath=',jLCPath)
|
||
print('jmLCPath=',jmLCPath)
|
||
print('difficulty=',difficulty)
|
||
PortShapeProces(inFilePath, outpath, mLCPath, jLCPath, jmLCPath, difficulty_value=difficulty)
|
||
exit(2)
|
||
except Exception as e:
|
||
print(e)
|
||
exit(3)
|