Distance and Bearing: Difference between revisions

Content added Content deleted

Inline

Revision as of 08:44, 14 November 2022

It is very important in aviation to have knowledge of the nearby airports at any time in flight.

Task

Determine the distance and bearing from an Airplane to the 20 nearest Airports whenever requested. Use the non-commercial data from openflights.org airports.dat as reference.

A request comes from an airplane at position ( latitude, longitude ): ( 51.514669, 2.198581 ).

Your report should contain the following information from table airports.dat (column shown in brackets):

Name(2), Country(4), ICAO(6), Distance and Bearing calculated from Latitude(7) and Longitude(8).

Distance is measured in nautical miles (NM). Resolution is 0.1 NM.

Bearing is measured in degrees (°). 0° = 360° = north then clockwise 90° = east, 180° = south, 270° = west. Resolution is 1°.

See

openflights.org/data: Airport, airline and route data
Movable Type Scripts: Calculate distance, bearing and more between Latitude/Longitude points

J

Implementation:

DAT=: {{ require'web/gethttp csv'
  if.0=#$fread F=.'airports.dat' do.
    (gethttp y) fwrite F
  end.
  readcsv F
}}'https://raw.githubusercontent.com/jpatokal/openflights/master/data/airports.dat'
LATLONG=: _".&>6 7{"1 DAT

