检测机场中的机场跑道

输入参数（飞机的起降需要的长，宽，机场的四点经纬度（共8个），爆炸点的经纬度（共两个），爆炸点的范围，爆炸点的角度）

可输入多个爆炸点进行检测

测试数据：450 30 144.9127545108085 13.5793801633106 144.9130426945087 13.57889946157555  144.9454464308249 13.59314957370197 144.9452087239189 13.59362516033664 144.9432503601635 13.59240186299597 5 20 144.9409857674664 13.59130151250747 5 30

import os
import sys
import cv2
import numpy as np
from haversine import haversine, Unit
import math


def detectRunway(Location_array,boomlist,brng,taskWidth,taskLength):

    airWidth = haversine((Location_array[1],Location_array[0]),(Location_array[3],Location_array[2]), unit=Unit.METERS)#39.995304, 116.308264
    airLength = haversine((Location_array[3],Location_array[2]),(Location_array[5],Location_array[4]), unit=Unit.METERS)
    if airWidth > airLength:
        airWidth,airLength = airLength,airWidth
        tmp1 = Location_array[0]
        tmp2 = Location_array[1]
        Location_array[0] = Location_array[6]
        Location_array[1] = Location_array[7]
        Location_array[6] = Location_array[4]
        Location_array[7] = Location_array[5]
        Location_array[4] = Location_array[2]
        Location_array[5] = Location_array[3]
        Location_array[2] = tmp1
        Location_array[3] = tmp2
        brng = geoDegree(Location_array[0],Location_array[1],Location_array[6],Location_array[7])# degree with north
        
    for i in range (0,boomlist.shape[0]):
        pointDis = haversine((Location_array[3],Location_array[2]),(boomlist[i,1],boomlist[i,0]), unit=Unit.METERS)
        pointAngle = float(boomlist[i,3])   
        pointDis_x = math.cos(math.radians(pointAngle))*pointDis
        pointDis_y = math.sin(math.radians(pointAngle))*pointDis
        boomlist_pixel[i,0] = int(pointDis_x)
        boomlist_pixel[i,1] = int(pointDis_y)
        boomlist_pixel[i,2] = int(boomlist[i,2])
        boomlist_pixel[i,3] = float(pointAngle)


        print("爆炸点的坐标是"+str(boomlist_pixel[i,0])+","+str(boomlist_pixel[i,1])+",爆炸点的的角度是："+str(boomlist_pixel[i,3]) +",爆炸点的范围是"+str(boomlist_pixel[i,2])+",爆炸点的相对距离："+str(pointDis))
    taskWidth = int(taskWidth)
    taskLength = int(taskLength)
    takeoff_pixel,flag_bool = runWayRange(boomlist_pixel,airWidth,airLength,taskWidth,taskLength)

    #################################Calculate airport longitude and latitude############################################
    out_str = ''
    if flag_bool == True:

        from geographiclib.geodesic import Geodesic
        geod = Geodesic.WGS84  # define the WGS84 ellipsoid
        geod = Geodesic(6378388, 1/297.0) # altanatively custom the ellipsoid
               
        for i in range(1,(int(len(takeoff_pixel)/8))+1):
            if i >1 :
                weishu = 1
            else :
                weishu = 0

            pixel_0 = takeoff_pixel[0+(i-1)*8]
            pixel_1 = takeoff_pixel[1+(i-1)*8]
            pixel_2 = takeoff_pixel[2+(i-1)*8]
            pixel_3 = takeoff_pixel[3+(i-1)*8]
            pixel_4 = takeoff_pixel[4+(i-1)*8]
            pixel_5 = takeoff_pixel[5+(i-1)*8] 
            pixel_6 = takeoff_pixel[6+(i-1)*8]
            pixel_7 = takeoff_pixel[7+(i-1)*8]

            dis_1 = math.sqrt(pixel_0**2+ pixel_1**2)
            dis_1_angle = cal_ang([pixel_0+1,pixel_1+1],[0,0],[pixel_0+1,0]) + brng
            g1 = geod.Direct(Location_array[1],Location_array[0],dis_1_angle,dis_1)
            

            dis_2 = math.sqrt(pixel_2**2+ pixel_3**2)
            dis_2_angle = cal_ang([pixel_2+1,pixel_3+1],[0,0],[pixel_2+1,0]) + brng
            g2 = geod.Direct(Location_array[1],Location_array[0],dis_2_angle,dis_2)

            dis_3 = math.sqrt(pixel_4**2+ pixel_5**2)
            dis_3_angle = cal_ang([pixel_4+1,pixel_5+1],[0,0],[pixel_4+1,0]) + brng
            g3 = geod.Direct(Location_array[1],Location_array[0],dis_3_angle,dis_3)

            dis_4 = math.sqrt(pixel_6**2+ pixel_7**2)
            dis_4_angle = cal_ang([pixel_6+1,pixel_7+1],[0,0],[pixel_6+1,0]) + brng
            g4 = geod.Direct(Location_array[1],Location_array[0],dis_4_angle,dis_4)
            
            outstr_tmp = '[Output]_'+str(g1['lon2'])+'_'+str(g1['lat2'])+'_'+str(g2['lon2'])+'_'+str(g2['lat2'])+'_'+str(g3['lon2'])+'_'+str(g3['lat2'])+'_'+str(g4['lon2'])+'_'+str(g4['lat2'])+'_'
            out_str = out_str +outstr_tmp
            #print('[Output]_'+str(g1['lon2'])+'_'+str(g1['lat2'])+'_'+str(g2['lon2'])+'_'+str(g2['lat2'])+'_'+str(g3['lon2'])+'_'+str(g3['lat2'])+'_'+str(g4['lon2'])+'_'+str(g4['lat2'])+'_')
        #print('[Output]_'+str(final_lon_1)+'_'+str(final_lat_1)+'_'+str(final_lon_2)+'_'+str(final_lat_2)+'_'+str(final_lon_3)+'_'+str(final_lat_3)+'_'+str(final_lon_4)+'_'+str(final_lat_4))
    if out_str == '':
        pass
    else:
        print(out_str)
    return

