Distance and Bearing
You are encouraged to solve this task according to the task description, using any language you may know.
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
C++
#include <cmath>
#include <cstdint>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <numbers>
#include <regex>
#include <sstream>
#include <string>
#include <vector>
constexpr double RADIUS = 6'371 / 1.852; // Mean radius of the Earth in nautical miles
constexpr double RADIAN_TO_DEGREE = 180.0 / std::numbers::pi;
class airport {
public:
airport(std::string aName, std::string aCountry, std::string aIcao, double aLatitude, double aLongitude)
: name(aName), country(aCountry), icao(aIcao), latitude(aLatitude), longitude(aLongitude) {}
std::string name;
std::string country;
std::string icao;
double latitude;
double longitude;
};
// Convert the given string to a double, which represents an angle,
// and then convert the angle from degrees to radians
double to_double_radians(const std::string& text) {
std::istringstream stream(text);
double decimal = 0.0;
stream >> decimal;
return decimal / RADIAN_TO_DEGREE;
}
std::string do_replace(const std::string& text, const std::string& original, const std::string& replacement) {
return std::regex_replace(text, std::regex(original), replacement);
}
std::vector<std::string> split(const std::string& line, const char& delimiter) {
std::stringstream stream(line);
std::string item;
std::vector<std::string> items;
while ( std::getline(stream, item, delimiter) ) {
items.push_back(std::move(item));
}
return items;
}
void read_file(const std::string& file_name, std::vector<airport>& airports) {
std::ifstream airports_file(file_name);
std::string line;
while ( std::getline(airports_file, line) ) {
std::vector<std::string> sections = split(line, ',');
airport air_port(do_replace(sections[1], "\"", ""), // Remove the double quotes from the string
do_replace(sections[3], "\"", ""),
do_replace(sections[5], "\"", ""),
to_double_radians(sections[6]),
to_double_radians(sections[7]));
airports.push_back(std::move(air_port));
}
airports_file.close();
}
// The given angles are in radians, and the result is in nautical miles.
double distance(double phi1, double lambda1, double phi2, double lambda2) {
double a = pow(sin((phi2 - phi1) * 0.5), 2) + cos(phi1) * cos(phi2) * pow(sin((lambda2 - lambda1) * 0.5), 2);
double c = 2 * atan2(sqrt(a), sqrt(1 - a));
return RADIUS * c;
}
// The given angles are in radians, and the result is in degrees in the range [0, 360).
double bearing(double phi1, double lambda1, double phi2, double lambda2) {
double delta = lambda2 - lambda1;
double result = atan2(sin(delta) * cos(phi2), cos(phi1) * sin(phi2) - sin(phi1) * cos(phi2) * cos(delta));
return std::fmod(result * RADIAN_TO_DEGREE + 360.0, 360.0);
}
int main() {
std::vector<airport> airports;
read_file("airports.dat", airports);
const double plane_latitude = 51.514669 / RADIAN_TO_DEGREE;
const double plane_longitude = 2.198581 / RADIAN_TO_DEGREE;
std::vector<std::pair<double, uint64_t>> distances;
for ( uint64_t i = 0; i < airports.size(); ++i ) {
double dist = distance(plane_latitude, plane_longitude, airports[i].latitude, airports[i].longitude);
distances.push_back(std::make_pair(dist, i));
}
std::sort(distances.begin(), distances.end(),
[](auto& left, auto& right) { return left.first < right.first; });
std::cout << "Distance" << std::setw(9) << "Bearing" << std::setw(11) << "ICAO"
<< std::setw(20) << "Country" << std::setw(40) << "Airport" << std::endl;
std::cout << std::string(88, '-') << std::endl;
for ( uint32_t i = 0; i < 20; ++i ) {
auto[distance, index] = distances[i];
airport air_port = airports[index];
double bear = bearing(plane_latitude, plane_longitude, air_port.latitude, air_port.longitude);
std::cout << std::setw(8) << std::fixed << std::setprecision(1) << distance
<< std::setw(9) << std::setprecision(0) << std::round(bear)
<< std::setw(11) << air_port.icao << std::setw(20) << air_port.country
<< std::setw(40) << air_port.name << std::endl;
}
}
- Output:
Distance Bearing ICAO Country Airport
----------------------------------------------------------------------------------------
30.7 146 EBFN Belgium Koksijde Air Base
31.3 127 EBOS Belgium Ostend-Bruges International Airport
33.5 252 EGMH United Kingdom Kent International Airport
34.4 196 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 163 LFQT France Merville-Calonne Airport
56.5 137 EBKT Belgium Wevelgem Airport
57.3 90 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 262 EGTO United Kingdom Rochester Airport
66.2 149 LFQQ France Lille-Lesquin Airport
68.4 272 EGMT United Kingdom Thurrock Airfield
Fortran
program code_translation
use iso_fortran_env , only:real64 , real32
use helpers
implicit none
integer :: i
character(len = 300) :: input_string
character(len = *) , parameter :: filnam = "airports.csv"
character(len = *) , parameter :: delimiter = ","
character(len = 100) :: tokens(20)
integer :: num_tokens
integer :: iostat , counter
type(airport) , allocatable :: airports(:)
type(whereme) :: location
integer :: ii , jj
real(real64) :: dist
location%latitude = (51.514669D0)!Rosetta co-ords
location%longitude = (2.19858D0)
!Nearest to latitude 51.51467,longitude 2.19858 degrees
!
!Open the file for reading
open(unit = 10 , file = 'airports.csv' , status = 'OLD' , &
&action = 'READ' , iostat = iostat , encoding = 'UTF-8')
if( iostat/=0 )stop 'Could not open file'
counter = 0
do
read(10 , '(a)' , end = 100)input_string
counter = counter + 1
end do
100 continue ! Now we know how many elements there are
rewind(unit = 10)
write(*, '(a,1x,i0,1x,a)') 'Scanning',counter,'lines'
allocate(airports(counter))
call system_clock(count = ii)
do i = 1 , counter
read(10 , '(a)' , end = 200)input_string
call tokenizes(input_string , tokens , num_tokens , ',')
read(tokens(1) , *)airports(i)%airportid
airports(i)%name = trim(adjustl(tokens(2)))
airports(i)%city = trim(adjustl(tokens(3)))
airports(i)%country = trim(adjustl(tokens(4)))
airports(i)%iata = trim(adjustl(tokens(5)))
airports(i)%icao = trim(adjustl(tokens(6)))
read(tokens(7) , '(f20.16)')airports(i)%latitude
read(tokens(8) , '(F20.16)')airports(i)%longitude
read(tokens(9) , *)airports(i)%altitude
read(tokens(10) , '(F5.2)' )airports(i)%timezone
airports(i)%dst = trim(adjustl(tokens(11)))
airports(i)%tzolson = trim(adjustl(tokens(12)))
airports(i)%typez = trim(adjustl(tokens(13)))
airports(i)%source = trim(adjustl(tokens(14)))
! Calculate the distance and bearing straight away
airports(i)%distance = haversine(location%latitude , &
& location%longitude , airports(i)%latitude ,&
& airports(i)%longitude)
airports(i)%bearing = bearing(location%latitude , &
& location%longitude , airports(i)%latitude , &
& airports(i)%longitude)
end do
200 continue
call system_clock(count = jj)
write(* , *) 'Read complete, time taken = ' , (jj - ii) , &
& 'milliseconds' // char(10) // char(13)
call sortem(airports) ! Sort the airports out
write(*, '(/,2x,a,t14,a,t75,a,t95,a,t108,a,t117,a)') 'Num' , 'Name' , &
&'Country' , 'ICAO' , 'Dist.' , 'Bearing'
write(*, '(a)') repeat('=' , 130)
do jj = 1 , 20!First 20 only
write(*, '(i5,t8,a,t75,a,t95,a,t105,f8.1,t117,i0)') airports(jj) &
& %airportid , airports(jj)%name , airports(jj)%country , &
& airports(jj)%icao , airports(jj)%distance , &
& nint(airports(jj)%bearing)
end do
stop 'Normal completion' // char(10) // char(13)
end program code_translation
module helpers
use iso_fortran_env , only:real32,real64
implicit none
real(real64) , parameter :: radius_in_km = 6372.8D0
real(real64) , parameter :: kilos_to_nautical = 0.5399568D0
type whereme
real(real64) :: latitude , longitude
end type whereme
type airport
integer :: airportid
character(len = 100) :: name
character(len = 50) :: city , country
character(len = 10) :: iata , icao
real(real64) :: latitude , longitude
integer :: altitude
real(real32) :: timezone
character(len = 10) :: dst
character(len = 100) :: tzolson
character(len = 20) :: typez , source
real(real64) :: distance , bearing
end type airport
contains ! We'll calculate them and store in each airport
!
! The given angles are in radians, and the result is in degrees in the range [0, 360).
function bearing(lat1 , lon1 , lat2 , lon2)
real(real64) , parameter :: toradians = acos( - 1.0D0)/180.0D0
real(real64) :: bearing
real(real64) , intent(in) :: lat1
real(real64) , intent(in) :: lon1
real(real64) , intent(in) :: lat2
real(real64) , intent(in) :: lon2
real(real64) :: dlat
real(real64) :: dlon
real(real64) :: rlat1
real(real64) :: rlat2
real(real64) :: x
real(real64) :: y
!
dlat = (lat2 - lat1)*toradians
dlon = toradians*(lon2 - lon1)
rlat1 = toradians*(lat1)
rlat2 = toradians*(lat2)
!
x = cos(rlat2)*sin(dlon)
y = cos(rlat1)*sin(rlat2) - sin(rlat1)*cos(rlat2)*cos(dlon)
bearing = atan2(x , y)
bearing = to_degrees(bearing)
bearing = mod(bearing + 360.0D0 , 360.0D0)
end function bearing
!
!
function to_radian(degree) result(rad)
real(real64) , parameter :: deg_to_rad = atan(1.0D0)/45.0D0 ! exploit intrinsic atan to generate pi/180 runtime constant
!
real(real64) :: rad
real(real64) , intent(in) :: degree
! degrees to radians
rad = degree*deg_to_rad
end function to_radian
!
function to_degrees(radians) result(deg)
real(real64) , parameter :: radian_to_degree = 180.0D0/acos( - 1.0D0)
!
real(real64) :: deg
real(real64) , intent(in) :: radians
deg = radians*radian_to_degree
end function to_degrees
!
function haversine(deglat1 , deglon1 , deglat2 , deglon2) result(dist)
real(real64) :: dist
real(real64) , intent(in) :: deglat1
real(real64) , intent(in) :: deglon1
real(real64) , intent(in) :: deglat2
real(real64) , intent(in) :: deglon2
real(real64) :: a
real(real64) :: c
real(real64) :: dlat
real(real64) :: dlon
real(real64) :: lat1
real(real64) :: lat2
! great circle distance
dlat = to_radian(deglat2 - deglat1)
dlon = to_radian(deglon2 - deglon1)
lat1 = to_radian(deglat1)
lat2 = to_radian(deglat2)
a = (sin(dlat/2.0D0)**2) + cos(lat1)*cos(lat2) &
& *((sin(dlon/2.0D0))**2)
c = 2.0D0*asin(sqrt(a))
dist = radius_in_km*c*kilos_to_nautical
end function haversine
subroutine sortem(airports)!Bubble sort them, nice and easy
type(airport) , intent(inout) , dimension(:) :: airports
integer :: i
integer :: k
logical :: swapped
type(airport) :: temp
swapped = .true.
k = size(airports%distance)
do while ( swapped )
swapped = .false.
do i = 1 , k - 1
if( airports(i)%distance>airports(i + 1)%distance )then
temp = airports(i)
airports(i) = airports(i + 1)
airports(i + 1) = temp
swapped = .true.
end if
end do
end do
return
end subroutine sortem
!
subroutine tokenizes(input_string , tokens , num_tokens , delimiter)
character(*) , intent(in) :: input_string
character(*) , intent(out) , dimension(:) :: tokens
integer , intent(out) :: num_tokens
character(1) , intent(in) :: delimiter
integer :: end_idx
integer :: i
integer :: start_idx
character(100) :: temp_string
num_tokens = 0
temp_string = trim(input_string)
start_idx = 1
do i = 1 , len_trim(temp_string)
if( (temp_string(i:i)==delimiter) )then
end_idx = i - 1
if( end_idx>=start_idx )then
num_tokens = num_tokens + 1
tokens(num_tokens) = ''
tokens(num_tokens) = temp_string(start_idx:end_idx)
end if
start_idx = i + 1
end if
end do
! Handle the last token
if( start_idx<=len_trim(temp_string) )then
num_tokens = num_tokens + 1
tokens(num_tokens) = temp_string(start_idx:len_trim( &
& temp_string))
end if
end subroutine tokenizes
!
end module helpers
!
- Output:
Scanning 7698 lines Read complete, time taken = 93 milliseconds Num Name Country ICAO Dist. Bearing ================================================================================================================================== 306 Koksijde Air Base Belgium EBFN 30.7 146 310 Ostend-Bruges International Airport Belgium EBOS 31.3 127 510 Kent International Airport United Kingdom EGMH 33.6 252 1254 Calais-Dunkerque Airport France LFAC 34.4 196 7810 Westkapelle heliport Belgium EBKW 42.6 105 6886 Lympne Airport United Kingdom EGMK 51.6 240 314 Ursel Air Base Belgium EBUL 52.8 114 508 Southend Airport United Kingdom EGMC 56.2 274 1400 Merville-Calonne Airport France LFQT 56.4 163 308 Wevelgem Airport Belgium EBKT 56.5 137 7836 Midden-Zeeland Airport Netherlands EHMZ 57.3 90 509 Lydd Airport United Kingdom EGMD 58.0 235 558 RAF Wattisham United Kingdom EGUW 59.0 309 8965 Beccles Airport United Kingdom EGSM 59.3 339 10089 Lille/Marcq-en-Baroeul Airport France LFQO 59.7 146 10745 Lashenden (Headcorn) Airfield United Kingdom EGKH 62.2 250 1259 Le Touquet-Côte d'Opale Airport France LFAT 63.7 200 8894 Rochester Airport United Kingdom EGTO 64.2 262 1399 Lille-Lesquin Airport France LFQQ 66.3 149 10747 Thurrock Airfield United Kingdom EGMT 68.4 272
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}}°
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
Java
Author Sjharper79 (talk)
Description
This java program contains two classes. One class holds the airport information, and the other class does all of the work to perform calculations and print results. Some of the logic adapted from Python below.
Implementation
// The Airport class holds each airport object
package distanceAndBearing;
public class Airport {
private String airport;
private String country;
private String icao;
private double lat;
private double lon;
public String getAirportName() { return this.airport; }
public void setAirportName(String airport) { this.airport = airport; }
public String getCountry() { return this.country; }
public void setCountry(String country) { this.country = country; }
public String getIcao() { return this.icao; }
public void setIcao(String icao) { this.icao = icao; }
public double getLat() { return this.lat; }
public void setLat(double lat) { this.lat = lat; }
public double getLon() { return this.lon; }
public void setLon(double lon) { this.lon = lon; }
@Override
public String toString() {return "Airport: " + getAirportName() + ": ICAO: " + getIcao();}
}
// The DistanceAndBearing class does all the work.
package distanceAndBearing;
import java.io.File;
import java.text.DecimalFormat;
import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Scanner;
import java.util.stream.Collectors;
public class DistanceAndBearing {
private final double earthRadius = 6371;
private File datFile;
private List<Airport> airports;
public DistanceAndBearing() { this.airports = new ArrayList<Airport>(); }
public boolean readFile(String filename) {
this.datFile = new File(filename);
try {
Scanner fileScanner = new Scanner(datFile);
String line;
while (fileScanner.hasNextLine()) {
line = fileScanner.nextLine();
line = line.replace(", ", "; "); // There are some commas in airport names in the dat
// file
line = line.replace(",\",\"", "\",\""); // There are some airport names that
// end in a comma in the dat file
String[] parts = line.split(",");
Airport airport = new Airport();
airport.setAirportName(parts[1].replace("\"", "")); // Remove the double quotes from the string
airport.setCountry(parts[3].replace("\"", "")); // Remove the double quotes from the string
airport.setIcao(parts[5].replace("\"", "")); // Remove the double quotes from the string
airport.setLat(Double.valueOf(parts[6]));
airport.setLon(Double.valueOf(parts[7]));
this.airports.add(airport);
}
fileScanner.close();
return true; // Return true if the read was successful
} catch (Exception e) {
e.printStackTrace();
return false; // Return false if the read was not successful
}}
public double[] calculate(double lat1, double lon1, double lat2, double lon2) {
double[] results = new double[2];
double dLat = Math.toRadians(lat2 - lat1);
double dLon = Math.toRadians(lon2 - lon1);
double rlat1 = Math.toRadians(lat1);
double rlat2 = Math.toRadians(lat2);
double a = Math.pow(Math.sin(dLat / 2), 2)
+ Math.pow(Math.sin(dLon / 2), 2) * Math.cos(rlat1) * Math.cos(rlat2);
double c = 2 * Math.asin(Math.sqrt(a));
double distance = earthRadius * c;
DecimalFormat df = new DecimalFormat("#0.00");
distance = Double.valueOf(df.format(distance));
results[0] = distance;
double X = Math.cos(rlat2) * Math.sin(dLon);
double Y = Math.cos(rlat1) * Math.sin(rlat2) - Math.sin(rlat1) * Math.cos(rlat2) * Math.cos(dLon);
double heading = Math.atan2(X, Y);
heading = Math.toDegrees(heading);
results[1] = heading;
return results;
}
public Airport searchByName(final String name) {
Airport airport = new Airport();
List<Airport> results = this.airports.stream().filter(ap -> ap.getAirportName().contains(name))
.collect(Collectors.toList());
airport = results.get(0);
return airport;
}
public List<Airport> findClosestAirports(double lat, double lon) {
// TODO Auto-generated method stub
Map<Double, Airport> airportDistances = new HashMap<>();
Map<Double, Airport> airportHeading = new HashMap<>();
List<Airport> closestAirports = new ArrayList<Airport>();
// This loop finds the distance and heading for every airport and saves them
// into two separate Maps
for (Airport ap : this.airports) {
double[] result = calculate(lat, lon, ap.getLat(), ap.getLon());
airportDistances.put(result[0], ap);
airportHeading.put(result[1], ap);
}
// Get the keyset from the distance map and sort it.
ArrayList<Double> distances = new ArrayList<>(airportDistances.keySet());
Collections.sort(distances);
// Get the first 20 airports by finding the value in the distance map for
// each distance in the sorted Arraylist. Then get the airport name, and
// use that to search for the airport in the airports List.
// Save that into a new List
for (int i = 0; i < 20; i++) { closestAirports.add(searchByName((airportDistances.get(distances.get(i)).getAirportName())));}
// Find the distance and heading for each of the top 20 airports.
Map<String, Double> distanceMap = new HashMap<>();
for (Double d : airportDistances.keySet()) { distanceMap.put(airportDistances.get(d).getAirportName(), d);}
Map<String, Double> headingMap = new HashMap<>();
for (Double d : airportHeading.keySet()) {
double d2 = d;
if(d2<0){d2+=360'}
headingMap.put(airportHeading.get(d).getAirportName(), d2); }
// Print the results.
System.out.printf("%-4s %-40s %-25s %-6s %12s %15s\n", "Num", "Airport", "Country", "ICAO", "Distance", "Bearing");
System.out.println("-----------------------------------------------------------------------------------------------------------");
int i = 0;
for (Airport a : closestAirports) {
System.out.printf("%-4s %-40s %-25s %-6s %12.1f %15.0f\n", ++i, a.getAirportName(), a.getCountry(), a.getIcao(), distanceMap.get(a.getAirportName())*0.5399568, headingMap.get(a.getAirportName()));
}
return closestAirports;
}
}
Usage
import distanceAndBearing.DistanceAndBearing;
public class MyClass {
public static void main(String[] args) {
DistanceAndBearing dandb = new DistanceAndBearing();
dandb.readFile("airports.txt");
dandb.findClosestAirports(51.514669,2.198581);
}
}
Output when lat = 51.514669 and lon = 2.198581
Num Airport Country ICAO Distance Bearing ----------------------------------------------------------------------------------------------------------- 1 Koksijde Air Base Belgium EBFN 30.7 146 2 Ostend-Bruges International Airport 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.4 163 10 Wevelgem Airport Belgium EBKT 56.5 137 11 Midden-Zeeland Airport Netherlands EHMZ 57.3 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
jq
Adapted from Wren
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
Nim
import std/[algorithm, math, parsecsv, strformat, strutils]
const
R = 6_371 / 1.852 # Mean radius of Earth in nautic miles.
PlaneLat = 51.514669
PlaneLong = 2.198581
type
Airport = object
name: string
country: string
icao: string
latitude: float
longitude: float # In radians.
altitude: float # In radians.
func distance(𝜑₁, λ₁, 𝜑₂, λ₂: float): float =
## Return the distance computed using the Haversine formula.
## Angles are in radians. The result is expressed in nautic miles.
let a = sin((𝜑₂ - 𝜑₁) * 0.5)^2 + cos(𝜑₁) * cos(𝜑₂) * sin((λ₂ - λ₁) * 0.5)^2
let c = 2 * arctan2(sqrt(a), sqrt(1 - a))
result = R * c
func bearing(𝜑₁, λ₁, 𝜑₂, λ₂: float): float =
## Return the bearing.
## Angles are in radians. The result is expressed in degrees in range [0..360[.
let Δλ = λ₂ - λ₁
result = arctan2(sin(Δλ) * cos(𝜑₂), cos(𝜑₁) * sin(𝜑₂) - sin(𝜑₁) * cos(𝜑₂) * cos(Δλ))
result = (result.radToDeg + 360) mod 360
# Parse the "airports.dat" file.
var parser: CsvParser
var airports: seq[Airport]
parser.open("airports.dat")
while parser.readRow():
assert parser.row.len == 14
airports.add Airport(name: parser.row[1],
country: parser.row[3],
icao: parser.row[5],
latitude: parser.row[6].parseFloat().degToRad,
longitude: parser.row[7].parseFloat().degToRad)
# Compute the distances then sort them keeping the airport index.
type Distances = tuple[value: float; index: int]
var distances : seq[Distances]
let 𝜑₁ = PlaneLat.degToRad
let λ₁ = PlaneLong.degToRad
for i, airport in airports:
distances.add (distance(𝜑₁, λ₁, airport.latitude, airport.longitude), i)
distances.sort(Ascending)
# Display the result for the 20 nearest airports.
echo &"""{"Airport":<40}{"Country":<20}{"ICAO":<9}{"Distance":<9}{"Bearing":>9}"""
echo repeat("─", 88)
for i in 0..19:
let (d, idx) = distances[i]
let ap = airports[idx]
let b = bearing(𝜑₁, λ₁, ap.latitude, ap.longitude)
echo &"{ap.name:<40}{ap.country:<20}{ap.icao:<11}{d:4.1f}{b.toInt:10}"
- Output:
Airport 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
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.
Pascal
Free Pascal
program Dist_Bearing;
{$IFDEF FPC} {$Mode DELPHI}{$Optimization ON,ALL} {$ENDIF}
{$IFDEF WINDOWS}{$APPTYPE CONSOLE}{$ENDIF}
uses
SysUtils,Math;
const
cDegToRad = pi / 180; cRadToDeg = 180 / pi;
//One nautical mile ( 1" of earth circumfence )
cOneNmInKm = (12742*pi)/360/60;
DiaEarth = 12742/cOneNmInKm;
type
tLatLon = record
lat,lon:double;
sinLat,cosLat:double;
sinLon,cosLon:double;
end;
tDist_Dir = record
distance,bearing:double;
end;
tDst_Bear = record
Koor1,
Koor2 : tLatLon;
Dist_Dir : tDist_Dir;
end;
tmyName = String; //string[63-8] experiment
tmyCountry = String; //string[31]
tmyICAO = String; //string[7]
tSolution = record
Sol_Name : tmyName;
Sol_Country :tmyCountry;
Sol_ICAO : tmyICAO;
Sol_Koor : tLatLon;
Sol_dist_dir:tDist_Dir;
end;
tIdxDist = record
Distance: double;
AirportIdx :Int32;
end;
tMinSols = record
sols : array of tIdxDist;
maxValue: double;
actIdx,
maxidx :Int32;
end;
var
Airports: array of tSolution;
MinSols :tMinSols;
cntInserts : Cardinal;
procedure GetSolData(const OneAirport: String;
var TestSol :tSolution);
var
p1,p2,i1,i2,idx,l : integer;
begin
p1:=1;
idx := 0;
l := length(OneAirport);
repeat
p2 := p1;
i1 := p1;
IF OneAirport[p1] <>'"' then
Begin
repeat
p2 +=1;
until (p2>l) OR (OneAirport[p2]=',');
i2 := p2;
end
else
begin
repeat
p2 +=1;
until (p2>l) OR (OneAirport[p2]='"');
i2 := p2;
i1 +=1;
while (p2<l) do
begin
p2 +=1;
IF (OneAirport[p2]=',')then
break;
end;
end;
idx += 1;
with TestSol do
case idx of
2: Sol_Name := copy(OneAirport,i1,i2-i1);
4: Sol_Country := copy(OneAirport,i1,i2-i1);
6: Sol_ICAO := copy(OneAirport,i1,i2-i1);
7: Begin
With Sol_Koor do begin
lat := StrtoFloat(copy(OneAirport,i1,i2-i1))*cDegToRad;
sincos(lat,sinLat,cosLat);
end;
end;
8: Begin
With Sol_Koor do begin
lon := StrtoFloat(copy(OneAirport,i1,i2-i1))*cDegToRad;
sincos(lon,sinLon,cosLon);
end;
end;
end;
p1:= p2+1;
until (idx>7) OR (p1>l);
end;
function ReadAirports(fileName:String):boolean;
var
TF_Buffer : array[0..1 shl 14 -1] of byte;
AirportsFile: TextFile;
OneAirport : String;
l,cnt : UInt32;
begin
Assign(AirportsFile,fileName);
settextbuf(AirportsFile,TF_Buffer);
{$I-}
reset(AirportsFile);
{$I+}
IF ioResult <> 0 then
Begin
Close(AirportsFile);
EXIT(false);
end;
cnt := 0;
l := 100;
setlength(Airports,l);
while Not(EOF(AirportsFile)) do
Begin
Readln(AirportsFile,OneAirport);
GetSolData(OneAirport,Airports[cnt]);
inc(cnt);
if cnt >= l then
Begin
l := l*13 div 8;
setlength(Airports,l);
end;
end;
setlength(Airports,cnt);
Close(AirportsFile);
exit(true);
end;
procedure Out_MinSol;
var
i: integer;
begin
writeln(' ICAO Distance Bearing Country Airport');
writeln(' ---- -------- ------- -------------- -----------------------------------');
For i := 0 to minSols.actidx do
with AirPorts[minSols.sols[i].AirportIdx] do
writeln(Format(' %4s %8.1f %7.0f %-14s %-35s',
[Sol_ICAO,
Sol_dist_dir.distance*DiaEarth,
Sol_dist_dir.bearing*cRadToDeg,
Sol_Country,Sol_Name]));
writeln;
writeln(cntInserts,' inserts to find them');
end;
procedure Init_MinSol(MaxSolCount:Int32);
begin
setlength(MinSols.sols,MaxSolCount+1);
MinSols.actIdx := -1;
MinSols.maxIdx := MaxSolCount-1;
MinSols.MaxValue := maxdouble;
cntInserts := 0;
end;
procedure Insert_Sol(var sol:tDst_Bear;nrAirport:Int32);
var
dist : double;
idx : Int32;
begin
with MinSols do
begin
idx := actIdx;
dist := sol.Dist_Dir.distance;
if Idx >= maxIdx then
IF MaxValue < dist then
Exit;
if idx >= 0 then
begin
inc(idx);
inc(cntInserts);
while sols[idx-1].Distance >dist do
begin
sols[idx]:= sols[idx-1];
dec(idx);
If idx<=0 then
BREAK;
end;
with sols[idx] do
begin
AirportIdx := nrAirport;
Distance := dist;
end;
//update AirPorts[nrAirport] with right distance/bearing
AirPorts[nrAirport].Sol_dist_dir := sol.Dist_Dir;
if actIdx < maxIdx then
actIdx +=1;
end
else
begin
with sols[0] do
begin
AirportIdx := nrAirport;
Distance := dist;
end;
AirPorts[nrAirport].Sol_dist_dir := sol.Dist_Dir;
MinSols.actIdx := 0;
end;
MaxValue := sols[actIdx].Distance;
end;
end;
procedure Calc_Dist_bear(var Dst_Bear:tDst_Bear);
var
dLonSin,dLonCos,x,y : double;
begin
with Dst_Bear do
Begin
If (Koor1.Lat = Koor2.Lat) AND (Koor1.Lon = Koor2.Lon) then
Begin
Dist_Dir.distance := 0;
Dist_Dir.bearing := 0;
Exit;
end;
sincos(Koor1.lon - Koor2.lon,dLonSin,dLonCos);
//distance
Dist_Dir.distance := arcsin(sqrt(sqr(dLonCos * Koor1.Coslat
- Koor2.Coslat) + sqr(dLonSin* Koor1.Coslat)
+ sqr(Koor1.sinlat - Koor2.sinlat)) / 2);
x := dLonSin*Koor2.Coslat;
y := Koor1.Coslat*Koor2.sinlat - Koor1.sinlat*Koor2.Coslat*dLonCos;
//bearing dLonSin as tmp
dLonSin := ArcTan2(x,y);
if dLonSin < 0 then
dLonSin := -dLonSin
else
dLonSin := 2*pi-dLonSin;
Dist_Dir.bearing := dLonSin;
end;
end;
procedure FindNearest(var testKoors : tDst_Bear;cntAirports,cntNearest:Integer);
var
i : Int32;
begin
Init_MinSol(cntNearest);
For i := 0 to cntAirports-1 do
Begin
testKoors.Koor2 := AirPorts[i].Sol_Koor;
Calc_Dist_bear(testKoors);
Insert_Sol(testKoors,i);
end;
end;
const
rounds = 100;
cntNearest = 20;//128;//8000;
var
T1,T0 : Int64;
testKoors : tDst_Bear;
myKoor : tLatLon;
i,cntAirports : integer;
begin
T0 := Gettickcount64;
IF NOT(ReadAirports('airports.dat')) then
HALT(129);
T1 := Gettickcount64;
Writeln((T1-T0),' ms for reading airports.dat');
cntAirports := length(AirPorts);
with myKoor do
begin
lat := 51.514669*cDegToRad;
lon := 2.198581*cDegToRad;
sincos(lat,sinLat,cosLat);
sincos(lon,sinLon,cosLon);
end;
randomize;
T0 := Gettickcount64;
For i := rounds-2 downto 0 do
Begin
testKoors.Koor1 := AirPorts[random(cntAirports)].Sol_Koor;
FindNearest(testKoors,cntAirports,cntNearest);
end;
testKoors.Koor1 := myKoor;
FindNearest(testKoors,cntAirports,cntNearest);
T1 := Gettickcount64;
Writeln((T1-T0),' ms for searching ',rounds,' times of '
,cntNearest,' nearest out of ',cntAirports,' airports');
writeln(cntInserts,' inserts to find them');
writeln;
FindNearest(testKoors,cntAirports,20);
with myKoor do
writeln(Format('Nearest to latitude %7.5f,longitude %7.5f degrees',
[cRadToDeg*lat,cRadToDeg*lon]));
writeln;
Out_MinSol;
end.
- Output:
7 ms for reading airports.dat 125 ms for searching 100 times of 20 nearest out of 7698 airports 144 inserts to find them Nearest to latitude 51.51467,longitude 2.19858 degrees ICAO Distance Bearing Country Airport ---- -------- ------- -------------- ----------------------------------- EBFN 30.6 146 Belgium Koksijde Air Base EBOS 31.3 127 Belgium Ostend-Bruges International Airport EGMH 33.5 252 United Kingdom Kent International Airport LFAC 34.4 196 France Calais-Dunkerque Airport EBKW 42.5 105 Belgium Westkapelle heliport EGMK 51.6 240 United Kingdom Lympne Airport EBUL 52.8 114 Belgium Ursel Air Base EGMC 56.2 274 United Kingdom Southend Airport LFQT 56.3 163 France Merville-Calonne Airport EBKT 56.4 137 Belgium Wevelgem Airport EHMZ 57.2 90 Netherlands Midden-Zeeland Airport EGMD 58.0 235 United Kingdom Lydd Airport EGUW 58.9 309 United Kingdom RAF Wattisham EGSM 59.3 339 United Kingdom Beccles Airport LFQO 59.6 146 France Lille/Marcq-en-Baroeul Airport EGKH 62.2 250 United Kingdom Lashenden (Headcorn) Airfield LFAT 63.7 200 France Le Touquet-Côte d'Opale Airport EGTO 64.2 262 United Kingdom Rochester Airport LFQQ 66.2 149 France Lille-Lesquin Airport EGMT 68.4 272 United Kingdom Thurrock Airfield real 0m0.134s //test nearest 128 7 ms for reading airports.dat 131 ms for searching 100 times of 128 nearest out of 7698 airports 602 inserts to find them //test nearest of all -> sort all for distance 7 ms for reading airports.dat 1440 ms for searching 100 times of 8000 nearest out of 7698 airports 7697 inserts to find them
Perl
use v5.36;
use utf8;
binmode STDOUT, ':utf8';
use Text::CSV 'csv';
use Math::Trig;
use constant EARTH_RADIUS_IN_NAUTICAL_MILES => 6372.8 / 1.852;
use constant TAU => 2 * 2 * atan2(1, 0);
sub degrees ( $rad ) { $rad / TAU * 360 }
sub radians ( $deg ) { $deg * TAU / 360 }
sub haversine ( $x ) { sin($x / 2)**2 }
sub arc_haver ( $x ) { asin(sqrt($x)) * 2 }
sub great_circle_distance ( $p1, $l1, $p2, $l2 ) {
arc_haver(
haversine($p2 - $p1)
+ haversine($l2 - $l1) * cos($p1) * cos($p2)
);
}
sub great_circle_bearing ( $p1, $l1, $p2, $l2 ) {
atan2( cos($p2) * sin($l2 - $l1),
cos($p1)*sin($p2) - sin($p1) * cos($p2) * cos($l2 - $l1),
);
}
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);
$dist * EARTH_RADIUS_IN_NAUTICAL_MILES, degrees( $theta < 0 ? $theta + TAU : $theta);
}
sub find_nearest_airports ( $latitude, $longitude, $csv_path ) {
my $airports = csv(
in => $csv_path,
headers => [<ID Name City Country IATA ICAO Latitude Longitude>],
);
for my $row (@$airports) {
($$row{'Distance'},$$row{'Bearing'}) = distance_and_bearing( $latitude, $longitude, $$row{'Latitude'}, $$row{'Longitude'} );
}
sort { $a->{'Distance'} <=> $b->{'Distance'} } @$airports;
}
my($lat, $lon, $wanted, $csv) = (51.514669, 2.198581, 20, 'ref/airports.dat');
printf "%7s\t%7s\t%-7s\t%-15s\t%s\n", <Dist Bear ICAO Country Name>;
for my $airport (find_nearest_airports($lat, $lon, $csv)) {
printf "%7.1f\t %03d\t%-7s\t%-15s\t%s\n", map { $airport->{$_} } <Distance Bearing ICAO Country Name>;
last unless --$wanted
}
- 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 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 72.5 313 EGXH United Kingdom RAF Honington
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
Ruby
require 'open-uri'
require 'csv'
include Math
RADIUS = 6372.8 # rough radius of the Earth, in kilometers
def spherical_distance(start_coords, end_coords)
lat1, long1 = deg2rad(*start_coords)
lat2, long2 = deg2rad(*end_coords)
2 * RADIUS * asin(sqrt(sin((lat2-lat1)/2)**2 + cos(lat1) * cos(lat2) * sin((long2 - long1)/2)**2))
end
def bearing(start_coords, end_coords)
lat1, long1 = deg2rad(*start_coords)
lat2, long2 = deg2rad(*end_coords)
dlon = long2 - long1
atan2(sin(dlon) * cos(lat2), cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dlon))
end
def deg2rad(lat, long)
[lat * PI / 180, long * PI / 180]
end
uri = "https://raw.githubusercontent.com/jpatokal/openflights/master/data/airports.dat"
headers = %i(airportID
name
city
country
iata
icao
latitude
longitude
altitude
timezone
dst
tzOlson
type
source)
data = CSV.parse(URI.open(uri), headers: headers, converters: :numeric)
position = [51.514669, 2.198581]
data.each{|r| r[:dist] = (spherical_distance(position, [r[:latitude], r[:longitude]])/1.852).round(1)}
closest = data.min_by(20){|row| row[:dist] }
closest.each do |r|
bearing = (bearing(position,[r[:latitude], r[:longitude]])*180/PI).round % 360
puts "%-40s %-25s %-6s %12.1f %15.0f" % (r.values_at(:name, :country, :ICAO, :dist) << bearing)
end
- Output:
Koksijde Air Base Belgium 30.7 146 Ostend-Bruges International Airport Belgium 31.3 127 Kent International Airport United Kingdom 33.6 252 Calais-Dunkerque Airport France 34.4 196 Westkapelle heliport Belgium 42.6 105 Lympne Airport United Kingdom 51.6 240 Ursel Air Base Belgium 52.8 114 Southend Airport United Kingdom 56.2 274 Merville-Calonne Airport France 56.4 163 Wevelgem Airport Belgium 56.5 137 Midden-Zeeland Airport Netherlands 57.3 90 Lydd Airport United Kingdom 58.0 235 RAF Wattisham United Kingdom 59.0 309 Beccles Airport United Kingdom 59.3 339 Lille/Marcq-en-Baroeul Airport France 59.7 146 Lashenden (Headcorn) Airfield United Kingdom 62.2 250 Le Touquet-Côte d'Opale Airport France 63.7 200 Rochester Airport United Kingdom 64.2 262 Lille-Lesquin Airport France 66.3 149 Thurrock Airfield United Kingdom 68.4 272
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
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)