fmtm=: {{
 9!:7]9 1 1#'  -'
 r=._1 _1}.1 1}.":boxopen y
 ~.r[9!:7]9 1 1#'+|-'
}}
AIRPORTS=: 1 3 5{"1 DAT

deg=: %&180p_1     NB. convert from degrees
R=:  21600 % 2p1   NB. radius of spherical  earth in nautical miles

dist=: {{R*2*_1 o.1<.%:(*:1 o.0.5*x-y)p.x*&(2 o.{.)y}}&deg
bear=: {{
  re=. (2 o.{.y)*1 o.dlon=.y-&{:x
  2p1|1r2p1-{:*.re j. ((*|.)/+.^j.x,&{.y) p.-2 o.dlon
}}&.deg

NB. round y to nearest whole number under multiplying by x
pkg=: {{ <"0 (<.0.5+x*y)%x }}

LABELS=: ;:'Airport Country ICAO Distance Bearing'

nearest=: {{
  ndx=. 20{./: d=. y dist"1 LATLONG     NB. 20 nearest airports
  fmtm LABELS,(ndx{AIRPORTS),"1 (10 pkg ndx{d),.1 pkg y bear"1 ndx{LATLONG
}}

For this task, we assume that the earth is a sphere rather than an oblate ellipsoid. Also, we report the initial bearing (at the requested location). For the task example, these issues are irrelevant (the distances involved are too small for those details to be significant), but they might become significant if code from this implementation was used in other contexts.

Task example:

   nearest 51.514669 2.198581
Airport                             Country        ICAO Distance Bearing 
----------------------------------- -------------- ---- -------- ------- 
Koksijde Air Base                   Belgium        EBFN 30.6     146     
Ostend-Bruges International Airport Belgium        EBOS 31.3     127     
Kent International Airport          United Kingdom EGMH 33.5     252     
Calais-Dunkerque Airport            France         LFAC 34.4     196     
Westkapelle heliport                Belgium        EBKW 42.5     105     
Lympne Airport                      United Kingdom EGMK 51.6     240     
Ursel Air Base                      Belgium        EBUL 52.8     114     
Southend Airport                    United Kingdom EGMC 56.2     274     
Merville-Calonne Airport            France         LFQT 56.3     163     
Wevelgem Airport                    Belgium        EBKT 56.4     137     
Midden-Zeeland Airport              Netherlands    EHMZ 57.2     90      
Lydd Airport                        United Kingdom EGMD 58       235     
RAF Wattisham                       United Kingdom EGUW 58.9     309     
Beccles Airport                     United Kingdom EGSM 59.3     339     
Lille/Marcq-en-Baroeul Airport      France         LFQO 59.6     146     
Lashenden (Headcorn) Airfield       United Kingdom EGKH 62.2     250     
Le Touquet-Côte d'Opale Airport     France         LFAT 63.7     200    
Rochester Airport                   United Kingdom EGTO 64.2     262     
Lille-Lesquin Airport               France         LFQQ 66.2     149     
Thurrock Airfield                   United Kingdom EGMT 68.4     272

jq

Adapted from Wren

Works with: jq

Works with gojq, the Go implementation of jq.

jq does not natively support the complexities of CSV, so to focus on the task at hand, this entry assumes that the data file is in TSV format. However, a jq program for converting CSV to JSON arrays may be found here on RC at Convert_CSV_records_to_TSV; it has been tested against the airlines.dat file.

def pi: 1|atan * 4;

# Read a TSV file one line at a time,
# and emit a stream of arrays, one array per line:
def readTSV: [inputs | sub("  *$";"") / "\t" ];

def airports:
  [readTSV
   |    {airportID: (.[0] | tonumber),
         name      : .[1],
         city      : .[2],
         country   : .[3],
         iata      : .[4],
         icao      : .[5],
         latitude  : (.[6] | tonumber),
         longitude : (.[7] | tonumber),
         altitude  : (.[8] | tonumber),
         timezone  : .[9],
         dst       : .[10],
         tzOlson   : .[11],
         type      : .[12],
         source    : .[13]} ];
    
def calculateDistance(lat1; lon1; lat2; lon2; units):
  if (lat1 == lat2 and lon1 == lon2) then 0
  else pi as $pi
  | ($pi * lat1 / 180) as $radlat1
  | ($pi * lat2 / 180) as $radlat2
  | (lon1 - lon2) as $theta
  | ($pi * $theta / 180) as $radtheta
  | {}
  | .dist = ($radlat1|sin) * ($radlat2|sin) + 
            ($radlat1|cos) * ($radlat2|cos) * ($radtheta|cos)
  | if .dist > 1 then .dist = 1 else . end
  | .dist |= (acos * 180 / $pi) * 60 * 1.1515576  # distance in statute miles
  | if   (units == "K") then .dist *= 1.609344          # distance in kilometers
    elif (units == "N") then .dist *= 0.868976          # distance in nautical miles
    else .
    end
  | .dist
  end ;

def calculateBearing(lat1; lon1; lat2; lon2):
  if (lat1 == lat2 and lon1 == lon2) then 0
  else pi as $pi
  | ($pi * lat1 / 180) as $radlat1
  | ($pi * lat2 / 180) as $radlat2
  | ($pi * (lon2 - lon1) / 180) as $raddlon
  | ( ($raddlon|sin) * ($radlat2|cos) ) as $y
  | (($radlat1|cos) * ($radlat2|sin) -
     ($radlat1|sin) * ($radlat2|cos) * ($raddlon|cos)) as $x
  | ((atan2($y;$x) * 180 / $pi) + 360) % 360
  end;

# request from airplane at position $X; $Y
def request($X; $Y):
  airports as $airports
  | reduce $airports[] as $a ([];
        ((calculateDistance($X; $Y; $a.latitude; $a.longitude; "N") * 10 | round)
         / 10) as $dist
      | (calculateBearing($X; $Y; $a.latitude; $a.longitude)|round ) as $bear
      | . + [ [$a.name, $a.country, $a.icao, $dist, $bear] ] ) ;

def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;

def fmt:
  def l($n): lpad($n);
  "\(.[0] |l(36)) | \(.[1]|l(14)) | \(.[2]|l(4)) |\(.[3]|l(15)) |\(.[4]|l(10))";

# Example:
def display($n):
  "Closest \($n) destinations relative to position \(.[0]), \(.[1]):",
  "",
  "                Name                 |    Country     | ICAO | Distance in NM | Bearing in ° ",
  "-------------------------------------+----------------+------+----------------+--------------",
  (request(.[0]; .[1])
   | sort_by( .[3] )
   | .[0:$n][] | fmt ) ;


[51.514669, 2.198581] | display(20)

Invocation: jq -nrRf distance-and-bearing.jq airports.tsv

Output:

Closest 20 destinations relative to position 51.514669, 2.198581:

                Name                 |    Country     | ICAO | Distance in NM | Bearing in ° 
-------------------------------------+----------------+------+----------------+--------------
                   Koksijde Air Base |        Belgium | EBFN |           30.7 |       146
 Ostend-Bruges International Airport |        Belgium | EBOS |           31.3 |       127
          Kent International Airport | United Kingdom | EGMH |           33.5 |       252
            Calais-Dunkerque Airport |         France | LFAC |           34.4 |       195
                Westkapelle heliport |        Belgium | EBKW |           42.6 |       105
                      Lympne Airport | United Kingdom | EGMK |           51.6 |       240
                      Ursel Air Base |        Belgium | EBUL |           52.8 |       114
                    Southend Airport | United Kingdom | EGMC |           56.2 |       274
            Merville-Calonne Airport |         France | LFQT |           56.4 |       162
                    Wevelgem Airport |        Belgium | EBKT |           56.5 |       137
              Midden-Zeeland Airport |    Netherlands | EHMZ |           57.3 |        89
                        Lydd Airport | United Kingdom | EGMD |             58 |       235
                       RAF Wattisham | United Kingdom | EGUW |             59 |       309
                     Beccles Airport | United Kingdom | EGSM |           59.3 |       339
      Lille/Marcq-en-Baroeul Airport |         France | LFQO |           59.7 |       146
       Lashenden (Headcorn) Airfield | United Kingdom | EGKH |           62.2 |       250
     Le Touquet-Côte d'Opale Airport |         France | LFAT |           63.7 |       200
                   Rochester Airport | United Kingdom | EGTO |           64.2 |       261
               Lille-Lesquin Airport |         France | LFQQ |           66.2 |       149
                   Thurrock Airfield | United Kingdom | EGMT |           68.4 |       271

Julia

using DataFrames
using CSV

const EARTH_RADIUS_KM = 6372.8
const TASK_CONVERT_NM = 0.0094174
const AIRPORT_DATA_FILE = "airports.dat.txt"

const QUERY = (latitude =  51.514669, longitude = 2.198581)

"""
Given two latitude, longitude pairs in degrees for two points on the Earth,
get distance (nautical miles) and initial direction of travel (degrees)
for travel from lat1, lon1 to lat2, lon2
"""
function haversine(lat1, lon1, lat2, lon2)
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = sind(dlat / 2)^2 + cosd(lat1) * cosd(lat2) * sind(dlon / 2)^2
    c = 2.0 * asind(sqrt(a))
    theta = atand(sind(dlon) * cosd(lat2),
        cosd(lat1) * sind(lat2) - sind(lat1) * cosd(lat2) * cosd(dlon))
    theta = (theta + 360) % 360
    return EARTH_RADIUS_KM * c * TASK_CONVERT_NM, theta
end

""" Given latitude and longitude, find `wanted` closest airports in database file csv. """
function find_nearest_airports(latitude, longitude, wanted = 20, csv = AIRPORT_DATA_FILE)
    airports = DataFrame(CSV.File(csv))[:, [2, 4, 6, 7, 8]]
    airports[!, "d"] .= 0.0
    airports[!, "b"] .= 0
    rename!(airports, [:Name, :Country, :ICAO, :Latitude_deg, :Longitude_deg, :Distance, :Bearing])
    for r in eachrow(airports)
        distance, bearing = haversine(latitude, longitude, r.Latitude_deg, r.Longitude_deg)
        r.Distance, r.Bearing = round(distance, digits = 1), Int(round(bearing))
    end
    sort!(airports, :Distance)
    return airports[1:wanted, [:Name, :Country, :ICAO, :Distance, :Bearing]]
end

println(find_nearest_airports(QUERY.latitude, QUERY.longitude))

Output:

20×5 DataFrame
 Row │ Name                               Country         ICAO     Distance  Bearing 
     │ String                             String          String7  Float64   Int64   
─────┼───────────────────────────────────────────────────────────────────────────────
   1 │ Koksijde Air Base                  Belgium         EBFN         30.6      146 
   2 │ Ostend-Bruges International Airp…  Belgium         EBOS         31.3      127 
   3 │ Kent International Airport         United Kingdom  EGMH         33.5      252 
   4 │ Calais-Dunkerque Airport           France          LFAC         34.4      196 
   5 │ Westkapelle heliport               Belgium         EBKW         42.6      105 
   6 │ Lympne Airport                     United Kingdom  EGMK         51.6      240 
   7 │ Ursel Air Base                     Belgium         EBUL         52.8      114 
   8 │ Southend Airport                   United Kingdom  EGMC         56.2      274 
   9 │ Merville-Calonne Airport           France          LFQT         56.3      163 
  10 │ Wevelgem Airport                   Belgium         EBKT         56.4      137 
  11 │ Midden-Zeeland Airport             Netherlands     EHMZ         57.2       90 
  12 │ Lydd Airport                       United Kingdom  EGMD         58.0      235 
  13 │ RAF Wattisham                      United Kingdom  EGUW         59.0      309 
  14 │ Beccles Airport                    United Kingdom  EGSM         59.3      339 
  15 │ Lille/Marcq-en-Baroeul Airport     France          LFQO         59.7      146 
  16 │ Lashenden (Headcorn) Airfield      United Kingdom  EGKH         62.2      250 
  17 │ Le Touquet-Côte d'Opale Airport    France          LFAT         63.7      200
  18 │ Rochester Airport                  United Kingdom  EGTO         64.2      262
  19 │ Lille-Lesquin Airport              France          LFQQ         66.2      149
  20 │ Thurrock Airfield                  United Kingdom  EGMT         68.4      272

PARI/GP

/* copy CSV aiports.dat from URL */
csv=externstr("curl 'https://raw.githubusercontent.com/jpatokal/openflights/master/data/airports.dat'");

/* We want an easy convert from CSV to matrix. */
/* Remove all those backslashes in airports.dat */
{
    for(i = 1 , #csv, tmp = Vecsmall(csv[i]);
        for(j = 1 , #tmp, if(tmp[j]==92, tmp[j]=78));
        csv[i] = strchr(tmp)
    );
};
 
{  
convert(csv,mat)=
    write1(mat,"[");
    for(i=1, #csv-1,
        write1(mat,csv[i]); write1(mat,";")
    );
    write1(mat,csv[#csv]); write(mat,"]")
};
 
convert(csv,airports);
M = read(airports);

{
distance(lat1, lon1, lat2, lon2, unit = "mi") =
    my(dist = 0, rlat1, rlat2, rtheta);
    if(lat1!=lat2 || lon1!=lon2,
        rlat1 = lat1*Pi/180; rlat2 = lat2*Pi/180;
        rtheta = (lon1-lon2)*Pi/180;
        dist = sin(rlat1)*sin(rlat2) + cos(rlat1)*cos(rlat2)*cos(rtheta);
        dist = acos(min(dist,1))*180/Pi*60*1.1515576;
        if(unit=="Km", dist*= 1.609344, unit=="NM", dist*= 0.868976)
    );dist
};

{
bearing(lat1, lon1, lat2, lon2)=
    my(bear, rlat1, rlat2, rdlon, x, y);
    if (lat1!=lat2 || lon1!=lon2,
        rlat1 = lat1*Pi/180; rlat2 = lat2*Pi/180;
        rdlon = (lon2-lon1)*Pi/180;
        y = sin(rdlon)*cos(rlat2);
        x = cos(rlat1)*sin(rlat2) - sin(rlat1)*cos(rlat2)*cos(rdlon);
        /* Pari/GP has no built-in atan2() function. */
        if( x>0, bear=atan(y/x), x<0 && y>0, bear=atan(y/x)+Pi,
            x<0 && y==0, bear=Pi, x<0 && y<0, bear=atan(y/x)-Pi,
            x==0 && y>0, bear=Pi/2, x==0 && y<0, bear=-Pi/2);
        bear = (bear*180/Pi+360)%360;
    );bear
};

{
request(lat, lon, r = 20 , unit = "NM")=
    my(D = vector(#M[,1]), R = matrix(r+1, 5));
    R[1,1] = "NAME"; R[1,2] = "COUNTRY"; R[1,3] = "ICAO"; R[1,4] = "DISTANCE"; R[1,5] = "BEARING";
    printf("%-38s\t%-16s\t%-4s\t%s\t%s\n", R[1,1], R[1,2], R[1,3], R[1,4], R[1,5]);  
    for(i = 1, #M[,1],
        D[i] = [distance(lat, lon, M[i,7], M[i,8], unit) , i];
    );
    D=vecsort(D);
    for(i = 1, r,
        R[i+1,1] = M[D[i][2],2]; R[i+1,2] = M[D[i][2],4]; R[i+1,3] = M[D[i][2],6];
        R[i+1,4] = D[i][1]; R[i+1,5] = bearing(lat, lon, M[D[i][2],7], M[D[i][2],8]);
        printf("%-38s\t%-16s\t%-4s\t%8.1f\t%7.0f\n", R[i+1,1], R[i+1,2], R[i+1,3], R[i+1,4], R[i+1,5]);    
    );
} 

request(51.514669, 2.198581)

Output:

NAME                                  	COUNTRY         	ICAO	DISTANCE	BEARING
Koksijde Air Base                     	Belgium         	EBFN	    30.7	    146
Ostend-Bruges International Airport   	Belgium         	EBOS	    31.3	    127
Kent International Airport            	United Kingdom  	EGMH	    33.5	    252
Calais-Dunkerque Airport              	France          	LFAC	    34.4	    196
Westkapelle heliport                  	Belgium         	EBKW	    42.6	    105
Lympne Airport                        	United Kingdom  	EGMK	    51.6	    240
Ursel Air Base                        	Belgium         	EBUL	    52.8	    114
Southend Airport                      	United Kingdom  	EGMC	    56.2	    274
Merville-Calonne Airport              	France          	LFQT	    56.4	    163
Wevelgem Airport                      	Belgium         	EBKT	    56.5	    137
Midden-Zeeland Airport                	Netherlands     	EHMZ	    57.3	     90
Lydd Airport                          	United Kingdom  	EGMD	    58.0	    235
RAF Wattisham                         	United Kingdom  	EGUW	    59.0	    309
Beccles Airport                       	United Kingdom  	EGSM	    59.3	    339
Lille/Marcq-en-Baroeul Airport        	France          	LFQO	    59.7	    146
Lashenden (Headcorn) Airfield         	United Kingdom  	EGKH	    62.2	    250
Le Touquet-Côte d'Opale Airport      	France          	LFAT	    63.7	    200
Rochester Airport                     	United Kingdom  	EGTO	    64.2	    262
Lille-Lesquin Airport                 	France          	LFQQ	    66.2	    149
Thurrock Airfield                     	United Kingdom  	EGMT	    68.4	    272
cpu time = 183 ms, real time = 183 ms.

Phix

without js -- file i/o
enum Airport_ID, Name, City, Country, IATA, ICAO, Latitude, Longitude, Altitude, Timezone, DST, Tz_Olson, Type, Source
function tabulate(sequence s)
    s = split(trim(s),"\n")
    for i,ri in s do
        ri = split(ri,',')
        if length(ri)!=14 then
            ri[2] &= ","&ri[3]
            ri[3..3] = {}
        end if
        for j in {Airport_ID, Latitude, Longitude, Altitude} do
            ri[j] = to_number(ri[j])
        end for
        for j in {Name,City,Country,ICAO} do
            ri[j] = trim(ri[j],`"`)
        end for
        s[i] = ri
    end for
    return s    
end function

constant airports = tabulate(get_text("airports.dat"))

function distance(atom lat1, lon1, lat2, lon2, integer units)
    atom dist = 0
    if lat1!=lat2 or lon1!=lon2 then
        atom radlat1 = lat1*PI/180,
             radlat2 = lat2*PI/180,
             radtheta = (lon1-lon2)*PI/180
        dist = sin(radlat1)*sin(radlat2) + cos(radlat1)*cos(radlat2)*cos(radtheta)
        dist = arccos(min(dist,1))*180/PI * 60*1.1515576 -- in Statute Miles
        if units = 'K' then dist *= 1.609344 end if     -- in Kilometres
        if units = 'N' then dist *= 0.868976 end if    -- in Nautical Miles
    end if
    return dist
end function

function bearing(atom lat1, lon1, lat2, lon2)
    atom bear = NULL
    if lat1!=lat2 or lon1!=lon2 then
        atom radlat1 = lat1*PI/180,
             radlat2 = lat2*PI/180,
             raddlon = (lon2-lon1)*PI/180,
             y = sin(raddlon)*cos(radlat2),
             x = cos(radlat1)*sin(radlat2) - sin(radlat1)*cos(radlat2)*cos(raddlon)
        bear = remainder(atan2(y, x)*180/PI+360,360)
    end if
    return bear;
end function

procedure query(atom lat,lon, integer rows=20)
    sequence r = {}
    for a in airports do
        atom lat2 = a[Latitude],
             lon2 = a[Longitude],
             dist = round(distance(lat,lon,lat2,lon2,'N'),10),
             bear = round(bearing(lat,lon,lat2,lon2))
        r = append(r,{a[Name],a[Country],a[ICAO],dist,bear})
    end for
    r = sort_columns(r,{4})[1..rows]
    printf(1,"                Name                 |    Country     | ICAO | Distance in NM | Bearing in ° \n")
    printf(1,"-------------------------------------+----------------+------+----------------+--------------\n")
    printf(1,"%s\n",{join(r,"\n",fmt:=" %-36s| %-15s| %4s |           %4.1f |          %3d")})
    printf(1,"(%d rows)\n",rows)
end procedure

query(51.514669, 2.198581)

Output:

                Name                 |    Country     | ICAO | Distance in NM | Bearing in ° 
-------------------------------------+----------------+------+----------------+--------------
 Koksijde Air Base                   | Belgium        | EBFN |           30.7 |          146
 Ostend-Bruges International Airport | Belgium        | EBOS |           31.3 |          127
 Kent International Airport          | United Kingdom | EGMH |           33.5 |          252
 Calais-Dunkerque Airport            | France         | LFAC |           34.4 |          196
 Westkapelle heliport                | Belgium        | EBKW |           42.6 |          105
 Lympne Airport                      | United Kingdom | EGMK |           51.6 |          240
 Ursel Air Base                      | Belgium        | EBUL |           52.8 |          114
 Southend Airport                    | United Kingdom | EGMC |           56.2 |          274
 Merville-Calonne Airport            | France         | LFQT |           56.4 |          163
 Wevelgem Airport                    | Belgium        | EBKT |           56.5 |          137
 Midden-Zeeland Airport              | Netherlands    | EHMZ |           57.3 |           90
 Lydd Airport                        | United Kingdom | EGMD |           58.0 |          235
 RAF Wattisham                       | United Kingdom | EGUW |           59.0 |          309
 Beccles Airport                     | United Kingdom | EGSM |           59.3 |          339
 Lille/Marcq-en-Baroeul Airport      | France         | LFQO |           59.7 |          146
 Lashenden (Headcorn) Airfield       | United Kingdom | EGKH |           62.2 |          250
 Le Touquet-Côte d'Opale Airport     | France         | LFAT |           63.7 |          200
 Rochester Airport                   | United Kingdom | EGTO |           64.2 |          262
 Lille-Lesquin Airport               | France         | LFQQ |           66.2 |          149
 Thurrock Airfield                   | United Kingdom | EGMT |           68.4 |          272
(20 rows)

Python

''' Rosetta Code task Distance_and_Bearing '''

from math import radians, degrees, sin, cos, asin, atan2, sqrt
from pandas import read_csv


EARTH_RADIUS_KM = 6372.8
TASK_CONVERT_NM =  0.0094174
AIRPORT_DATA_FILE = 'airports.dat.txt'

QUERY_LATITUDE, QUERY_LONGITUDE = 51.514669, 2.198581


def haversine(lat1, lon1, lat2, lon2):
    '''
    Given two latitude, longitude pairs in degrees for two points on the Earth,
    get distance (nautical miles) and initial direction of travel (degrees)
    for travel from lat1, lon1 to lat2, lon2
    '''
    rlat1, rlon1, rlat2, rlon2 = [radians(x) for x in [lat1, lon1, lat2, lon2]]
    dlat = rlat2 - rlat1
    dlon = rlon2 - rlon1
    arc = sin(dlat / 2) ** 2 + cos(rlat1) * cos(rlat2) * sin(dlon / 2) ** 2
    clen = 2.0 * degrees(asin(sqrt(arc)))
    theta = atan2(sin(dlon) * cos(rlat2),
                  cos(rlat1) * sin(rlat2) - sin(rlat1) * cos(rlat2) * cos(dlon))
    theta = (degrees(theta) + 360) % 360
    return EARTH_RADIUS_KM * clen * TASK_CONVERT_NM, theta


def find_nearest_airports(latitude, longitude, wanted=20, csv=AIRPORT_DATA_FILE):
    ''' Given latitude and longitude, find `wanted` closest airports in database file csv. '''
    airports = read_csv(csv, header=None, usecols=[1, 3, 5, 6, 7], names=[
                        'Name', 'Country', 'ICAO', 'Latitude', 'Longitude'])
    airports['Distance'] = 0.0
    airports['Bearing'] = 0
    for (idx, row) in enumerate(airports.itertuples()):
        distance, bearing = haversine(
            latitude, longitude, row.Latitude, row.Longitude)
        airports.at[idx, 'Distance'] = round(distance, ndigits=1)
        airports.at[idx, 'Bearing'] = int(round(bearing))

    airports.sort_values(by=['Distance'], ignore_index=True, inplace=True)
    return airports.loc[0:wanted-1, ['Name', 'Country', 'ICAO', 'Distance', 'Bearing']]


print(find_nearest_airports(QUERY_LATITUDE, QUERY_LONGITUDE))

Output:

                                   Name         Country  ICAO  Distance  Bearing
0                     Koksijde Air Base         Belgium  EBFN      30.7      146
1   Ostend-Bruges International Airport         Belgium  EBOS      31.3      127
2            Kent International Airport  United Kingdom  EGMH      33.5      252
3              Calais-Dunkerque Airport          France  LFAC      34.4      196
4                  Westkapelle heliport         Belgium  EBKW      42.6      105
5                        Lympne Airport  United Kingdom  EGMK      51.6      240
6                        Ursel Air Base         Belgium  EBUL      52.8      114
7                      Southend Airport  United Kingdom  EGMC      56.2      274
8              Merville-Calonne Airport          France  LFQT      56.4      163
9                      Wevelgem Airport         Belgium  EBKT      56.5      137
10               Midden-Zeeland Airport     Netherlands  EHMZ      57.3       90
11                         Lydd Airport  United Kingdom  EGMD      58.0      235
12                        RAF Wattisham  United Kingdom  EGUW      59.0      309
13                      Beccles Airport  United Kingdom  EGSM      59.3      339
14       Lille/Marcq-en-Baroeul Airport          France  LFQO      59.7      146
15        Lashenden (Headcorn) Airfield  United Kingdom  EGKH      62.2      250
16      Le Touquet-Côte d'Opale Airport          France  LFAT      63.7      200
17                    Rochester Airport  United Kingdom  EGTO      64.2      262
18                Lille-Lesquin Airport          France  LFQQ      66.2      149
19                    Thurrock Airfield  United Kingdom  EGMT      68.4      272

Raku

use Text::CSV;

constant $EARTH_RADIUS_IN_NAUTICAL_MILES = 6372.8 / 1.852;

sub degrees   ( Real $rad ) { $rad / tau * 360 }
sub radians   ( Real $deg ) { $deg * tau / 360 }
sub haversine ( Real $x   ) { ($x / 2).sin.²   }
sub arc_haver ( Real $x   ) { $x.sqrt.asin * 2 }

sub great_circle_distance ( \φ1, \λ1, \φ2, \λ2 ) {
    # https://en.wikipedia.org/wiki/Haversine_formula
    #   latitude (represented by φ) and longitude (represented by λ)
    #   hav(θ) = hav(φ₂ − φ₁) + cos(φ₁) cos(φ₂) hav(λ₂ − λ₁)
    arc_haver(
        haversine(φ2 - φ1)
      + haversine(λ2 - λ1) * cos(φ1) * cos(φ2)
    );
}

sub great_circle_bearing ( \φ1, \λ1, \φ2, \λ2 ) {
    atan2(                          φ2.cos * (λ2 - λ1).sin,
        φ1.cos * φ2.sin - φ1.sin * φ2.cos * (λ2 - λ1).cos,
    );
}

sub distance_and_bearing ( $lat1, $lon1, $lat2, $lon2 ) {
    my @ll = map &radians, $lat1, $lon1, $lat2, $lon2;

    my $dist  = great_circle_distance(|@ll);
    my $theta = great_circle_bearing( |@ll);

    return  $dist * $EARTH_RADIUS_IN_NAUTICAL_MILES,
            degrees( $theta + (tau if $theta < 0) )
}

sub find_nearest_airports ( $latitude, $longitude, $csv_path ) {
    my &d_and_b_from_here = &distance_and_bearing.assuming($latitude, $longitude);

    my @airports = csv(
        in => $csv_path,
        headers => [<ID Name City Country IATA ICAO Latitude Longitude>],
    );

    for @airports -> %row {
        %row<Distance Bearing> = d_and_b_from_here( +%row<Latitude>, +%row<Longitude> );
    }

    return @airports.sort(*.<Distance>);
}

sub MAIN ( $lat = 51.514669, $long = 2.198581, $wanted = 20, $csv = 'airports.dat' ) {
    printf "%7s\t%7s\t%-7s\t%-15s\t%s\n", <Dist Bear ICAO Country Name>;

    for find_nearest_airports($lat, $long, $csv).head($wanted) {
        printf "%7.1f\t    %03d\t%-7s\t%-15s\t%s\n",
            .<Distance Bearing ICAO Country Name>;
    }
}

Output:

   Dist	   Bear	ICAO   	Country        	Name
   30.7	    146	EBFN   	Belgium        	Koksijde Air Base
   31.3	    127	EBOS   	Belgium        	Ostend-Bruges International Airport
   33.6	    252	EGMH   	United Kingdom 	Kent International Airport
   34.4	    195	LFAC   	France         	Calais-Dunkerque Airport
   42.6	    105	EBKW   	Belgium        	Westkapelle heliport
   51.6	    240	EGMK   	United Kingdom 	Lympne Airport
   52.8	    114	EBUL   	Belgium        	Ursel Air Base
   56.2	    274	EGMC   	United Kingdom 	Southend Airport
   56.4	    162	LFQT   	France         	Merville-Calonne Airport
   56.5	    137	EBKT   	Belgium        	Wevelgem Airport
   57.3	    089	EHMZ   	Netherlands    	Midden-Zeeland Airport
   58.0	    235	EGMD   	United Kingdom 	Lydd Airport
   59.0	    309	EGUW   	United Kingdom 	RAF Wattisham
   59.3	    339	EGSM   	United Kingdom 	Beccles Airport
   59.7	    146	LFQO   	France         	Lille/Marcq-en-Baroeul Airport
   62.2	    250	EGKH   	United Kingdom 	Lashenden (Headcorn) Airfield
   63.7	    200	LFAT   	France         	Le Touquet-Côte d'Opale Airport
   64.2	    261	EGTO   	United Kingdom 	Rochester Airport
   66.3	    149	LFQQ   	France         	Lille-Lesquin Airport
   68.4	    271	EGMT   	United Kingdom 	Thurrock Airfield

SQL/PostgreSQL

Create table and copy from URL.

-- create table airports with 14 columns
CREATE TABLE airports (
    Airport_ID serial PRIMARY KEY,
    Name VARCHAR NOT NULL,
    City VARCHAR,
    Country VARCHAR NOT NULL,
    IATA VARCHAR,
    ICAO VARCHAR,
    Latitude double precision NOT NULL,
    Longitude double precision NOT NULL,
    Altitude SMALLINT,
    Timezone VARCHAR,
    DST VARCHAR,
    Tz_Olson VARCHAR,
    Type VARCHAR,
    Source VARCHAR
);   

-- copy CSV airports.dat from URL 
COPY airports FROM 
PROGRAM 'curl "https://raw.githubusercontent.com/jpatokal/openflights/master/data/airports.dat"' 
WITH (FORMAT csv);

Functions for distance and bearing.

-- calculate distance
CREATE OR REPLACE FUNCTION calculate_distance(lat1 float, lon1 float, lat2 float, lon2 float, units varchar)
RETURNS numeric AS $dist$
    DECLARE
        dist float = 0;
        radlat1 float;
        radlat2 float;
        theta float;
        radtheta float;
    BEGIN
        IF lat1 = lat2 AND lon1 = lon2
            THEN RETURN dist;
        ELSE
            radlat1 = pi() * lat1 / 180;
            radlat2 = pi() * lat2 / 180;
            theta = lon1 - lon2;
            radtheta = pi() * theta / 180;
            dist = sin(radlat1) * sin(radlat2) + cos(radlat1) * cos(radlat2) * cos(radtheta);

            IF dist > 1 THEN dist = 1; END IF;

            dist = acos(dist);
            dist = dist * 180 / pi();

            -- Distance in Statute Miles
            dist = dist * 60 * 1.1515576; 
            
            -- Distance in Kilometres
            IF units = 'K' THEN dist = dist * 1.609344; END IF;

            -- Distance in Nautical Miles
            IF units = 'N' THEN dist = dist * 0.868976; END IF;

            dist = dist::numeric;
            RETURN dist;
        END IF;
    END;
$dist$ LANGUAGE plpgsql;


-- calculate bearing
CREATE OR REPLACE FUNCTION calculate_bearing(lat1 float, lon1 float, lat2 float, lon2 float)
RETURNS numeric AS $bear$
    DECLARE
        bear float = NULL;
        radlat1 float;
        radlat2 float;
        raddlon float;
        y float;
        x float;
        
    BEGIN
        IF lat1 = lat2 AND lon1 = lon2
            THEN RETURN bear;
        ELSE
            radlat1 = pi() * lat1 / 180;
            radlat2 = pi() * lat2 / 180;
            raddlon = pi() * (lon2 - lon1) / 180;

            y = sin(raddlon) * cos(radlat2);
            x = cos(radlat1) * sin(radlat2) - sin(radlat1) * cos(radlat2) * cos(raddlon);
            
            bear = atan2(y, x) * 180 / pi();
            bear = (bear::numeric + 360) % 360;
       
            RETURN bear;
        END IF;
    END;
$bear$ LANGUAGE plpgsql;

Request from airplane at position ( 51.514669, 2.198581 ).

Select 
   Name "Name",
   Country "Country",
   ICAO "ICAO",
   ROUND(calculate_distance(51.514669, 2.198581, Latitude, Longitude, 'N'), 1) "Distance in NM",
   ROUND(calculate_bearing(51.514669, 2.198581, Latitude, Longitude), 0) "Bearing in °"
From 
    airports
ORDER BY "Distance in NM"
LIMIT 20;

Output:

                Name                 |    Country     | ICAO | Distance in NM | Bearing in ° 
-------------------------------------+----------------+------+----------------+--------------
 Koksijde Air Base                   | Belgium        | EBFN |           30.7 |          146
 Ostend-Bruges International Airport | Belgium        | EBOS |           31.3 |          127
 Kent International Airport          | United Kingdom | EGMH |           33.5 |          252
 Calais-Dunkerque Airport            | France         | LFAC |           34.4 |          196
 Westkapelle heliport                | Belgium        | EBKW |           42.6 |          105
 Lympne Airport                      | United Kingdom | EGMK |           51.6 |          240
 Ursel Air Base                      | Belgium        | EBUL |           52.8 |          114
 Southend Airport                    | United Kingdom | EGMC |           56.2 |          274
 Merville-Calonne Airport            | France         | LFQT |           56.4 |          163
 Wevelgem Airport                    | Belgium        | EBKT |           56.5 |          137
 Midden-Zeeland Airport              | Netherlands    | EHMZ |           57.3 |           90
 Lydd Airport                        | United Kingdom | EGMD |           58.0 |          235
 RAF Wattisham                       | United Kingdom | EGUW |           59.0 |          309
 Beccles Airport                     | United Kingdom | EGSM |           59.3 |          339
 Lille/Marcq-en-Baroeul Airport      | France         | LFQO |           59.7 |          146
 Lashenden (Headcorn) Airfield       | United Kingdom | EGKH |           62.2 |          250
 Le Touquet-Côte d'Opale Airport     | France         | LFAT |           63.7 |          200
 Rochester Airport                   | United Kingdom | EGTO |           64.2 |          262
 Lille-Lesquin Airport               | France         | LFQQ |           66.2 |          149
 Thurrock Airfield                   | United Kingdom | EGMT |           68.4 |          272
(20 rows)

Wren

Translation of: SQL

Library: Wren-dynamic

Library: Wren-str

Library: Wren-fmt

However, rather than use Wren-sql, we program it so that it can be run from Wren-CLI. Runtime is about 0.24 seconds.

import "io" for File
import "./dynamic" for Tuple
import "./str" for Str
import "./fmt" for Fmt

var airportFields = [
    "airportID",
    "name",
    "city",
    "country",
    "iata",
    "icao",
    "latitude",
    "longitude",
    "altitude",
    "timezone",
    "dst",
    "tzOlson",
    "type",
    "source"
]
var Airport = Tuple.create("Airport", airportFields)

var fileName = "airports.dat" // local copy
var lines = File.read(fileName).trimEnd().split("\n")
var lc = lines.count
var airports = List.filled(lc, null)
for (i in 0...lc) {
    var fields    = lines[i].split(",").map { |f| f.replace("\"", "") }.toList
    if (fields.count > 14) {  // there are field(s) with embedded comma(s)
        fields    = Str.splitCsv(lines[i])
    }
    var airportID = Num.fromString(fields[0])
    var name      = fields[1]
    var city      = fields[2]
    var country   = fields[3]
    var iata      = fields[4]
    var icao      = fields[5]
    var latitude  = Num.fromString(fields[6])
    var longitude = Num.fromString(fields[7])
    var altitude  = Num.fromString(fields[8])
    var timezone  = fields[9]
    var dst       = fields[10]
    var tzOlson   = fields[11]
    var type      = fields[12]
    var source    = fields[13]
    airports[i]   = Airport.new(airportID, name, city, country, iata, icao, latitude, longitude,
                                altitude, timezone, dst, tzOlson, type, source)
}

var calculateDistance = Fn.new { |lat1, lon1, lat2, lon2, units|
    if (lat1 == lat2 && lon1 == lon2) return 0
    var radlat1 = Num.pi * lat1 / 180
    var radlat2 = Num.pi * lat2 / 180
    var theta = lon1 - lon2
    var radtheta = Num.pi * theta / 180
    var dist = radlat1.sin * radlat2.sin + radlat1.cos * radlat2.cos * radtheta.cos
    if (dist > 1) dist = 1
    dist = dist.acos * 180 / Num.pi * 60 * 1.1515576  // distance in statute miles
    if (units == "K") dist = dist * 1.609344          // distance in kilometers
    if (units == "N") dist = dist * 0.868976          // distance in nautical miles
    return dist
}

var calculateBearing = Fn.new { |lat1, lon1, lat2, lon2|
    if (lat1 == lat2 && lon1 == lon2) return 0
    var radlat1 = Num.pi * lat1 / 180
    var radlat2 = Num.pi * lat2 / 180
    var raddlon = Num.pi * (lon2 - lon1) / 180
    var y = raddlon.sin * radlat2.cos
    var x = radlat1.cos * radlat2.sin - radlat1.sin * radlat2.cos * raddlon.cos
    var bear = y.atan(x) * 180 / Num.pi
    return (bear + 360) % 360
}

// request from airplane at position (51.514669, 2.198581)
var query = List.filled(airports.count, null)
for (i in 0...airports.count) {
    var a = airports[i]    
    var dist = calculateDistance.call(51.514669, 2.198581, a.latitude, a.longitude, "N")
    dist = (dist * 10).round / 10
    var bear = calculateBearing.call(51.514669, 2.198581, a.latitude, a.longitude).round
    query[i] = [a.name, a.country, a.icao, dist, bear]
}
query.sort { |q1, q2| q1[3] < q2[3] }
System.print("                Name                 |    Country     | ICAO | Distance in NM | Bearing in ° ")
System.print("-------------------------------------+----------------+------+----------------+--------------")
var fmt = " $-36s| $-15s| $4s |           $4.1f |          $3d"
for (i in 0..19) Fmt.lprint(fmt, query[i])
System.print("(20 rows)")

Output:

                Name                 |    Country     | ICAO | Distance in NM | Bearing in ° 
-------------------------------------+----------------+------+----------------+--------------
 Koksijde Air Base                   | Belgium        | EBFN |           30.7 |          146
 Ostend-Bruges International Airport | Belgium        | EBOS |           31.3 |          127
 Kent International Airport          | United Kingdom | EGMH |           33.5 |          252
 Calais-Dunkerque Airport            | France         | LFAC |           34.4 |          196
 Westkapelle heliport                | Belgium        | EBKW |           42.6 |          105
 Lympne Airport                      | United Kingdom | EGMK |           51.6 |          240
 Ursel Air Base                      | Belgium        | EBUL |           52.8 |          114
 Southend Airport                    | United Kingdom | EGMC |           56.2 |          274
 Merville-Calonne Airport            | France         | LFQT |           56.4 |          163
 Wevelgem Airport                    | Belgium        | EBKT |           56.5 |          137
 Midden-Zeeland Airport              | Netherlands    | EHMZ |           57.3 |           90
 Lydd Airport                        | United Kingdom | EGMD |           58.0 |          235
 RAF Wattisham                       | United Kingdom | EGUW |           59.0 |          309
 Beccles Airport                     | United Kingdom | EGSM |           59.3 |          339
 Lille/Marcq-en-Baroeul Airport      | France         | LFQO |           59.7 |          146
 Lashenden (Headcorn) Airfield       | United Kingdom | EGKH |           62.2 |          250
 Le Touquet-Côte d'Opale Airport     | France         | LFAT |           63.7 |          200
 Rochester Airport                   | United Kingdom | EGTO |           64.2 |          262
 Lille-Lesquin Airport               | France         | LFQQ |           66.2 |          149
 Thurrock Airfield                   | United Kingdom | EGMT |           68.4 |          272
(20 rows)