def runWayRange(boomlist_pixel,airWidth,airLength,taskWidth,taskLength):
    
    baseMap = np.zeros((int(airWidth),int(airLength),3))
    for i in range (0,boomlist.shape[0]):
        #print('爆炸点的圆心是:'+str(boomlist_pixel[i,0])+','+str(boomlist_pixel[i,1])+', 爆炸点的半径是'+str(boomlist_pixel[i,2]))
        cv2.circle(baseMap,(int(boomlist_pixel[i,0]),int(boomlist_pixel[i,1])),int(boomlist_pixel[i,2]),(255,255,255),-1)
    #print(type(baseMap))
    baseMap = np.flip(baseMap,axis=0)
#####################################爆炸点测试数据##########################
    #cv2.circle(baseMap,(10,10),10,(255,255,255),-1)
    #cv2.circle(baseMap,(1922,30),15,(255,255,255),-1)
    #cv2.circle(baseMap,(800,30),15,(255,255,255),-1)
#############################################################################
    cv2.imwrite('basemap.jpg',baseMap)
    #cv2.imshow('image',baseMap)
    #cv2.waitKey(0)
    baseMap = baseMap.astype(int)
    flag_bool = False
    takeoff_pixel = []
    for i in range(0,len(baseMap[0])-taskWidth):#x
        for j in range(0,len(baseMap)-taskLength):#y
            #baseMap的shape参数是y，x，z
            '''
450 30 144.9127545108085 13.5793801633106 144.9130426945087 13.57889946157555  144.9454464308249 13.59314957370197 144.9452087239189 13.59362516033664 144.9432503601635 13.59240186299597 5 144.9409857674664 13.59130151250747 5
400 50 115.0856523072000 29.2972547380000 115.0960060420000 29.2971791631000 115.0958171052000 29.3006744992000 115.0857845641000 29.3009012237000 115.0879195484000 29.2976137184000 26.3579722676000 4 115.0943245048000 29.3004099873000 49.0883069438000  8.9 115.0905004495000 29.3003804254000 27.8902681093000 15.4
            '''
            tmp_zuoshang_y = j
            tmp_zuoshang_x = i
            #print(str(i)+','+str(j))
#            if baseMap[j:taskLength+j+1,i:taskWidth+i+1].max() < 254.0:
            if baseMap[j:taskLength+j+1,i:taskWidth+i+1].max() < 254.0:
                tmp_youxia_x = taskWidth+i+1
                tmp_youxia_y = taskLength+j+1
                #print("矩阵范围是：("+str(i)+","+str(j)+")到("+str(tmp_youxia_x)+","+str(tmp_youxia_y)+");"+"图像边界是：("+str(len(baseMap[0]))+","+str(len(baseMap))+")")
                takeoff_pixel_tmp = [i,j,i,taskLength+j+1,i+taskWidth+1,taskLength+j+1,i+taskWidth+1,j]
                ######################################################长度搜索########################################################
                for k in range(tmp_youxia_x+1,len(baseMap[0])):
                    if baseMap[tmp_zuoshang_y:tmp_youxia_y,tmp_youxia_x:k].max() < 254.0:
                        tmp_youxia_x = k
                        takeoff_pixel_chang = [tmp_zuoshang_x,tmp_zuoshang_y,tmp_zuoshang_x,tmp_youxia_y,tmp_youxia_x,tmp_youxia_y,tmp_youxia_x,tmp_zuoshang_y]#逆时针坐标顺序  
                        takeoff_pixel_tmp = takeoff_pixel_chang
                    else:
                        break
                ######################################################宽度搜索########################################################
                for l in range(tmp_youxia_y+1,len(baseMap)):
                    if baseMap[tmp_youxia_y:l,tmp_zuoshang_x:tmp_youxia_x].max() < 254.0:
                        tmp_youxia_y = l
                        takeoff_pixel_kuan = [tmp_zuoshang_x,tmp_zuoshang_y,tmp_zuoshang_x,tmp_youxia_y,tmp_youxia_x,tmp_youxia_y,tmp_youxia_x,tmp_zuoshang_y]#逆时针坐标顺序
                        takeoff_pixel_tmp = takeoff_pixel_kuan
                    else:
                        break
                ######################################################搜索结束########################################################                    
                flag_bool =True
                #print(type(takeoff_pixel))
                cv2.rectangle(baseMap,(int(tmp_zuoshang_x),int(tmp_zuoshang_y)),(int(tmp_youxia_x),int(tmp_youxia_y)),(255,255,255),thickness=-1)
                takeoff_pixel.extend(takeoff_pixel_tmp)
                #cv2.imshow('image',baseMap)
                #cv2.waitKey(0)
                #break#搜索出单个机场范围
                #print("next")
    if flag_bool == True:
        print("[True]The Minimum takeoff and landing window coordinate is exist")
    else:
        print("[False]_The airport runway does not have a minimum takeoff and landing window")
    #cv2.imshow('image',baseMap)
    cv2.imwrite('testsave.jpg',baseMap)
    #cv2.waitKey(0)
    return takeoff_pixel,flag_bool

def cal_ang(point_1, point_2, point_3):
    """
    根据三点坐标计算夹角
    :param point_1: 点1坐标
    :param point_2: 点2坐标
    :param point_3: 点3坐标
    :return: 返回任意角的夹角值，这里只是返回点2的夹角
    """
    a = math.sqrt(
        (point_2[0] - point_3[0]) * (point_2[0] - point_3[0]) + (point_2[1] - point_3[1]) * (point_2[1] - point_3[1]))
    b = math.sqrt(
        (point_1[0] - point_3[0]) * (point_1[0] - point_3[0]) + (point_1[1] - point_3[1]) * (point_1[1] - point_3[1]))
    c = math.sqrt(
        (point_1[0] - point_2[0]) * (point_1[0] - point_2[0]) + (point_1[1] - point_2[1]) * (point_1[1] - point_2[1]))
    A = math.degrees(math.acos((a * a - b * b - c * c) / (-2 * b * c)))
    B = math.degrees(math.acos((b * b - a * a - c * c) / (-2 * a * c)))
    C = math.degrees(math.acos((c * c - a * a - b * b) / (-2 * a * b)))
    return B




def geoDegree(lng1,lat1,lng2,lat2):

    '''
    公式计算两点间方位角
    方位角：是与正北方向、顺时针之间的夹角
    '''
    lng1,lat1,lng2,lat2=map(math.radians,[float(lng1),float(lat1),float(lng2),float(lat2)]) 
    dlon=lng2-lng1
    y=math.sin(dlon)*math.cos(lat2)
    x=math.cos(lat1)*math.sin(lat2)-math.sin(lat1)*math.cos(lat2)*math.cos(dlon)
    brng=math.degrees(math.atan2(y,x))
    brng=(brng+360)%360
    return brng

if __name__ == '__main__':
    '''
450 30 144.9127545108085 13.5793801633106 144.9130426945087 13.57889946157555  144.9454464308249 13.59314957370197 144.9452087239189 13.59362516033664 144.9432503601635 13.59240186299597 5 20 144.9409857674664 13.59130151250747 5 30
    '''
    args = sys.argv
    args = input()
    string_list  = args.split()

    count_boom = ((len(string_list)-10)/4)
    if (count_boom % 1 == 0) == True:
        boomlist  = np.zeros((int(count_boom),4))
        boomlist_pixel  = np.zeros((int(count_boom),4))#[[0]*3]*int(count_boom)
        print("[Info]_Enter a total of {} explosion point".format(int(count_boom)))
        taskWidth = string_list[0]
        taskLength = string_list[1]
        
        lon1 = float(string_list[0+2])
        lat1 = float(string_list[1+2])
        lon2 = float(string_list[2+2])
        lat2 = float(string_list[3+2])
        lon3 = float(string_list[4+2])
        lat3 = float(string_list[5+2])
        lon4 = float(string_list[6+2])
        lat4 = float(string_list[7+2])
        brng = geoDegree(lon1,lat1,lon4,lat4)# degree with north
        x_label = 0 
        for i in range(10,len(string_list),4):
            for j in range(0,int(count_boom)):
                boomlist[x_label,0] = float(string_list[i])
                boomlist[x_label,1] = float(string_list[i+1])
                boomlist[x_label,2] = float(string_list[i+2])
                boomlist[x_label,3] = float(string_list[i+3])
            i + 4
            x_label += 1
        Location_array = [lon1,lat1,lon2,lat2,lon3,lat3,lon4,lat4]
        detectRunway(Location_array,boomlist,brng,taskWidth,taskLength)
        print('[Info]_Airplane minimum take-off and landing detection algorithm over')

    else:
        print("[Error]_The input of the explosion point parameters does not meet the requirements, Please check")